- import sys, re, getopt
- from array import array
- from ROOT import *
-
- fileName, fileNameGeo, outputFileName, sampleType, jobID = None, None, None, None, None
- maxevents = 0
- verbose = True
- # Energy deposit threshold for the liquid scintillator
- threshold = 0.045
-
- # Parse commandline inputs
- try:
- opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:j:", ['Ethr='])
- except getopt.GetoptError:
- # print help information and exit:
- print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents -v verbose (0 or 1) -j jobID (0 to ...) '
- sys.exit()
- for o, a in opts:
- if o in ("-f",):
- fileName = a
- if o in ("-o",):
- outputFileName = a
- if o in ("-n",):
- print "-n",a
- maxevents = int(a)
- if o in ("-j",):
- print "-j",a
- jobID = int(a)
- if o in ("-t",):
- sampleType = str(a)
- if o in ("-v",):
- verbose = bool(int(a))
- if o in ("--Ethr=",):
- threshold = float(a)
-
- fileNameGeo = fileName.replace('ship', 'geofile_full').replace("_rec","")
-
- def namestr(obj, namespace=globals()):
- return [name for name in namespace if namespace[name] is obj]
-
- for item in [fileName, fileNameGeo, outputFileName, sampleType]:
- if not item:
- print '\tFATAL! %s not defined!'%namestr(item)[0]
- sys.exit()
-
- for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose, threshold]:
- print '\tINFO: %s set to %s'%(namestr(item)[0], item)
-
- # Load SHiP (LHCb) style
- gROOT.ProcessLine(".x mystyle.C")
- # Obsolete debug flags
- oldGeo = False
- ON = True
- # Load PDG database
- pdg = TDatabasePDG.Instance()
- # Add elements to PDG database
- import pythia8_conf
- pythia8_conf.addHNLtoROOT()
- # TODO: MAKE THIS A DICTIONARY AND DE-VERBOSIZE IT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
- 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",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):
- tmpName = pdg.GetParticle(RPCpt.PdgCode())
- print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName()
- print
-
- # weight computation moved to offlineForBarbara.py
-
- 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
- foundStat = True
- eff_val = gRandom.Uniform(0.,1.)
- if eff_val<0.9:
- flag_eff = flag_eff or True
- foundStat_eff = True
- else:
- flag_eff = flag_eff or False
- if foundStat:
- nstat+=1
- if foundStat_eff:
- nstat_eff+=1
- particles = 0
- 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"
- 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):
- tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName()
- tmpName = convertion(tmpName)
- if hit.GetZ()>8000:
- continue
- 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()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]):
- continue
- # 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
-
-
- """ Main program """
-
- # Open file
- f = TFile(fileName)
- t = f.Get("cbmsim")
-
- totentries = t.GetEntries()
- if verbose:
- print
- print
- print "Program started, reading %s from %s"%(t.GetName(), fileName)
- if (maxevents>0) and (maxevents < totentries):
- entries = maxevents
- else:
- entries = totentries
- if verbose:
- print "Requested to process %s events out of %s in tree"%(entries, totentries)
-
- # Load geometry
- from lookAtGeo import *
- r = loadGeometry(fileNameGeo)
- fGeo = r['fGeo']
-
- # Functions by Elena for offline selection
- # (Also includes another dictionary of all geometry nodes)
- import offlineForBarbara as offline
- offline.loadNodes(fGeo)
- offline.ShipGeo = r['ShipGeo']
- offline.initBField(fileNameGeo)
- offline.sh = offline.StrawHits(t, offline.modules, offline.ShipGeo.straw.resol, 0, None, offline.ShipGeo)
- offline.pdg = pdg
-
- # Handle geometry
- # Names of useful geometry nodes
- myNodes_name = ["VetoTimeDet_1"]
- myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
- myNodes_name += ["Veto_5"]
- # Positions of selected nodes
- myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None)
- # Places where a neutrino can interact
- 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"]
- # Z positions of the VETO and tracking systems
- 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]]
- # Places where a neutrino can interact (some duplicates?)
- 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"]
- # Printout positions
- print "OPERApos: ",RPCstationsPos
- print "trackingPos: ",trackStationsPos
- print "strawVetoPos: ",strawVetoPos
- print "scintPos: ",vetoWall
- print "upstreamVetoPos: ",upstreamVetoPos
-
- # Prepare output ntuple
- nf = TFile(outputFileName,"RECREATE")
- nt = TTree("t","t")
- # Event number --------------------------------------------
- nt, event = addVar(nt, 'event','i')
- # Neutrino information ------------------------------------
- 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')
- # Interaction element code --------------------------------
- ## 0: last passive material opera-mu-system
- ## 1: two windows around liquid scintillator
- ## 2: 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')
- # Neutrino daughters --------------------------------------
- # Do the particles come from a 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')
- # Do charged particles come from a neutrino?
- nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i')
- nt, idChargedPart_fromNu = addVect(nt, 'idChargedPart_fromNu', 'int')
- #nt, idStrChargedPart_fromNu = addVect(nt, 'idStrChargedPart_fromNu', 'string')
- # Do neutral particles come from a neutrino?
- nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i')
- nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int')
- #nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string')
- # RPC -----------------------------------------------------
- # In case something (no matter if it comes from the nu-interaction kids) fired the RPC
- nt, RPCany = addVar(nt,'RPCany', 'i' )
- # accounting for an eff. of each station of 90%
- nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' )
- # upstreamVeto --------------------------------------------
- nt, upstreamVetoany = addVar(nt,'upstreamVetoAny', 'i' )
- # accounting for an eff. of each station of 90%
- nt, upstreamVetoany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' )
- # scintVeto -----------------------------------------------
- ## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls
- nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' )
- nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' )
- # strawVeto -----------------------------------------------
- # In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto
- 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')
- # VETO and tracking systems with thresholds ---------------
- nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i")
- nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i")
- nt, upstreamVetoAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i")
- nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i")
- nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i")
- # Decay information ---------------------------------------
- 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")
- # Is this a NC interaction? -------------------------------
- nt, NC = addVar(nt, "NC", "i")
- # Store where the neutrino interacted ---------------------
- #nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string")
- nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","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}
- # Was the event reconstructed? ----------------------------
- nt, recoed = addVar(nt, "recoed","i")
- # How many candidates in the reco event? ------------------
- nt, nRecoed = addVar(nt, "nRecoed","i")
- # Add stuff for online selection --------------------------
- offline.addOfflineToTree(nt)
-
- # Now loop on the tree
- # t = original tree
- # nt = new tree
- for entry in xrange(entries):
- if not (entry%1000):
- print "\tProcessing entry %s of %s..."%(entry,entries)
- t.GetEntry(entry)
- event[0] = entry+(jobID*entries)
- particles = t.MCTrack
- rpc = t.ShipRpcPoint
- scint = t.vetoPoint
- strawPoints = t.strawtubesPoint
- vetoScint = t.vetoPoint
- recoParts = t.Particles
- # initialize containers
- putToZero(nNu)
- putToZero(nPart_fromNu)
- putToZero(nChargedPart_fromNu)
- putToZero(nNeutrPart_fromNu)
- idPart_fromNu.clear()
- idChargedPart_fromNu.clear()
- idNeutrPart_fromNu.clear()
- #idStrPart_fromNu.clear()
- #idStrChargedPart_fromNu.clear()
- #idStrNeutrPart_fromNu.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
- # Which one is the "primary" particle?
- if sampleType=="nuBg" or sampleType == "cosmics":
- startPartRange = [0]
- primaryMum = -1
- elif sampleType=="sig":
- startPartRange = [1]
- primaryMum = 0
- # Look at the primary particle
- ip = startPartRange[0]
- skipEvent=False
- # Were there any hit in the tracking stations?
- TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, None, None, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] )
- # If no, skip this event
- if TrackSyst[0]==0:
- skipEvent=True
- #print "i am jumping this one"
- continue
- # exit from loop on particles
- #break
- # Select the primary MC particle
- part = particles[ip]
- pdgPart = pdg.GetParticle(part.GetPdgCode())
- if sampleType=="nuBg":
- assert("nu" in pdgPart.GetName())
- ## 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
- ## (and de-indented the following block)
- ## Starting the counter of how many particles were produced by the interaction of this specific nu
- for recoP in recoParts:
- nRecoed[0] += 1
- #offline.pushOfflineByParticle(t, recoP)
- if nRecoed[0]>0:
- recoed[0]=1
- else:
- skipEvent=True
- continue
- #break
- nNu[0]+=1
- if part.GetMotherId()==primaryMum:
- assert(primaryDone==False)
- primaryDone=True
- isPrimary[0] = int(True)
- assert(len(idPart_fromNu)==0)
- assert(len(idNeutrPart_fromNu)==0)
- else:
- continue
- # Where did this guy interact?
- tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName()
- #nuInteractionNode.push_back(tmpName)
- tmpName = convertion(tmpName)
- if not tmpName in dictNodeNames:
- tmpName = "others"
- nuIntNumSimpl[0] = dictNodeNames[tmpName]
- # Store interaction point coordinates
- startZ_nu[0] = part.GetStartZ()
- startY_nu[0] = part.GetStartY()
- startX_nu[0] = part.GetStartX()
- # Primary particle information
- nuE[0] = part.GetEnergy()
- # Looking for particles produced by the neutrino interaction
- interacted = False
- partKidsId = []
- NC[0] = 0
-
- # Add information for offline selection
- # (mainly signal normalisation and zeroing of particle-wise arrays)
- ntr, nref = offline.pushOfflineByEvent(t, vetoScint, sampleType, verbose, threshold)
- for recoP in recoParts:
- offline.pushOfflineByParticle(t, recoP, ntr, nref)
-
- # Loop on MCTracks
- 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)
- # Was this a neutral current interaction?
- if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0):
- NC[0] = 1
- 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)
- else:
- nChargedPart_fromNu[0]+=1
- idChargedPart_fromNu.push_back(part2Id)
- #idStrChargedPart_fromNu.push_back(part2Name)
- # How many are produced by the interaction in the last passive element of the "opera-mu system"?
- # Finding out the "interaction element code"
- part2Z = part2.GetStartZ()
- part2X = part2.GetStartX()
- part2Y = part2.GetStartY()
- somewhere = False
- nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName()
- if nodeName in lastPassive_nodeName:
- somewhere = True
- interactionElement[0] = 0
- intElName = 'OPERA'
- # To know if it was in the full OPERA-system (excluded last passive)
- if tmpName in OPERA_nodeName and somewhere==False:
- somewhere = True
- interactionElement[0] = 3
- intElName = 'OPERA'
- # To know if it was between the two windows
- if tmpName in entrance_nodeName:
- assert(somewhere==False)
- somewhere = True
- interactionElement[0] = 1
- # Vacuum tank outer window
- if tmpName in volumeIn_nodeName and somewhere==False:
- interactionElement[0] = 5
- somewhere = True
- # Vacuum tank inner window
- elif nodeName in volumeOut_nodeName and somewhere==False:
- interactionElement[0] = 6
- somewhere = True
- # Vacuum tank ribs
- if "Rib" in tmpName:
- interactionElement[0] = 7
- # Somewhere else
- if somewhere==False:
- interactionElement[0] = -1
- # End loop on MCTracks
- # Should be changed to avoid entering in the loop if there isn't anything in the tracking
- if skipEvent:
- continue
- 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])
- upstreamVetoany[0], upstreamVetoany_eff[0], dummy,dummy, upstreamVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, upstreamVetoPos, None, None, checkOn=False, pointVects=None)#[upstreamVetoPoint_x, upstreamVetoPoint_y, upstreamVetoPoint_z])
- 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])
- scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=threshold)
- #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 no primary nu interacting"
- PrintEventPart(particles,pdg)
- sys.exit()
- weight[0] = offline.findWeight(sampleType, NC[0], nuE[0], part, entries, pdgPart.GetName(), ON)#calcWeight(NC[0], nuE[0], part.GetWeight(), entries, pdgPart.GetName(), ON)
- nt.Fill()
- # End loop on tree
- nt.Write()
- nf.Save()
- nf.Close()
- f.Close()
- print "Output wrote to %s"%outputFileName
- print "Number of events with vertex outside Veto_5-500 - Tr4_4: %s"%offline.num_bad_z