from ROOT import *
from array import array
import re
import getopt

try:
    opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:", [])
except getopt.GetoptError:
    # print help information and exit:
    print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents'
    sys.exit()
for o, a in opts:
    if o in ("-f"):
        fileName = a
    if o in ("-o"):
        outputFileName = a
    if o in ("-n"):
        maxEvents = int(a)
    if o in ("-t"):
        sampleType = str(a)

fileNameGeo = fileName.replace('ship', 'geofile_full')
gROOT.ProcessLine(".x mystyle.C")
pdg = TDatabasePDG.Instance()

sampleType = 'nuBg'

oldGeo = False
## are the weights to be switched on?
ON = True
nu = True

#where = "../../data/74x_75x/"
#where = "/afs/cern.ch/work/b/bstora/data/"
#
#fileNameGeo = "%sgeofile_full.10.0.Genie-TGeant4.root"%where
#
#if nu: 
#    fileName = "%sship.10.0.Genie-TGeant4_rec_neutrino_newReco_nu_77_79.root"%where
#    outputFileName = "%sfromAnalysis/ntupleStudy_highFlux_nu_77_77_higherThr.root"%where
#else:
#    fileName = "%sship.10.0.Genie-TGeant4_rec_neutrino_newReco_antinu_76_78.root"%where
#    outputFileName = "%sfromAnalysis/ntupleStudy_highFlux_antinu_76_78_higherThr.root"%where


if not(pdg.GetParticle(1000010020)):
    pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020)
    pdg.AddAntiParticle("AntiDeuteron", - 1000010020)
    
if not(pdg.GetParticle(1000010030)):
    pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030);
    pdg.AddAntiParticle("AntiTriton", - 1000010030);

if not(pdg.GetParticle(1000020040) ):
    #pdg.AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,khShGev/(12.33*kYear2Sec),6,"Ion",ionCode);
    pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040);
    pdg.AddAntiParticle("AntiAlpha", - 1000020040);
        
if not(pdg.GetParticle(1000020030) ):
    pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030);
    pdg.AddAntiParticle("AntiHE3", -1000020030);

if not(pdg.GetParticle(1000030070) ):
    print "TO BE CHECKED the data insert for Li7Nucleus"
    pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070);
    
if not(pdg.GetParticle(1000060120) ):
    print "ERROR: random values insert for C12Nucleus"
    pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120);
        
if not(pdg.GetParticle(1000030060) ):
    print "ERROR: random values insert for Li6Nucleus"
    pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060);       

if not(pdg.GetParticle(1000070140) ):
    print "ERROR: random values insert for for N14"
    pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140);       

if not(pdg.GetParticle(1000050100) ):
    print "TO BE CHECKED the data insert for B10"
    pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100);

if not(pdg.GetParticle(1000020060) ):
    print "TO BE CHECKED the data insert for He6"
    pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060);

if not(pdg.GetParticle(1000040080) ):
    print "TO BE CHECKED the data insert for Be8"
    pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080);

if not(pdg.GetParticle(1000030080) ):
    print "TO BE CHECKED the data insert for Li8"
    pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080);    
 
if not(pdg.GetParticle(1000040170) ):
    print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!"
    pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170);    

if not(pdg.GetParticle(1000040100) ):
    print "TO BE CHECKED the data insert for Be10"
    pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100);    

if not(pdg.GetParticle(1000040070) ):
    print "TO BE CHECKED the data insert for Be7"
    pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070);    

if not(pdg.GetParticle(1000230470) ):
    print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!"
    pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470);    

if not(pdg.GetParticle(1000080170) ):
    print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!"
    pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170);    

if not(pdg.GetParticle(1000240500) ):
    print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!"
    pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500);    

if not(pdg.GetParticle(1000210450) ):
    print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!"
    pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450);    

if not(pdg.GetParticle(1000040090) ):
    print "TO BE CHECKED the data insert for Be9"
    pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090);    

if not(pdg.GetParticle(1000080160) ):
    print "TO BE CHECKED the data insert for O16"
    pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160);    

if not(pdg.GetParticle(1000220460) ):
    print "TO BE CHECKED the data insert for Ar40"
    pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460);   

if not(pdg.GetParticle(1000050110) ):
    print "TO BE CHECKED the data insert for B11"
    pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110);   



def PrintEventPart(particles,pdg):
    print "Particles in the event:"
    for (pid,part) in enumerate(particles):
        print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId()
    print

