Newer
Older
FairShipTools / ntuple_nuVeto_small_newGen_backup
@Ubuntu Ubuntu on 22 Mar 2015 28 KB software run on yandex
  1. import sys, re, getopt
  2. from array import array
  3. from ROOT import *
  4.  
  5. fileName, fileNameGeo, outputFileName, sampleType, jobID = None, None, None, None, None
  6. maxevents = 0
  7. verbose = True
  8. # Energy deposit threshold for the liquid scintillator
  9. threshold = 0.045
  10.  
  11. # Parse commandline inputs
  12. try:
  13. opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:j:", ['Ethr='])
  14. except getopt.GetoptError:
  15. # print help information and exit:
  16. 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 ...) '
  17. sys.exit()
  18. for o, a in opts:
  19. if o in ("-f",):
  20. fileName = a
  21. if o in ("-o",):
  22. outputFileName = a
  23. if o in ("-n",):
  24. print "-n",a
  25. maxevents = int(a)
  26. if o in ("-j",):
  27. print "-j",a
  28. jobID = int(a)
  29. if o in ("-t",):
  30. sampleType = str(a)
  31. if o in ("-v",):
  32. verbose = bool(int(a))
  33. if o in ("--Ethr=",):
  34. threshold = float(a)
  35. fileNameGeo = fileName.replace('ship', 'geofile_full').replace("_rec","")
  36.  
  37. def namestr(obj, namespace=globals()):
  38. return [name for name in namespace if namespace[name] is obj]
  39.  
  40. for item in [fileName, fileNameGeo, outputFileName, sampleType]:
  41. if not item:
  42. print '\tFATAL! %s not defined!'%namestr(item)[0]
  43. sys.exit()
  44.  
  45. for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose, threshold]:
  46. print '\tINFO: %s set to %s'%(namestr(item)[0], item)
  47.  
  48. # Load SHiP (LHCb) style
  49. gROOT.ProcessLine(".x mystyle.C")
  50. # Obsolete debug flags
  51. oldGeo = False
  52. ON = True
  53. # Load PDG database
  54. pdg = TDatabasePDG.Instance()
  55. # Add elements to PDG database
  56. import pythia8_conf
  57. pythia8_conf.addHNLtoROOT()
  58. # TODO: MAKE THIS A DICTIONARY AND DE-VERBOSIZE IT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
  59. if not(pdg.GetParticle(1000010020)):
  60. pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020)
  61. pdg.AddAntiParticle("AntiDeuteron", - 1000010020)
  62. if not(pdg.GetParticle(1000010030)):
  63. pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030);
  64. pdg.AddAntiParticle("AntiTriton", - 1000010030);
  65. if not(pdg.GetParticle(1000020040) ):
  66. pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040);
  67. pdg.AddAntiParticle("AntiAlpha", - 1000020040);
  68. if not(pdg.GetParticle(1000020030) ):
  69. pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030);
  70. pdg.AddAntiParticle("AntiHE3", -1000020030);
  71. if not(pdg.GetParticle(1000030070) ):
  72. print "TO BE CHECKED the data insert for Li7Nucleus"
  73. pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070);
  74. if not(pdg.GetParticle(1000060120) ):
  75. print "ERROR: random values insert for C12Nucleus"
  76. pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120);
  77. if not(pdg.GetParticle(1000030060) ):
  78. print "ERROR: random values insert for Li6Nucleus"
  79. pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060);
  80. if not(pdg.GetParticle(1000070140) ):
  81. print "ERROR: random values insert for for N14"
  82. pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140);
  83. if not(pdg.GetParticle(1000050100) ):
  84. print "TO BE CHECKED the data insert for B10"
  85. pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100);
  86. if not(pdg.GetParticle(1000020060) ):
  87. print "TO BE CHECKED the data insert for He6"
  88. pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060);
  89. if not(pdg.GetParticle(1000040080) ):
  90. print "TO BE CHECKED the data insert for Be8"
  91. pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080);
  92. if not(pdg.GetParticle(1000030080) ):
  93. print "TO BE CHECKED the data insert for Li8"
  94. pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080);
  95. if not(pdg.GetParticle(1000040170) ):
  96. print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!"
  97. pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170);
  98. if not(pdg.GetParticle(1000040100) ):
  99. print "TO BE CHECKED the data insert for Be10"
  100. pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100);
  101. if not(pdg.GetParticle(1000040070) ):
  102. print "TO BE CHECKED the data insert for Be7"
  103. pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070);
  104. if not(pdg.GetParticle(1000230470) ):
  105. print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!"
  106. pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470);
  107. if not(pdg.GetParticle(1000080170) ):
  108. print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!"
  109. pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170);
  110. if not(pdg.GetParticle(1000240500) ):
  111. print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!"
  112. pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500);
  113. if not(pdg.GetParticle(1000210450) ):
  114. print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!"
  115. pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450);
  116. if not(pdg.GetParticle(1000040090) ):
  117. print "TO BE CHECKED the data insert for Be9"
  118. pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090);
  119. if not(pdg.GetParticle(1000080160) ):
  120. print "TO BE CHECKED the data insert for O16"
  121. pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160);
  122. if not(pdg.GetParticle(1000220460) ):
  123. print "TO BE CHECKED the data insert for Ar40"
  124. pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460);
  125. if not(pdg.GetParticle(1000050110) ):
  126. print "TO BE CHECKED the data insert for B11"
  127. pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110);
  128.  
  129. def PrintEventPart(particles,pdg):
  130. print "Particles in the event:"
  131. for (pid,part) in enumerate(particles):
  132. print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId()
  133. print
  134.  
  135. def PrintRPChits(rpc,pdg):
  136. print "Hits in:"
  137. for (ix,RPCpt) in enumerate(rpc):
  138. tmpName = pdg.GetParticle(RPCpt.PdgCode())
  139. print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName()
  140. print
  141.  
  142. # weight computation moved to offlineForBarbara.py
  143.  
  144. def addVect(t,name,vectType):
  145. vect = vector(vectType)()
  146. t.Branch( name, vect )
  147. return t, vect
  148.  
  149. def addVar(t,name,varType):
  150. var = array(varType,[-999])
  151. t.Branch(name,var,name+"/"+varType.upper())
  152. return t, var
  153.  
  154. def putToZero(var):
  155. var[0] = 0
  156. def getPartName(partId):
  157. return pdg.GetParticle(partId).GetName()
  158.  
  159. def lookingForDecay(particles):
  160. decayMotherList = []
  161. for p in particles:
  162. mID = p.GetMotherId()
  163. if mID>=0 and not (mID in decayMotherList):
  164. decayMotherList.append(mID)
  165. decayPartIndex.push_back(mID)
  166. decayPartID.push_back(particles[mID].GetPdgCode())
  167. decayPartStrID.push_back(pdg.GetParticle(particles[mID].GetPdgCode()).GetName())
  168. decayPos_x.push_back(p.GetStartX())
  169. decayPos_y.push_back(p.GetStartY())
  170. decayPos_z.push_back(p.GetStartZ())
  171. decayMother.push_back(particles[mID].GetMotherId())
  172. decayP.push_back(p.GetP())
  173.  
  174. def wasFired(indexKids, detPoints, detPos, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):
  175. #
  176. def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId):
  177. charges = []
  178. trackids = []
  179. if hitCharges is None:
  180. assert(hitTrackId is None)
  181. for pos in detPos:
  182. foundStat = False
  183. foundStat_eff = False
  184. for hit in detPoints:
  185. if (indexKids is None) or (hit.GetTrackID() in indexKids):
  186. # check if it is in one of the considered active layer
  187. if pos[0]<=hit.GetZ()<=pos[1]:
  188. Eloss += hit.GetEnergyLoss()
  189. hitID = hit.GetTrackID()
  190. if not hitID in trackids and not hitCharges is None:
  191. if hit.PdgCode()>100000:
  192. charges.append(9)
  193. else:
  194. charges.append(int(pdg.GetParticle(hit.PdgCode()).Charge()))
  195. trackids.append(hitID)
  196. hitCharges.push_back(charges[-1])
  197. hitTrackId.push_back(trackids[-1])
  198. if pointVects is not None:
  199. pointVects[0].push_back(hit.GetX())
  200. pointVects[1].push_back(hit.GetY())
  201. pointVects[2].push_back(hit.GetZ())
  202. flag = True
  203. foundStat = True
  204. eff_val = gRandom.Uniform(0.,1.)
  205. if eff_val<0.9:
  206. flag_eff = flag_eff or True
  207. foundStat_eff = True
  208. else:
  209. flag_eff = flag_eff or False
  210. if foundStat:
  211. nstat+=1
  212. if foundStat_eff:
  213. nstat_eff+=1
  214. particles = 0
  215. flag_Ethr = Eloss>=Ethr
  216. return flag, flag_eff, nstat, nstat_eff, Eloss, flag_Ethr,particles
  217. #
  218. # Now in partKidTrackID I should have the trackID connected to my charged particles
  219. flag_eff = False
  220. flag = False
  221. nstat = 0
  222. nstat_eff = 0
  223. Eloss = 0
  224. flag,flag_eff,nstat,nstat_eff,Eloss,flag_Ethr,particles = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId)
  225. if flag==False and checkOn:
  226. print "To be study event %s"%entry
  227. return flag, flag_eff, nstat, nstat_eff, flag_Ethr, particles
  228.  
  229. def convertion(tmpName):
  230. if "Rib" in tmpName: tmpName = "Rib"
  231. elif "LiSc" in tmpName: tmpName= "LiSc"
  232. elif "Tr" in tmpName: tmpName = "Tr"
  233. elif "Startplate" in tmpName: tmpName = "Startplate"
  234. elif re.search("T\d+O", tmpName): tmpName = "TO"
  235. elif re.search("T\d+I", tmpName): tmpName = "TI"
  236. elif "Endplate" in tmpName: tmpName = "Endplate"
  237. elif "volCoil" in tmpName: tmpName = "volCoil"
  238. elif "volFeYoke" in tmpName: tmpName = "volFeYoke"
  239. elif "gas" in tmpName: tmpName = "gas"
  240. elif "wire" in tmpName: tmpName == "wire"
  241. elif "VetoTimeDet" in tmpName: tmpName = "upstreamVeto"
  242. elif "Veto" in tmpName: tmpName = "strawVeto"
  243. return tmpName
  244.  
  245. def wasFired_node(indexKids, detPoints, detVolNames, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):
  246. #
  247. def lookingForHits(detPoints,flag,hitCharges, hitTrackId):
  248. if hitCharges is None:
  249. assert(hitTrackId is None)
  250. for hit in detPoints:
  251. if (indexKids is None) or (hit.GetTrackID() in indexKids):
  252. tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName()
  253. tmpName = convertion(tmpName)
  254. if hit.GetZ()>8000:
  255. continue
  256. if tmpName in detVolNames:
  257. ## trick to remove the case of the tracking system the hits in the gas or straw from the straw-veto:
  258. if 'trackingSystem' in detVolNames and hit.GetZ()<0:
  259. continue
  260. if "strawVetoSystem" in detVolNames and hit.GetZ()>0:
  261. continue
  262. if "upstreamVeto" in detVolNames and not(hit.GetZ()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]):
  263. continue
  264. # check if it is in one of the considered active layer
  265. if True: #pos[0]<=hit.GetZ()<=pos[1]:
  266. flag = True
  267. return flag
  268. #
  269. # Now in partKidTrackID I should have the trackID connected to my charged particles
  270. flag = False
  271. flag = lookingForHits(detPoints,flag,hitCharges, hitTrackId)
  272. if flag==False and checkOn:
  273. print "To be study event %s"%entry
  274. return flag
  275.  
  276.  
  277. """ Main program """
  278.  
  279. # Open file
  280. f = TFile(fileName)
  281. t = f.Get("cbmsim")
  282.  
  283. totentries = t.GetEntries()
  284. if verbose:
  285. print
  286. print
  287. print "Program started, reading %s from %s"%(t.GetName(), fileName)
  288. if (maxevents>0) and (maxevents < totentries):
  289. entries = maxevents
  290. else:
  291. entries = totentries
  292. if verbose:
  293. print "Requested to process %s events out of %s in tree"%(entries, totentries)
  294.  
  295. # Load geometry
  296. from lookAtGeo import *
  297. r = loadGeometry(fileNameGeo)
  298. fGeo = r['fGeo']
  299.  
  300. # Functions by Elena for offline selection
  301. # (Also includes another dictionary of all geometry nodes)
  302. import offlineForBarbara as offline
  303. offline.loadNodes(fGeo)
  304. offline.ShipGeo = r['ShipGeo']
  305. offline.initBField(fileNameGeo)
  306. offline.sh = offline.StrawHits(t, offline.modules, offline.ShipGeo.straw.resol, 0, None, offline.ShipGeo)
  307. offline.pdg = pdg
  308.  
  309. # Handle geometry
  310. # Names of useful geometry nodes
  311. myNodes_name = ["VetoTimeDet_1"]
  312. myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
  313. myNodes_name += ["Veto_5"]
  314. # Positions of selected nodes
  315. myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None)
  316. # Places where a neutrino can interact
  317. lastPassive_nodeName = ["volIron_23"]
  318. 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"]
  319. entrance_nodeName = ["T1Lid"]#['T1lid_1']
  320. volumeIn_nodeName = ["TI"]#['T%sI'%i for i in [1,2,3,5]]
  321. volumeOut_nodeName = ["TO"]#['T%sO'%i for i in [1,2,3,5]]
  322. OPERA_volNames = ["volIron","volRpc","volHPT"]
  323. strawVeto_volNames = ["Veto"]
  324. # Z positions of the VETO and tracking systems
  325. Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']]
  326. upstreamVeto = [myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]
  327. 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)]
  328. strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]]
  329. 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.]]
  330. upstreamVetoPos = [[myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]]
  331. RPCstationsPos = [[-3500,myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ']-0.5]]
  332. # Places where a neutrino can interact (some duplicates?)
  333. trackStationsPos_node = ["Tr", "gas", "straw", 'trackingSystem']
  334. ## should be that this does not work since it could be that hits are also assigned to gas or straw
  335. strawVetoPos_node = ["strawVeto", "gas", "straw", "strawVetoSystem"]
  336. vetoWall_node = ["LiSc","cave","Rib","TI","TO","Tr","strawVeto"]
  337. upstreamVetoPos_node = ["upstreamVeto","cave"]
  338. RPCstationsPos_node = ["volRpc","volMagneticSpectrometer","volIron","volHPT"]
  339. # Printout positions
  340. print "OPERApos: ",RPCstationsPos
  341. print "trackingPos: ",trackStationsPos
  342. print "strawVetoPos: ",strawVetoPos
  343. print "scintPos: ",vetoWall
  344. print "upstreamVetoPos: ",upstreamVetoPos
  345.  
  346. # Prepare output ntuple
  347. nf = TFile(outputFileName,"RECREATE")
  348. nt = TTree("t","t")
  349. # Event number --------------------------------------------
  350. nt, event = addVar(nt, 'event','i')
  351. # Neutrino information ------------------------------------
  352. nt, nNu = addVar(nt, 'nNu', 'i')
  353. nt, isPrimary = addVar(nt, 'isPrimary','i')
  354. nt, startZ_nu = addVar(nt, 'startZ_nu', 'f')
  355. nt, startY_nu = addVar(nt, 'startY_nu', 'f')
  356. nt, startX_nu = addVar(nt, 'startX_nu', 'f')
  357. nt, nuE = addVar(nt, 'nuE', 'f')
  358. nt, weight = addVar(nt, 'weight', 'f')
  359. # Interaction element code --------------------------------
  360. ## 0: last passive material opera-mu-system
  361. ## 1: two windows around liquid scintillator
  362. ## 2: along vaccum tank
  363. ## 3: along all OPERA
  364. ## 4: between the two entrance windows
  365. ## 5: along vacuum tank outer window
  366. ## 6: along vacuum tank inner window
  367. ## 999: wrong OPERA place ## needed until we have the new production
  368. ## -1: anything else, but what???? #it should never been present
  369. nt, interactionElement = addVar(nt, 'interactionElement','i')
  370. # Neutrino daughters --------------------------------------
  371. # Do the particles come from a neutrino?
  372. nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i')
  373. nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int')
  374. #nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string')
  375. # Do charged particles come from a neutrino?
  376. nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i')
  377. nt, idChargedPart_fromNu = addVect(nt, 'idChargedPart_fromNu', 'int')
  378. #nt, idStrChargedPart_fromNu = addVect(nt, 'idStrChargedPart_fromNu', 'string')
  379. # Do neutral particles come from a neutrino?
  380. nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i')
  381. nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int')
  382. #nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string')
  383. # RPC -----------------------------------------------------
  384. # In case something (no matter if it comes from the nu-interaction kids) fired the RPC
  385. nt, RPCany = addVar(nt,'RPCany', 'i' )
  386. # accounting for an eff. of each station of 90%
  387. nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' )
  388. # upstreamVeto --------------------------------------------
  389. nt, upstreamVetoany = addVar(nt,'upstreamVetoAny', 'i' )
  390. # accounting for an eff. of each station of 90%
  391. nt, upstreamVetoany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' )
  392. # scintVeto -----------------------------------------------
  393. ## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls
  394. nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' )
  395. nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' )
  396. # strawVeto -----------------------------------------------
  397. # In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto
  398. nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' )
  399. nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' )
  400. # TrackSyst -----------------------------------------------
  401. # only the case of anything fired it is considered, obviously
  402. nt, TrackSyst = addVar(nt, "TrackSyst", 'i')
  403. nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i')
  404. # VETO and tracking systems with thresholds ---------------
  405. nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i")
  406. nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i")
  407. nt, upstreamVetoAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i")
  408. nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i")
  409. nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i")
  410. # Decay information ---------------------------------------
  411. nt, decayPartIndex = addVect(nt, "decayPartIndex", "float")
  412. nt, decayPartID = addVect(nt, "decayPartID", "float")
  413. nt, decayPartStrID = addVect(nt, "decayPartStrID", "string")
  414. nt, decayPos_x = addVect(nt, "decayPos_x", "float")
  415. nt, decayPos_y = addVect(nt, "decayPos_y", "float")
  416. nt, decayPos_z = addVect(nt, "decayPos_z", "float")
  417. nt, decayMother = addVect(nt, "decayMother", "float")
  418. nt, decayP = addVect(nt,"decayP","float")
  419. # Is this a NC interaction? -------------------------------
  420. nt, NC = addVar(nt, "NC", "i")
  421. # Store where the neutrino interacted ---------------------
  422. #nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string")
  423. nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","i")
  424. 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}
  425. # Was the event reconstructed? ----------------------------
  426. nt, recoed = addVar(nt, "recoed","i")
  427. # How many candidates in the reco event? ------------------
  428. nt, nRecoed = addVar(nt, "nRecoed","i")
  429. # Add stuff for online selection --------------------------
  430. offline.addOfflineToTree(nt)
  431.  
  432. # Now loop on the tree
  433. # t = original tree
  434. # nt = new tree
  435. for entry in xrange(entries):
  436. if not (entry%1000):
  437. print "\tProcessing entry %s of %s..."%(entry,entries)
  438. t.GetEntry(entry)
  439. event[0] = entry+(jobID*entries)
  440. particles = t.MCTrack
  441. rpc = t.ShipRpcPoint
  442. scint = t.vetoPoint
  443. strawPoints = t.strawtubesPoint
  444. vetoScint = t.vetoPoint
  445. recoParts = t.Particles
  446. # initialize containers
  447. putToZero(nNu)
  448. putToZero(nPart_fromNu)
  449. putToZero(nChargedPart_fromNu)
  450. putToZero(nNeutrPart_fromNu)
  451. idPart_fromNu.clear()
  452. idChargedPart_fromNu.clear()
  453. idNeutrPart_fromNu.clear()
  454. #idStrPart_fromNu.clear()
  455. #idStrChargedPart_fromNu.clear()
  456. #idStrNeutrPart_fromNu.clear()
  457. decayPartIndex.clear()
  458. decayPartID.clear()
  459. decayPartStrID.clear()
  460. decayPos_x.clear()
  461. decayPos_y.clear()
  462. decayPos_z.clear()
  463. decayMother.clear()
  464. decayP.clear()
  465. #nuInteractionNode.clear()
  466. primaryDone=False
  467. isPrimary[0] = int(False)
  468. lookingForDecay(particles)
  469. nRecoed[0] = 0
  470. recoed[0] = 0
  471. # Which one is the "primary" particle?
  472. if sampleType=="nuBg" or sampleType == "cosmics":
  473. startPartRange = [0]
  474. primaryMum = -1
  475. elif sampleType=="sig":
  476. startPartRange = [1]
  477. primaryMum = 0
  478. # Look at the primary particle
  479. ip = startPartRange[0]
  480. skipEvent=False
  481. # Were there any hit in the tracking stations?
  482. 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] )
  483. # If no, skip this event
  484. if TrackSyst[0]==0:
  485. skipEvent=True
  486. #print "i am jumping this one"
  487. continue
  488. # exit from loop on particles
  489. #break
  490. # Select the primary MC particle
  491. part = particles[ip]
  492. pdgPart = pdg.GetParticle(part.GetPdgCode())
  493. if sampleType=="nuBg":
  494. assert("nu" in pdgPart.GetName())
  495. ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother
  496. #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena
  497. ## (and de-indented the following block)
  498. ## Starting the counter of how many particles were produced by the interaction of this specific nu
  499. for recoP in recoParts:
  500. nRecoed[0] += 1
  501. #offline.pushOfflineByParticle(t, recoP)
  502. if nRecoed[0]>0:
  503. recoed[0]=1
  504. else:
  505. skipEvent=True
  506. continue
  507. #break
  508. nNu[0]+=1
  509. if part.GetMotherId()==primaryMum:
  510. assert(primaryDone==False)
  511. primaryDone=True
  512. isPrimary[0] = int(True)
  513. assert(len(idPart_fromNu)==0)
  514. assert(len(idNeutrPart_fromNu)==0)
  515. else:
  516. continue
  517. # Where did this guy interact?
  518. tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName()
  519. #nuInteractionNode.push_back(tmpName)
  520. tmpName = convertion(tmpName)
  521. if not tmpName in dictNodeNames:
  522. tmpName = "others"
  523. nuIntNumSimpl[0] = dictNodeNames[tmpName]
  524. # Store interaction point coordinates
  525. startZ_nu[0] = part.GetStartZ()
  526. startY_nu[0] = part.GetStartY()
  527. startX_nu[0] = part.GetStartX()
  528. # Primary particle information
  529. nuE[0] = part.GetEnergy()
  530. # Looking for particles produced by the neutrino interaction
  531. interacted = False
  532. partKidsId = []
  533. NC[0] = 0
  534.  
  535. # Add information for offline selection
  536. # (mainly signal normalisation and zeroing of particle-wise arrays)
  537. ntr, nref = offline.pushOfflineByEvent(t, vetoScint, sampleType, verbose, threshold)
  538. for recoP in recoParts:
  539. offline.pushOfflineByParticle(t, recoP, ntr, nref)
  540.  
  541. # Loop on MCTracks
  542. for ip2 in xrange(0,len(particles)):
  543. part2 = particles[ip2]
  544. # exit if we have reached the empty part of the array
  545. if not (type(part2)==type(ShipMCTrack())):
  546. break
  547. if part2.GetMotherId()==ip:
  548. interacted = True
  549. part2Id = part2.GetPdgCode()
  550. partKidsId.append(part2Id)
  551. # Was this a neutral current interaction?
  552. if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0):
  553. NC[0] = 1
  554. nPart_fromNu[0]+=1
  555. idPart_fromNu.push_back(part2Id)
  556. if pdg.GetParticle(part2.GetPdgCode()) == None:
  557. part2Name = str(part2Id)
  558. else:
  559. part2Name = getPartName(part2Id)
  560. #idStrPart_fromNu.push_back(part2Name)
  561. # Counting how many particles produced by the nu-interaction are charged or neutral
  562. if (pdg.GetParticle(part2.GetPdgCode())) == None or (int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0)):
  563. nNeutrPart_fromNu[0]+=1
  564. idNeutrPart_fromNu.push_back(part2Id)
  565. #idStrNeutrPart_fromNu.push_back(part2Name)
  566. else:
  567. nChargedPart_fromNu[0]+=1
  568. idChargedPart_fromNu.push_back(part2Id)
  569. #idStrChargedPart_fromNu.push_back(part2Name)
  570. # How many are produced by the interaction in the last passive element of the "opera-mu system"?
  571. # Finding out the "interaction element code"
  572. part2Z = part2.GetStartZ()
  573. part2X = part2.GetStartX()
  574. part2Y = part2.GetStartY()
  575. somewhere = False
  576. nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName()
  577. if nodeName in lastPassive_nodeName:
  578. somewhere = True
  579. interactionElement[0] = 0
  580. intElName = 'OPERA'
  581. # To know if it was in the full OPERA-system (excluded last passive)
  582. if tmpName in OPERA_nodeName and somewhere==False:
  583. somewhere = True
  584. interactionElement[0] = 3
  585. intElName = 'OPERA'
  586. # To know if it was between the two windows
  587. if tmpName in entrance_nodeName:
  588. assert(somewhere==False)
  589. somewhere = True
  590. interactionElement[0] = 1
  591. # Vacuum tank outer window
  592. if tmpName in volumeIn_nodeName and somewhere==False:
  593. interactionElement[0] = 5
  594. somewhere = True
  595. # Vacuum tank inner window
  596. elif nodeName in volumeOut_nodeName and somewhere==False:
  597. interactionElement[0] = 6
  598. somewhere = True
  599. # Vacuum tank ribs
  600. if "Rib" in tmpName:
  601. interactionElement[0] = 7
  602. # Somewhere else
  603. if somewhere==False:
  604. interactionElement[0] = -1
  605. # End loop on MCTracks
  606. # Should be changed to avoid entering in the loop if there isn't anything in the tracking
  607. if skipEvent:
  608. continue
  609. 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])
  610. 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])
  611. 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])
  612. scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=threshold)
  613. #assert(len(idStrPart_fromNu)==len(idPart_fromNu))
  614. #assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu))
  615. #assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu))
  616. #assert(len(idStrPart_fromNu)==nChargedPart_fromNu[0]+nNeutrPart_fromNu[0])
  617. #assert(len(idStrChargedPart_fromNu)==nChargedPart_fromNu[0])
  618. #assert(len(idStrNeutrPart_fromNu)==nNeutrPart_fromNu[0])
  619. if not primaryDone:
  620. print "It looks like there is no primary nu interacting"
  621. PrintEventPart(particles,pdg)
  622. sys.exit()
  623. 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)
  624. nt.Fill()
  625. # End loop on tree
  626. nt.Write()
  627. nf.Save()
  628. nf.Close()
  629. f.Close()
  630. print "Output wrote to %s"%outputFileName
  631. print "Number of events with vertex outside Veto_5-500 - Tr4_4: %s"%offline.num_bad_z