Newer
Older
FairShipTools / funcsByBarbara.py
  1. from elena import *
  2. from array import array
  3. from collections import OrderedDict
  4. import ROOT
  5.  
  6. # Custom sorting
  7. def sortNodesBy(nodes, coord='z', what='pos', printdict=True):
  8. d = OrderedDict(sorted(nodes.items(), key = lambda item: item[1][coord][what]))
  9. if printdict:
  10. print
  11. print
  12. print "Nodes by %s %s:"%(coord, what)
  13. for node in d.keys():
  14. print node, '\t\t', (d[node][coord]['pos']-d[node][coord]['dim'])/100., ' m', '\t', 2*d[node][coord]['dim']/100., ' m'
  15. return d
  16.  
  17. def retrieveMCParticleInfo(sTree, key):
  18. mcPartKey = sTree.fitTrack2MC[key]
  19. mcPart = sTree.MCTrack[mcPartKey]
  20. if not mcPart:
  21. print "Error in retrieveMCParticleInfo"
  22. sys.exit()
  23. mother = sTree.MCTrack[mcPart.GetMotherId()]
  24. return int(mcPart.GetPdgCode()), int(mother.GetPdgCode()), mcPart.GetStartX(), mcPart.GetStartY(), mcPart.GetStartZ()
  25.  
  26.  
  27. def AddVect(t,name,vectType):
  28. vect = ROOT.vector(vectType)()
  29. t.Branch( name, vect )
  30. return t, vect
  31.  
  32. def AddVar(t,name,varType):
  33. var = array(varType[0],[-999])
  34. t.Branch(name,var,name+"/"+varType.upper())
  35. return t, var
  36.  
  37. def PutToZero(var):
  38. if not isinstance(var, array):
  39. var.clear()
  40. else:
  41. var[0] = 0
  42.  
  43. def Push(leaf, value):
  44. if not isinstance(leaf, array):
  45. leaf.push_back(value)
  46. else:
  47. leaf[0] = value
  48.  
  49.  
  50. def wasFired(indexKids, detPoints, detPos, pointsVects=None, Ethr=0):
  51. def lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr):
  52. for hit in detPoints:
  53. if (indexKids is None) or (hit.GetTrackID() in indexKids):
  54. #print RPChit.GetTrackID(), t_ID
  55. # check if it is in one of the considered active layer
  56. for pos in detPos:
  57. if pos[0]<=hit.GetZ()<=pos[1]:
  58. Eloss += hit.GetEnergyLoss()
  59. #print pos, hit.GetZ()
  60. if pointsVects is not None:
  61. pointsVects[0].push_back(hit.GetX())
  62. pointsVects[1].push_back(hit.GetY())
  63. pointsVects[2].push_back(hit.GetZ())
  64. flag = True
  65. nstat+=1
  66. #eff_val = gRandom.Uniform(0.,1.)
  67. #if eff_val<0.9:
  68. # flag_eff = flag_eff or True
  69. # nstat_eff+=1
  70. #else:
  71. # flag_eff = flag_eff or False
  72. flag_Ethr = Eloss>=Ethr
  73. return flag, nstat, nstat_eff, Eloss, flag_Ethr
  74. # Now in partKidTrackID I should have the trackID connected to my charged particles
  75. #flag_eff = False
  76. flag = False
  77. nstat = 0
  78. nstat_eff = 0
  79. Eloss = 0
  80. flag,nstat,nstat_eff,Eloss,flag_Ethr = lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr)
  81. #if flag==False and checkOn:
  82. #print "To be study event %s"%entry
  83. #PrintEventPart(particles,pdg)
  84. #PrintRPChits(detPoints,pdg)
  85. #print
  86. #print
  87. #assert(False)
  88. #return flag, flag_eff, nstat, nstat_eff, flag_Ethr
  89. return flag, flag_Ethr
  90.  
  91.  
  92. #ff = ROOT.TFile("histoForWeights.root","read")
  93. #h_GioHans = ff.Get("h_Gio")
  94.  
  95. def calcWeight(x,y,z, E, nodes, entries, weightHist, datatype):
  96. if datatype == 'sig':
  97. return -999
  98. params = {'OPERA': {'material':'Fe','lenght':60.,},
  99. 'window1': {'material':'Fe','lenght':nodes["lidT1O_1"]['z']['dim']},
  100. 'window2': {'material':'Al','lenght':nodes["lidT1I_1"]['z']['dim']},
  101. 'volumeOut': {'material':'Fe','lenght':30.,},
  102. 'volumeIn': {'material':'Al','lenght':8.}}
  103. materialPars = {'Fe':{'A':55.847, #g/mol # fuori
  104. 'rho': 7874 * 1.e+3/1.e+6,#g/cm3
  105. },
  106. 'Al':{'A':26.98, #g/mol # dentro
  107. 'rho': 2700 * 1.e+3/1.e+6 #g/cm3
  108. }}
  109. perc = {'OPERA':0.33,
  110. 'window1':0.015,
  111. 'window2':0.015,
  112. 'volumeOut':0.32,
  113. 'volumeIn':0.32}
  114. intElName = interactionElement(x,y,z, nodes)
  115. NA = 6.022e+23
  116. material = params[intElName]['material']
  117. crossSec = 8.8e-39*E #1.5e-38*fabs(E) ##3./2*(1.2+0.35)*1.e-38*fabs(E)
  118. flux = 2.e+18#5.e+16
  119. L = params[intElName]['lenght']
  120. rho = materialPars[material]['rho']
  121. n_element = perc[intElName]*entries
  122. binN = weightHist.GetXaxis().FindBin(E)
  123. weight = crossSec * NA * rho * L * flux / n_element * weightHist.GetBinContent(binN)
  124. return weight
  125.  
  126. def interactionElement(x,y,z, nodes):
  127. # window1
  128. if nodes["lidT1O_1"]['z']['pos']-nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim']:
  129. return "window1"
  130. # window2
  131. if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT1I_1"]['z']['pos']+nodes["lidT1I_1"]['z']['dim']:
  132. return "window2"
  133. # vacuum tank borders
  134. if nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT6O_1"]['z']['pos']+nodes["lidT6O_1"]['z']['dim']: # include le lid x la parte fuori ma non x la parte dentro
  135. if inInnerLid(x,y,z, nodes):
  136. return "volumeIn"
  137. else:
  138. return "volumeOut"
  139. # opera
  140. try:
  141. minz, maxz = nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'], nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim']
  142. except KeyError:
  143. try:
  144. minz, maxz = nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'], nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim']
  145. except KeyError:
  146. # Pazienza esaurita, ora lo hanno chiamato "volMagneticSpectrometer_1" ed e' lungo piu' di 9m!!!
  147. try:
  148. minz, maxz = nodes["volMagneticSpectrometer_1"]['z']['pos']-nodes["volMagneticSpectrometer_1"]['z']['dim'], nodes["volMagneticSpectrometer_1"]['z']['pos']+nodes["volMagneticSpectrometer_1"]['z']['dim']
  149. except KeyError:
  150. print "Cannot find OPERA"
  151. sys.exit()
  152. #if nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'] < z < nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim']:
  153. #if nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'] < z < nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim']:
  154. if minz < z < maxz:
  155. return "OPERA"
  156.  
  157.  
  158. def inInnerLid(x,y,z, nodes):
  159. if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT6I_1"]['z']['pos']+nodes["lidT6I_1"]['z']['dim']:
  160. if z < nodes["Veto_5"]['z']['pos']-nodes["Veto_5"]['z']['dim']:
  161. a,b = nodes['T1I_1']['x']['dim'], nodes['T1I_1']['y']['dim'] # smaller part of the tank (first 5 m)
  162. else:
  163. a,b = nodes['T2I_1']['x']['dim'], nodes['T2I_1']['y']['dim'] # larger part of the tank
  164. if inEllipse(x,y,a,b):
  165. return True
  166. return False
  167.  
  168.  
  169. def inEllipse(x,y,a,b):
  170. if ( (x*x)/(a*a) + (y*y)/(b*b) ) < 1.:
  171. return True
  172. return False