def PrintRPChits(rpc,pdg):
    print "Hits in:"
    for (ix,RPCpt) in enumerate(rpc):
        #if pdg.GetParticle(RPCpt.PdgCode())<10000:
        #    tmpName = pdg.GetParticle(RPCpt.PdgCode()).GetName()
        #else:
        tmpName = pdg.GetParticle(RPCpt.PdgCode())
        #if RPCpt.GetZ()>=-2484.0 and RPCpt.GetZ()<=-2482.0:
        print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName()
    print

KsID = 310
KLID = 130
NID = 2112

# weight computation moved to offlineForBarbara.py
"""
ff = TFile("histoForWeights.root")
h_GioHans = ff.Get("h_Gio")
# t.Draw("interactionElement", "(isPrimary==1)*weight")

def calcWeight(NC, E, w, entries, nuName, ON=True):
    if not ON:
        return 1

    if "bar" in nuName:
        reduction = 0.5
        #flux = 1.88e+11 * 2.e+20 / 5.e+13
        flux = 6.98e+11 * 2.e+20 / 5.e+13
    else: 
        reduction = 1.
        #flux = 2.88e+11 * 2.e+20/ 5.e+13
        flux = 1.09e+12 * 2.e+20/ 5.e+13

        #reduction = 1.
        #flux = 2.e+18    
        #entries = 99999+99999.
    crossSec = 6.7e-39*E * reduction

    
        
        #flux = 2.e+18
    NA = 6.022e+23
    
    binN = h_GioHans.GetXaxis().FindBin(E)    
    return crossSec * flux * h_GioHans.GetBinContent(binN) * w * NA / entries
"""

def addVect(t,name,vectType):
    vect =  vector(vectType)()
    t.Branch( name, vect )
    return t, vect

def addVar(t,name,varType):
    var = array(varType,[-999])
    t.Branch(name,var,name+"/"+varType.upper())
    return t, var

def putToZero(var):
    var[0] = 0
    
def getPartName(partId):
    return pdg.GetParticle(partId).GetName()

def lookingForDecay(particles):
    decayMotherList = []
    
    for p in particles:
        mID = p.GetMotherId()
        if mID>=0 and not (mID in decayMotherList): 
            decayMotherList.append(mID)
            decayPartIndex.push_back(mID)
            decayPartID.push_back(particles[mID].GetPdgCode())
            decayPartStrID.push_back(pdg.GetParticle(particles[mID].GetPdgCode()).GetName())
            decayPos_x.push_back(p.GetStartX())
            decayPos_y.push_back(p.GetStartY())
            decayPos_z.push_back(p.GetStartZ())
            decayMother.push_back(particles[mID].GetMotherId())
            decayP.push_back(p.GetP())

def wasFired(indexKids, detPoints, detPos, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):

    def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId):
        charges = []
        trackids = []
        if hitCharges is None:
            assert(hitTrackId is None)
            
        for pos in detPos:
            foundStat = False
            foundStat_eff = False
            
            for hit in detPoints:
                if (indexKids is None) or (hit.GetTrackID() in indexKids):
                # check if it is in one of the considered active layer
                    if pos[0]<=hit.GetZ()<=pos[1]:
                        Eloss += hit.GetEnergyLoss()
                        hitID = hit.GetTrackID()
                        if not hitID in trackids and not hitCharges is None:
                            if  hit.PdgCode()>100000:
                                charges.append(9)
                            else:
                                charges.append(int(pdg.GetParticle(hit.PdgCode()).Charge()))
                            trackids.append(hitID)
                            hitCharges.push_back(charges[-1])
                            hitTrackId.push_back(trackids[-1])
                        if pointVects is not None:
                            pointVects[0].push_back(hit.GetX())
                            pointVects[1].push_back(hit.GetY())
                            pointVects[2].push_back(hit.GetZ())
                            
                        flag = True
                        #nstat+=1
                        foundStat = True
                        eff_val = gRandom.Uniform(0.,1.)
                        if eff_val<0.9:
                            flag_eff = flag_eff or True
                            foundStat_eff = True
                            #nstat_eff+=1
                        else:
                            flag_eff = flag_eff or False
            if foundStat:
                nstat+=1
            if foundStat_eff:
                nstat_eff+=1
                
        #charges.sort()
        particles = 0
        #for i in xrange(int(len(charges)/2.)):
        #    if charges[0+i]*charges[-1-i]<0:
        #        particles+=1
        #    else:
        #        break
        flag_Ethr = Eloss>=Ethr
        return flag, flag_eff, nstat, nstat_eff, Eloss, flag_Ethr,particles


    # Now in partKidTrackID I should have the trackID connected to my charged particles
    
    flag_eff = False
    flag = False
    nstat = 0
    nstat_eff = 0
    Eloss = 0
    
    flag,flag_eff,nstat,nstat_eff,Eloss,flag_Ethr,particles = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId)
    
    if flag==False and checkOn:
        print "To be study event %s"%entry        
        
    return flag, flag_eff, nstat, nstat_eff, flag_Ethr, particles

