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