def convertion(tmpName):
    if "Rib" in tmpName: tmpName = "Rib"
    elif "LiSc" in tmpName: tmpName= "LiSc"
    elif "Tr" in tmpName: tmpName = "Tr"
    elif "Startplate" in tmpName: tmpName = "Startplate"
    elif re.search("T\d+O", tmpName): tmpName = "TO"
    elif re.search("T\d+I", tmpName): tmpName = "TI"
    elif "Endplate" in tmpName: tmpName = "Endplate"
        #if "VetoTimeDet" in tmpName: tmpName = "upstreamVeto"
    elif "volCoil" in tmpName: tmpName = "volCoil"
    elif "volFeYoke" in tmpName: tmpName = "volFeYoke"
    elif "gas" in tmpName: tmpName = "gas"
    elif "wire" in tmpName: tmpName == "wire"
    elif "VetoTimeDet" in tmpName: tmpName = "upstreamVeto"
    elif "Veto" in tmpName: tmpName = "strawVeto"
    return tmpName

def wasFired_node(indexKids, detPoints, detVolNames, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):

    def lookingForHits(detPoints,flag,hitCharges, hitTrackId):
        if hitCharges is None:
            assert(hitTrackId is None)
        
        for hit in detPoints:
            if (indexKids is None) or (hit.GetTrackID() in indexKids):
                #print hit.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName(), detVolNames
                             
                tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName()
                tmpName = convertion(tmpName)
                if hit.GetZ()>8000:
                    continue
                #print tmpName, detVolNames
                if tmpName in detVolNames:
                    ## trick to remove the case of the tracking system the hits in the gas or straw from the straw-veto:
                    if 'trackingSystem' in detVolNames and hit.GetZ()<0:
                        continue
                    if "strawVetoSystem" in detVolNames and hit.GetZ()>0:
                        continue
                    #if "upstreamVeto" in detVolNames and not(hit.GetZ()>=-2484.0 and hit.GetZ()<=-2482.0):
                        #upstreamVetoPos
                    if "upstreamVeto" in detVolNames and not(hit.GetZ()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]):    
                        continue
                #print RPChit.GetTrackID(), t_ID
                # check if it is in one of the considered active layer
                    if True: #pos[0]<=hit.GetZ()<=pos[1]:
                        flag = True
                
        return flag


    # Now in partKidTrackID I should have the trackID connected to my charged particles
    
    flag = False
    
    flag = lookingForHits(detPoints,flag,hitCharges, hitTrackId)
    
    if flag==False and checkOn:
        print "To be study event %s"%entry        
        
    return flag



f = TFile(fileName)
t = f.Get("cbmsim")

entries = t.GetEntries()
from lookAtGeo import *


r = loadGeometry(fileNameGeo)
fGeo = r['fGeo']

# stuff by elena
from offlineForBarbara import *

myNodes_name = ["VetoTimeDet_1"]
myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
myNodes_name += ["Veto_5"]

myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None)

lastPassive_nodeName = ["volIron_23"]
OPERA_nodeName = ["volIron", "volRPC", "volHPT","volArm2MS","volFeYokes","volCoil","volFeYoke"]#["volIron_%s"%i for i in xrange(12,24)]+["volRpc_%s"%i for i in xrange(11,22)]+["volArm2MS_1"] 
entrance_nodeName = ["T1Lid"]#['T1lid_1']
volumeIn_nodeName = ["TI"]#['T%sI'%i for i in [1,2,3,5]]
volumeOut_nodeName = ["TO"]#['T%sO'%i for i in [1,2,3,5]]

OPERA_volNames = ["volIron","volRpc","volHPT"]
strawVeto_volNames = ["Veto"]

top=fGeo.GetTopVolume()

Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']]
upstreamVeto = [myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]

trackStationsPos = [[myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ'],myGeoEl["Tr%s_%s"%(i,i)]['z']+myGeoEl["Tr%s_%s"%(i,i)]['dimZ']] for i in xrange(1,5)]
strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]]
vetoWall = [[myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+0.001,myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ']]]#myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+6000.]]
upstreamVetoPos = [[myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]]
RPCstationsPos = [[-3500,myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ']-0.5]]


trackStationsPos_node = ["Tr", "gas", "straw", 'trackingSystem']
## should be that this does not work since it could be that hits are also assigned to gas or straw
strawVetoPos_node = ["strawVeto", "gas", "straw", "strawVetoSystem"]
vetoWall_node = ["LiSc","cave","Rib","TI","TO","Tr","strawVeto"]
upstreamVetoPos_node = ["upstreamVeto","cave"]
RPCstationsPos_node = ["volRpc","volMagneticSpectrometer","volIron","volHPT"]

print "OPERApos: ",RPCstationsPos
print "trackingPos: ",trackStationsPos
print "strawVetoPos: ",strawVetoPos
print "scintPos: ",vetoWall
print "upstreamVetoPos: ",upstreamVetoPos


### Ntuple elements:
nf = TFile(outputFileName,"RECREATE")

nt = TTree("t","t")
nt, event = addVar(nt, 'event','i')
nt, nNu = addVar(nt, 'nNu', 'i')
nt, isPrimary = addVar(nt, 'isPrimary','i')
nt, startZ_nu = addVar(nt, 'startZ_nu', 'f')
nt, startY_nu = addVar(nt, 'startY_nu', 'f')
nt, startX_nu = addVar(nt, 'startX_nu', 'f')
nt, nuE = addVar(nt, 'nuE', 'f')
nt, weight = addVar(nt, 'weight', 'f')
#nt, hansW = addVar(nt, 'hansW', 'f')

## 0: last passive material opera-mu-system
## 1: two windows around liquid scintillator
## : along vaccum tank
## 3: along all OPERA
## 4: between the two entrance windows
## 5: along vacuum tank outer window
## 6: along vacuum tank inner window
## 999: wrong OPERA place ## needed until we have the new production
## -1: anything else, but what???? #it should never been present

nt, interactionElement = addVar(nt, 'interactionElement','i')
### From Neutrino
nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i')
nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int')
nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string')

nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i')
nt, idChargedPart_fromNu =  addVect(nt, 'idChargedPart_fromNu', 'int')
nt, idStrChargedPart_fromNu =  addVect(nt, 'idStrChargedPart_fromNu', 'string')

nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i')
nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int')
nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string')

'''
nt, rpcPoint_x = addVect(nt, 'rpcPoint_x', 'float')
nt, rpcPoint_y = addVect(nt, 'rpcPoint_y', 'float')
nt, rpcPoint_z = addVect(nt, 'rpcPoint_z', 'float')

nt, strawPoint_x = addVect(nt, 'strawPoint_x', 'float')
nt, strawPoint_y = addVect(nt, 'strawPoint_y', 'float')
nt, strawPoint_z = addVect(nt, 'strawPoint_z', 'float')

nt, strawVetoPoint_x = addVect(nt, 'strawVetoPoint_x', 'float')
nt, strawVetoPoint_y = addVect(nt, 'strawVetoPoint_y', 'float')
nt, strawVetoPoint_z = addVect(nt, 'strawVetoPoint_z', 'float')


nt, scintWPoint_x = addVect(nt, 'upstreamVetoPoint_x', 'float')
nt, scintWPoint_y = addVect(nt, 'upstreamVetoPoint_y', 'float')
nt, scintWPoint_z = addVect(nt, 'upstreamVetoPoint_z', 'float')
'''
### RPC ###
# was at least an RPC fired by something coming from the nu interaction? 
nt, RPC = addVar(nt,'RPC', 'i' )
# accounting for an eff. of each station of 90%
nt, RPC_eff = addVar(nt,'RPC_eff', 'i' )

#nt, nRPCstations = addVar(nt, 'nRPCstations', 'i')
#nt, nRPCstations_eff = addVar(nt, 'nRPCstations_eff', 'i')

## In case something (no matter if it comes from the nu-interaction kids) fired the RPC
#nt, nRPCstationsAny = addVar(nt, 'nRPCstationsAny', 'i')
nt, RPCany = addVar(nt,'RPCany', 'i' )
#nt, nRPCstationsAny_eff = addVar(nt, 'nRPCstationsAny_eff', 'i')
nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' )
############

### scintW ###
# was at least an scintW fired by something coming from the nu interaction? 
nt, scintW = addVar(nt,'upstreamVeto', 'i' )
nt, scintWany = addVar(nt,'upstreamVetoAny', 'i' )

# accounting for an eff. of each station of 90%
nt, scintW_eff = addVar(nt,'upstreamVeto_eff', 'i' )
nt, scintWany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' )

# was at least an scintW fired by something coming from the nu interaction? 
nt, scintVeto = addVar(nt,'scintVeto', 'i' )
# accounting for an eff. of each station of 90%
nt, scintVeto_eff = addVar(nt,'scintVeto_eff', 'i' )


## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls
#nt, nScintVetoStationsAny = addVar(nt, 'nScintVetoStationsAny', 'i')
nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' )
nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' )

## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls
nt, scintVetoAny = addVar(nt,'upstreamVetoAny', 'i' )
nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' )

nt, hitCharges = addVect(nt, "hitCharges", 'int')
nt, hitTrackId = addVect(nt, "hitTrackId", 'int')

###############

### strawVeto ###
# was at least an strawVeto fired by something coming from the nu interaction? 
nt, strawVeto = addVar(nt,'strawVeto', 'i' )
# accounting for an eff. of each station of 90%
nt, strawVeto_eff = addVar(nt,'strawVeto_eff', 'i' )


## In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto
#nt, nStrawVetoStationsAny = addVar(nt, 'nStrawVetoStationsAny', 'i')
nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' )
nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' )
###############


### TrackSyst ###
## only the case of anything fired it is considered, obviously
nt, TrackSyst = addVar(nt, "TrackSyst", 'i')
nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i')

##################


nt, decayPartIndex = addVect(nt, "decayPartIndex", "float")
nt, decayPartID = addVect(nt, "decayPartID", "float")
nt, decayPartStrID = addVect(nt, "decayPartStrID", "string")
nt, decayPos_x = addVect(nt, "decayPos_x", "float")
nt, decayPos_y = addVect(nt, "decayPos_y", "float")
nt, decayPos_z = addVect(nt, "decayPos_z", "float")
nt, decayMother = addVect(nt, "decayMother", "float")
nt, decayP = addVect(nt,"decayP","float")

nt, strawVeto_Ethr = addVar(nt,"strawVeto_Ethr","i")
nt, scintVeto_Ethr = addVar(nt,"scintVeto_Ethr","i")
nt, scintW_Ethr = addVar(nt,"upstreamVeto_Ethr","i")
nt, RPC_Ethr = addVar(nt,"RPC_Ethr","i")
nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i")
nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i")
nt, scintWAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i")
nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i")
nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i")

nt, NC = addVar(nt, "NC", "i")

#nt, nParticles = addVar(nt, "nParticles", "i")

nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string")
#nt, nuInteractionNode_sel = addVar(nt, "nuInteractionNode_sel","i")
nt, nuInteractionNodeSimpl = addVect(nt, "nuInteractionNodeSimpl", "string")

nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","i")

nt, TrackSyst_node = addVar(nt, "TrackSyst_node", "i")
nt, strawVetoAny_node = addVar(nt, "strawVetoAny_node", "i")
nt, scintWany_node = addVar(nt, "scintWany_node", "i")
nt, RPCany_node = addVar(nt, "RPCany_node", "i")
nt, scintVetoAny_node = addVar(nt, "scintVetoAny_node", "i")

nt, recoed = addVar(nt, "recoed","i")
nt, nRecoed = addVar(nt, "nRecoed","i")

dictNodeNames = {'volIron':0, 'cave':1, 'LiSc':2, 'Startplate':3, 'TI':4, 'rockD':5, 'Endplate':6, 'Rib':7, 'volFeYoke':8, 'volHPT':9, 'TO':10, 'volRpc':11, 'volCoil':12, 'T1Lid':13, 'straw':14, 'strawVeto':15, 'volBase':16, 'Tr':17, 'gas':18, 'wire':19, 'others':20}

entries = 10000
for entry in xrange(entries):
    if not (entry%1000):
        print "%s / %s"%(entry,entries)
        
    t.GetEntry(entry)

    event[0] = entry
    particles = t.MCTrack
    rpc = t.ShipRpcPoint
    scint = t.vetoPoint
    strawPoints = t.strawtubesPoint
    vetoScint = t.vetoPoint
    recoParts = t.Particles

    ## initialization
    putToZero(nNu)
    putToZero(nPart_fromNu)
    putToZero(nChargedPart_fromNu)
    putToZero(nNeutrPart_fromNu)
    
    '''putToZero(nNeutrPart_withKs)
    putToZero(nChargedPart_withKs)
    putToZero(nPart_withKs)
    
    putToZero(nNeutrPart_withKL)
    putToZero(nChargedPart_withKL)
    putToZero(nPart_withKL)
    
    putToZero(nNeutrPart_withN)
    putToZero(nChargedPart_withN)
    putToZero(nPart_withN)
    putToZero(hasKs)
    putToZero(hasKL)
    putToZero(hasN)
    '''
    hitCharges.clear() 
    hitTrackId.clear()
    
    idPart_fromNu.clear()
    idChargedPart_fromNu.clear()
    idNeutrPart_fromNu.clear()
    idStrPart_fromNu.clear()
    idStrChargedPart_fromNu.clear()
    idStrNeutrPart_fromNu.clear()

    '''
    strawVetoPoint_x.clear()
    strawVetoPoint_y.clear()
    strawVetoPoint_z.clear()

    strawPoint_x.clear()
    strawPoint_y.clear()
    strawPoint_z.clear()

    rpcPoint_x.clear()
    rpcPoint_y.clear()
    rpcPoint_z.clear()

    scintWPoint_x.clear()
    scintWPoint_y.clear()
    scintWPoint_z.clear()
    '''
    decayPartIndex.clear()
    decayPartID.clear()
    decayPartStrID.clear()
    decayPos_x.clear()
    decayPos_y.clear()
    decayPos_z.clear()
    decayMother.clear()
    decayP.clear()

    nuInteractionNode.clear()
    
    primaryDone=False
    isPrimary[0] = int(False)

    lookingForDecay(particles)
    nRecoed[0] = 0
    recoed[0] = 0
    #for (ip,part) in enumerate(particles):
    
    if sampleType=="nuBg" or sampleType == "cosmics":
        starPartRange = [0] 
    elif sampleType=="sig":
        starPartRange = [1] 

    for ip in startPartRange:

        TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, hitCharges, hitTrackId, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] )
        if TrackSyst[0]==0:
            skipEvent=True
            continue
        skipEvent=False
        
        part = particles[ip]

        pdgPart = pdg.GetParticle(part.GetPdgCode())
        if sampleType=="nuBg":
            assert("nu" in pdgPart.GetName())        

        pushOfflineByEvent(t, sampleType)
        ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother
        #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena
        if True:
            ## Starting the counter of how many particles were produced by the interaction of this specific nu
            #fromThisNu.push_back(0)
            for recoP in recoParts:
                nRecoed[0] += 1
                pushOfflineByParticle(t, recoP)

            if nRecoed[0]>0:
                recoed[0]=1
                #else:
                #skipEvent=True
                #continue
            
            nNu[0]+=1
            if part.GetMotherId()==-1:
                assert(primaryDone==False)
                primaryDone=True
                isPrimary[0] = int(True) 
                assert(len(idStrPart_fromNu)==0)
                assert(len(idStrNeutrPart_fromNu)==0)
            else: 
                continue


            tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName()
            nuInteractionNode.push_back(tmpName)
            
            tmpName = convertion(tmpName)
            
            nuInteractionNodeSimpl.push_back(tmpName)
            if not tmpName in dictNodeNames: tmpName = "others"
                
            nuIntNumSimpl[0] = dictNodeNames[tmpName]
            
            startZ_nu[0] = part.GetStartZ()
            startY_nu[0] = part.GetStartY()
            startX_nu[0] = part.GetStartX()
                
            nuE[0] = part.GetEnergy()
            # Re-initialization for each interacting neutrino
            ## ---> <--- ##
            
            ## Looking for particles produced by the neutrino interaction
            interacted = False
            partKidsId = []
            #indexPartKids = []
            
            #maxPcharged[0] = 0
            #maxPneutr[0] = 0
            NC[0] = 0
            for ip2 in xrange(0,len(particles)):
                part2 = particles[ip2]
                ## exit if we have reached the empty part of the array
                if not (type(part2)==type(ShipMCTrack())):
                    break
                if part2.GetMotherId()==ip:
                    interacted = True 
                    part2Id = part2.GetPdgCode()
                    partKidsId.append(part2Id)
                    #indexPartKids.append(ip2)

                    if not pdg.GetParticle(part2.GetPdgCode()) == None and  "nu" in pdg.GetParticle(part2.GetPdgCode()).GetName() and part2.GetMotherId()==0:
                        NC[0] = 1
                    #                    kid_decay_x.push_back(part2.GetX())
                    #kid_decay_y.push_back(part2.GetY())
                    #kid_decay_z.push_back(part2.GetZ())
                
                    
                    tmpP = part2.GetP()
                    #print "pID",part2Id, pdg.GetParticle(part2Id).GetName()
                    
                    nPart_fromNu[0]+=1
                    idPart_fromNu.push_back(part2Id)
                    if pdg.GetParticle(part2.GetPdgCode()) == None:
                        part2Name = str(part2Id)
                    else:
                        part2Name = getPartName(part2Id)
                    idStrPart_fromNu.push_back(part2Name)
                    
                    ## Counting how many particles produced by the nu-interaction are charged or neutral
                    if pdg.GetParticle(part2.GetPdgCode()) == None or int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0):
                        nNeutrPart_fromNu[0]+=1
                        idNeutrPart_fromNu.push_back(part2Id)
                        idStrNeutrPart_fromNu.push_back(part2Name)
                        #if fabs(maxPcharged[0])<fabs(tmpP):
                        #    maxPcharged[0]=tmpP
                    else:
                        nChargedPart_fromNu[0]+=1
                        idChargedPart_fromNu.push_back(part2Id)
                        idStrChargedPart_fromNu.push_back(part2Name)
                        #if fabs(maxPneutr[0])<fabs(tmpP):
                        #    maxPneutr[0]=tmpP

                    ## How many are produced by the interaction in the last passive element of the "opera-mu system"?
                    part2Z = part2.GetStartZ()
                    part2X = part2.GetStartX()
                    part2Y = part2.GetStartY()
                    #nuKid_z.push_back(part2Z)
                    somewhere = False
                    nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName()
                    #if part2Z>lastPassiveEl[0] and part2Z<lastPassiveEl[1]:
                    #print "----> ",nodeName, lastPassive_nodeName
                    
                    if  nodeName in lastPassive_nodeName:
                        somewhere = True
                        interactionElement[0] = 0 
                        intElName = 'OPERA'
                        #RpcPassiveFlag = True
                        #lastPassive[0] = int(True)
                        
                    ## To know if it was in the full OPERA-system (excluded last passive)
                    #if part2Z>OPERA[0] and part2Z<OPERA[1] and somewhere==False:
                    if tmpName in OPERA_nodeName and somewhere==False:
                        somewhere = True
                        interactionElement[0] = 3
                        intElName = 'OPERA'
                        ## To account for things that are both in the correct and wrong OPERA place
                        #if part2Z>OPERA_wrong[0] and part2Z<OPERA_wrong[1]:
                        #   interactionElement.push_back(999)

                    

                    ## To know if it was between the two windows
                    #if part2Z>entranceWindows[0][1] and part2Z<entranceWindows[1][0]:
                    if tmpName in entrance_nodeName:
                        assert(somewhere==False)
                        somewhere = True
                        interactionElement[0] = 1

                    if tmpName in volumeIn_nodeName and somewhere==False:
                        interactionElement[0] = 5
                        somewhere = True 
                    elif nodeName in volumeOut_nodeName and somewhere==False:
                        interactionElement[0] = 6
                        somewhere = True 

                    if "Rib" in tmpName:
                        interactionElement[0] = 7
                        
                    if somewhere==False:
                        interactionElement[0] = -1

    if skipEvent:
        continue
    #strawVeto[0],strawVeto_eff[0],dummy,dummy,strawVeto_Ethr[0],dummy = wasFired(indexPartKids, strawPoints, strawVetoPos, None, None, checkOn=False )
    strawVetoAny[0],strawVetoAny_eff[0],dummy,dummy,strawVetoAny_Ethr[0],dummy  = wasFired(None, strawPoints, strawVetoPos,None,None, checkOn=False, pointVects=None)#[strawVetoPoint_x, strawVetoPoint_y, strawVetoPoint_z])

    #scintW[0], scintW_eff[0], dummy,dummy,scintW_Ethr[0],dummy = wasFired(indexPartKids, vetoScint, upstreamVetoPos, None, None, checkOn=False)
    scintWany[0], scintWany_eff[0], dummy,dummy, scintWAny_Ethr[0],dummy = wasFired(None, vetoScint, upstreamVetoPos, None, None, checkOn=False, pointVects=None)#[scintWPoint_x, scintWPoint_y, scintWPoint_z])
    #RPC[0], RPC_eff[0], dummy,dummy, RPC_Ethr[0],dummy = wasFired(indexPartKids, rpc, RPCstationsPos, None, None, checkOn=False)
    RPCany[0], RPCany_eff[0],dummy,dummy, RPCany_Ethr[0],dummy = wasFired(None, rpc, RPCstationsPos, None, None, checkOn=False, pointVects=None)#[rpcPoint_x, rpcPoint_y, rpcPoint_z])
    #scintVeto[0], scintVeto_eff[0], dummy,dummy,scintVeto_Ethr[0],dummy = wasFired(indexPartKids, vetoScint, vetoWall, None, None, checkOn=False, Ethr=0.045)
    scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=0.045)
  
    assert(len(idStrPart_fromNu)==len(idPart_fromNu))
    assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu))
    assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu))
    assert(len(idStrPart_fromNu)==nChargedPart_fromNu[0]+nNeutrPart_fromNu[0])
    assert(len(idStrChargedPart_fromNu)==nChargedPart_fromNu[0])
    assert(len(idStrNeutrPart_fromNu)==nNeutrPart_fromNu[0])
    if not primaryDone:
        print "It looks like there is not primary nu interacting"
        PrintEventPart(particles,pdg)
        assert(False)
  
    weight[0] = findWeight(sampleType, NC[0], nuE[0], part, entries, pdgPart.GetName(), ON)#calcWeight(NC[0], nuE[0], part.GetWeight(), entries, pdgPart.GetName(), ON)
    #hansW[0] = part.GetWeight()
    nt.Fill()
        

        
nt.Write()
nf.Save()
nf.Close()
print "Wrote file %s"%outputFileName
'''f2 = TFile("prova.root")
t2 = f2.Get("t")
t2.Print()

t2.Draw("nuE", "nNu==1 && isPrimary==1")
#t2.Draw("startZ_nu", "nNu==1 && isPrimary==1 && startZ_nu>=-2501.8 && startZ_nu<=-2478.8")
'''


'''### Looking for the first RPC station that can be fired (i.e. after the interaction point)   
                leftRPC = None#tmp_firstRPCStation = None
                rightRPC = None
                for (indexRPC,RPCs) in enumerate(RPCstationsPos):
                    if part2Z<RPCs[1]:
                        #tmp_firstRPCStation = indexRPC
                        leftRPC = indexRPC-1
                        rightRPC = indexRPC
                        break
                    #if tmp_firstRPCStation is None:
                if leftRPC is None:
                    print "ERROR: need to debug event %s"%entry
                    print "part2z = ",part2Z
                    print "RPCstationPos = ",RPCstationsPos
                    sys.exit(0)
                    #if tmp_firstRPCStation is not 0:
                    # print "Did you change and you are not looking anymore at the last passive element? If not it should be a bug in the code, if yes remove this part (used only as a safety check)"
                    #sys.exit(0)
                print "###", leftRPC, rightRPC, RPCstationsPos[leftRPC], RPCstationsPos[rightRPC], part2Z
                # I used a list in case we would like to include more than the two RPC layers around the passive element where we had the interaction
                potentiallyFiredRPCindex = [leftRPC,rightRPC]

                ## checking if I have one of the two active stations around the interaction point with hits
                RPCflag = False
                RPCflag_eff = False
                for RPCpt in rpc:
                    print RPCpt.GetZ(), pdg.GetParticle(RPCpt.PdgCode()).GetName(), RPCpt.GetTrackID()
                assert(False)
                for RPCpt in rpc:
                    for potRPCindex in potentiallyFiredRPCindex:
                        print RPCpt.GetZ(), RPCstationsPos[potRPCindex]
                        if RPCpt.GetZ()<RPCstationsPos[potRPCindex][1] and RPCpt.GetZ()>RPCstationsPos[potRPCindex][0]:
                            RPCflag = True
                            
                            ## For each station after the interaction point I simulate a 90% efficiency of the detector
                            ## only if one station is still fired I put a positive flag
                            eff_flag = gRandom.Uniform(0.,1.)
                            if eff_flag<0.9:
                                RPCflag_eff = True
                        print RPCflag, RPCflag_eff

                for (pid,part) in enumerate(particles):
                    print pid, pdg.GetParticle(part.GetName(), part.GetMotherId()
                assert(False)
'''     
