SND@LHC Software
Loading...
Searching...
No Matches
MufluxDigiReco.py
Go to the documentation of this file.
1from __future__ import print_function
2from builtins import range
3import global_variables
4from global_variables import h
5import os
6import ROOT
7import MufluxPatRec
8import shipunit as u
9import rootUtils as ut
10
11from array import array
12
13import matplotlib
14matplotlib.use('pdf')
15import matplotlib.pyplot as plt
16import numpy as np
17
18stop = ROOT.TVector3()
19start = ROOT.TVector3()
20
21
22# function for calculating the strip number from a coordinate, for MuonTagger / RPC
23def StripX(x):
24 # defining constants for rpc properties
25 STRIP_XWIDTH = 0.8625 # internal STRIP V, WIDTH, in cm
26 EXT_STRIP_XWIDTH_L = 0.9625 # nominal (R&L) and Left measured external STRIP V, WIDTH, in cm (beam along z, out from the V plane)
27 EXT_STRIP_XWIDTH_R = 0.86 # measured Right external STRIP V, WIDTH,in cm (beam along z, out from the V plane)
28 V_STRIP_OFF = 0.200
29 NR_VER_STRIPS = 184
30 total_width = (NR_VER_STRIPS - 2) * STRIP_XWIDTH + EXT_STRIP_XWIDTH_L + EXT_STRIP_XWIDTH_R + (NR_VER_STRIPS - 1) * V_STRIP_OFF
31 x_start = (total_width - EXT_STRIP_XWIDTH_R + EXT_STRIP_XWIDTH_L) / 2
32 # calculating strip as an integer
33 strip_x = (x_start - EXT_STRIP_XWIDTH_L + 1.5 * STRIP_XWIDTH + V_STRIP_OFF - x)//(STRIP_XWIDTH + V_STRIP_OFF)
34 if not (0 < strip_x <= NR_VER_STRIPS):
35 print("WARNING: X strip outside range!")
36 strip_x = 0
37 return int(strip_x)
38
39
40def StripY(y):
41 STRIP_YWIDTH = 0.8625 # internal STRIP H, WIDTH, in cm
42 EXT_STRIP_YWIDTH = 0.3 # measured external STRIP H, WIDTH, in cm
43 H_STRIP_OFF = 0.1983
44 NR_HORI_STRIPS = 116
45 total_height = (NR_HORI_STRIPS - 2) * STRIP_YWIDTH + 2 * EXT_STRIP_YWIDTH + (NR_HORI_STRIPS - 1) * H_STRIP_OFF
46 y_start = total_height / 2
47 strip_y = (y_start - EXT_STRIP_YWIDTH + 1.5 * STRIP_YWIDTH + H_STRIP_OFF - y)//(STRIP_YWIDTH + H_STRIP_OFF)
48 if not (0 < strip_y <= NR_HORI_STRIPS):
49 print("WARNING: Y strip outside range!")
50 strip_y = 0
51 return int(strip_y)
52
53class MufluxDigiReco:
54 " convert FairSHiP MC hits / digitized hits to measurements"
55 def __init__(self,fout):
56
57 self.iEvent = 0
58
59 outdir=os.getcwd()
60 outfile=outdir+"/"+fout
61 self.fn = ROOT.TFile(fout,'update')
62 self.sTree = self.fn.cbmsim
63
64 if self.sTree.GetBranch("FitTracks"):
65 print("remove RECO branches and rerun reconstruction")
66 self.fn.Close()
67 # make a new file without reco branches
68 f = ROOT.TFile(fout)
69 sTree = f.cbmsim
70 if sTree.GetBranch("FitTracks"): sTree.SetBranchStatus("FitTracks",0)
71 if sTree.GetBranch("Particles"): sTree.SetBranchStatus("Particles",0)
72 if sTree.GetBranch("fitTrack2MC"): sTree.SetBranchStatus("fitTrack2MC",0)
73 if sTree.GetBranch("FitTracks_PR"): sTree.SetBranchStatus("FitTracks_PR",0)
74 if sTree.GetBranch("Particles_PR"): sTree.SetBranchStatus("Particles_PR",0)
75 if sTree.GetBranch("fitTrack2MC_PR"): sTree.SetBranchStatus("fitTrack2MC_PR",0)
76 if sTree.GetBranch("MufluxSpectrometerPoint"):
77 if sTree.GetBranch("Digi_MufluxSpectrometerHits"):
78 sTree.SetBranchStatus("Digi_MufluxSpectrometerHits",0)
79
80 rawFile = fout.replace("_rec.root","_raw.root")
81 recf = ROOT.TFile(rawFile,"recreate")
82 newTree = sTree.CloneTree(0)
83 for n in range(sTree.GetEntries()):
84 sTree.GetEntry(n)
85 rc = newTree.Fill()
86 sTree.Clear()
87 newTree.AutoSave()
88 f.Close()
89 recf.Close()
90 os.system('cp '+rawFile +' '+fout)
91 self.fn = ROOT.TFile(fout,'update')
92 self.sTree = self.fn.cbmsim
93
94 # prepare for output
95 if simulation:
96 # event header
97 self.header = ROOT.FairEventHeader()
98 self.eventHeader = self.sTree.Branch("ShipEventHeader",self.header,32000,-1)
99 self.fitTrack2MC = ROOT.std.vector('int')()
100 self.mcLink = self.sTree.Branch("fitTrack2MC" + global_variables.realPR, self.fitTrack2MC, 32000, -1)
101 self.digiMufluxSpectrometer = ROOT.TClonesArray("MufluxSpectrometerHit")
102 self.digiMufluxSpectrometerBranch = self.sTree.Branch("Digi_MufluxSpectrometerHits",self.digiMufluxSpectrometer,32000,-1)
103 #muon taggger
104 if self.sTree.GetBranch("MuonTaggerPoint"):
105 self.digiMuonTagger = ROOT.TClonesArray("MuonTaggerHit")
106 self.digiMuonTaggerBranch = self.sTree.Branch("Digi_MuonTagger", self.digiMuonTagger, 32000, -1)
107 # setup random number generator
108 self.random = ROOT.TRandom()
109 ROOT.gRandom.SetSeed()
110 # fitted tracks
111 self.fGenFitArray = ROOT.TClonesArray("genfit::Track")
112 self.fGenFitArray.BypassStreamer(ROOT.kFALSE)
113 self.fitTracks = self.sTree.Branch("FitTracks" + global_variables.realPR, self.fGenFitArray, 32000, -1)
114
115 self.PDG = ROOT.TDatabasePDG.Instance()
116 # for the digitizing and reconstruction step
117 self.v_drift = global_variables.modules["MufluxSpectrometer"].TubeVdrift()
118 self.sigma_spatial = global_variables.modules["MufluxSpectrometer"].TubeSigmaSpatial()
119 self.viewangle = global_variables.modules["MufluxSpectrometer"].ViewAngle()
120
121 # access ShipTree
122 self.sTree.GetEvent(0)
123 self.geoMat = ROOT.genfit.TGeoMaterialInterface()
124 # init geometry and mag. field
125 self.gMan = ROOT.gGeoManager
126 self.bfield = ROOT.genfit.FairShipFields()
127 self.fM = ROOT.genfit.FieldManager.getInstance()
128 self.fM.init(self.bfield)
129
130 ROOT.genfit.MaterialEffects.getInstance().init(self.geoMat)
131
132 # init fitter, to be done before importing shipPatRec
133 self.fitter = ROOT.genfit.DAF()
134 self.fitter.setMaxIterations(20)
135
136 if global_variables.debug:
137
138 # self.fitter.setDebugLvl(1) # produces lot of printout
139 output_dir = 'pics/'
140 if os.path.exists(output_dir):
141 print('The directory already exists. It is OK.')
142 else:
143 os.mkdir(output_dir)
144
145 def reconstruct(self):
146 ntracks = self.findTracks()
147
148 def digitize(self):
149
150 self.sTree.t0 = self.random.Rndm()*1*u.microsecond
151 self.header.SetEventTime( self.sTree.t0 )
152 self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() )
153 self.header.SetMCEntryNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1
154 self.eventHeader.Fill()
155 self.digiMufluxSpectrometer.Delete()
156 self.digitizeMufluxSpectrometer()
157 self.digiMufluxSpectrometerBranch.Fill()
158 #self.digiMuonTagger.Delete() # produces a lot of warnings, rpc station 0
159 #self.digitizeMuonTagger()
160 #self.digiMuonTaggerBranch.Fill()
161
162 def digitizeMuonTagger(self, fake_clustering=True):
163
164 station = 0
165 strip = 0
166 DetectorID = set() # set of detector ids - already deduplicated
167 for MuonTaggerHit in self.sTree.MuonTaggerPoint:
168 # getting rpc nodes, name and matrix
169 rpc_box = self.gMan.FindNode(
170 MuonTaggerHit.GetX(),
171 MuonTaggerHit.GetY(),
172 MuonTaggerHit.GetZ())
173 rpc = rpc_box.GetName()
174 master_matrix = rpc_box.GetMatrix()
175
176 # getting muon box volume (lower level)
177 muon_box = self.gMan.GetTopVolume().FindNode("VMuonBox_1")
178 muonbox_matrix = muon_box.GetMatrix()
179
180 # translation from top to MuonBox_1
181 point = array('d', [
182 MuonTaggerHit.GetX(),
183 MuonTaggerHit.GetY(),
184 MuonTaggerHit.GetZ(), 1])
185 point_muonbox = array('d', [0, 0, 0, 1])
186 muonbox_matrix.MasterToLocal(point, point_muonbox)
187
188 # translation to local frame
189 point_local = array('d', [0, 0, 0, 1])
190 master_matrix.MasterToLocal(point_muonbox, point_local)
191
192 xcoord = point_local[0]
193 ycoord = point_local[1]
194
195 # identify individual rpcs
196 station = int(rpc[-1])
197 if station not in range(1, 6): # limiting the range of rpcs
198 print("WARNING: Invalid RPC number, something's wrong with the geometry ",station)
199
200 # calculate strip
201 # x gives vertical direction
202 direction = 1
203 strip = StripX(xcoord)
204 if not strip:
205 continue
206 # sampling number of strips around the exact strip for emulating clustering
207 if fake_clustering:
208 s = np.random.poisson(3)
209 strip = strip - int(s/2)
210 for i in range(0, s):
211 detectorid = station*10000 + direction*1000 + strip + i
212 DetectorID.add(detectorid)
213 else:
214 detectorid = station*10000 + direction*1000 + strip
215 DetectorID.add(detectorid)
216
217 # y gives horizontal direction
218 direction = 0
219 strip = StripY(ycoord)
220 if not strip:
221 continue
222 # sampling number of strips around the exact strip for emulating clustering
223 if fake_clustering:
224 s = np.random.poisson(3)
225 strip = strip - int(s/2)
226 for i in range(0, s):
227 detectorid = station*10000 + direction*1000 + strip + i
228 DetectorID.add(detectorid)
229 else:
230 detectorid = station*10000 + direction*1000 + strip
231 DetectorID.add(detectorid)
232
233 self.digiMuonTagger.Expand(len(DetectorID))
234 for index, detID in enumerate(DetectorID):
235 hit = ROOT.MuonTaggerHit(detID, 0)
236 self.digiMuonTagger[index] = hit
237
238 if fake_clustering:
239 # cluster size loop - plotting the cluster size distribution
240 cluster_size = list()
241 DetectorID_list = list(DetectorID) # turn set into list to allow indexing
242 DetectorID_list.sort() # sorting the list
243 if len(DetectorID_list) > 1:
244 clusters = [[DetectorID_list[0]]]
245 for x in DetectorID_list[1:]:
246 if abs(x - clusters[-1][-1]) <= 1:
247 clusters[-1].append(x)
248 else:
249 clusters.append([x])
250 cluster_size = [len(x) for x in clusters]
251 for i in cluster_size:
252 h['muontagger_clusters'].Fill(i)
253
254
255 def digitizeMufluxSpectrometer(self):
256
257 # digitize FairSHiP MC hits
258 index = 0
259 hitsPerDetId = {}
260
261 for aMCPoint in self.sTree.MufluxSpectrometerPoint:
262 aHit = ROOT.MufluxSpectrometerHit(aMCPoint,self.sTree.t0)
263 if self.digiMufluxSpectrometer.GetSize() == index: self.digiMufluxSpectrometer.Expand(index+1000)
264 self.digiMufluxSpectrometer[index]=aHit
265 detID = aHit.GetDetectorID()
266 if detID in hitsPerDetId:
267 if self.digiMufluxSpectrometer[hitsPerDetId[detID]].tdc() > aHit.tdc():
268 # second hit with smaller tdc
269 self.digiMufluxSpectrometer[hitsPerDetId[detID]].setInvalid()
270 hitsPerDetId[detID] = index
271 else:
272 hitsPerDetId[detID] = index
273 index+=1
274
275 T1_entries_px = {}
276 T4_entries_px = {}
277 nMufluxHits = self.sTree.MufluxSpectrometerPoint.GetEntriesFast()
278 for i in range(nMufluxHits):
279 MufluxHit = self.sTree.MufluxSpectrometerPoint[i]
280 detector = self.gMan.FindNode(MufluxHit.GetX(),MufluxHit.GetY(),MufluxHit.GetZ()).GetName()
281 MufluxTrackId = MufluxHit.GetTrackID()
282 pid = MufluxHit.PdgCode()
283 xcoord = MufluxHit.GetX()
284 ycoord = MufluxHit.GetY()
285 if abs(pid)==13:
286 if (detector[0:8]=="gas_12_1"):
287 rc=h['hits-T1'].Fill(xcoord,ycoord)
288 if (detector[0:9]=="gas_12_10"):
289 rc=h['hits-T1x'].Fill(xcoord,ycoord)
290 if (detector[0:9]=="gas_12_11"):
291 rc=h['hits-T1u'].Fill(xcoord,ycoord)
292 if (detector[0:8]=="gas_12_2"):
293 rc=h['hits-T2'].Fill(xcoord,ycoord)
294 if (detector[0:9]=="gas_12_20"):
295 rc=h['hits-T2v'].Fill(xcoord,ycoord)
296 if (detector[0:9]=="gas_12_21"):
297 rc=h['hits-T2x'].Fill(xcoord,ycoord)
298 if (detector[0:5]=="gas_3"):
299 rc=h['hits-T3'].Fill(xcoord,ycoord)
300 if (detector[0:5]=="gas_4"):
301 rc=h['hits-T4'].Fill(xcoord,ycoord)
302
303 if (detector[0:9]=="gas_12_10"):
304 if MufluxTrackId in T1_entries_px:
305 continue
306 else:
307 if abs(pid)==13 :
308 T1_entries_px[MufluxTrackId]=[MufluxHit.GetPx()]
309
310 if (detector[0:5]=="gas_4"):
311 if MufluxTrackId in T4_entries_px:
312 continue
313 else:
314 pid = MufluxHit.PdgCode()
315 if abs(pid)==13 :
316 T4_entries_px[MufluxTrackId]=[MufluxHit.GetPx()]
317
318 for i in range(nMufluxHits):
319 MufluxHit = self.sTree.MufluxSpectrometerPoint[i]
320 MufluxTrackId = MufluxHit.GetTrackID()
321 if (T1_entries_px.get(MufluxTrackId) is None or T4_entries_px.get(MufluxTrackId) is None) :
322 continue
323 else:
324 rc=h['pt-kick'].Fill(T1_entries_px.get(MufluxTrackId)[0]-T4_entries_px.get(MufluxTrackId)[0])
325
326 def withT0Estimate(self):
327 # loop over all straw tdcs and make average, correct for ToF
328 n = 0
329 t0 = 0.
330 key = -1
331 SmearedHits = []
332 v_drift = global_variables.modules["MufluxSpectrometer"].TubeVdrift()
333 z1 = stop.z()
334 for aDigi in self.digiMufluxSpectrometer:
335 key+=1
336 if not aDigi.isValid: continue
337 detID = aDigi.GetDetectorID()
338 # don't use hits from straw veto
339 station = int(detID/10000000)
340 if station > 4 : continue
341 aDigi.MufluxSpectrometerEndPoints(start,stop)
342 # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop)
343 delt1 = (start[2]-z1)/u.speedOfLight
344 t0+=aDigi.GetDigi()-delt1
345 SmearedHits.append( {'digiHit':key,'xtop':stop.x(),'ytop':stop.y(),'z':stop.z(),'xbot':start.x(),'ybot':start.y(),'dist':aDigi.GetDigi(), 'detID':detID} )
346 n+=1
347 if n>0:
348 t0 = t0/n - 73.2*u.ns
349 print("t0 ",t0)
350 for s in SmearedHits:
351 delt1 = (s['z']-z1)/u.speedOfLight
352 s['dist'] = (s['dist'] -delt1 -t0)*v_drift
353 print("s['dist']",s['dist'])
354 return SmearedHits
355
356 def smearHits(self,no_amb=None):
357 # smear strawtube points
358 SmearedHits = []
359 key = -1
360 for ahit in self.sTree.MufluxSpectrometerPoint:
361 key+=1
362 detID = ahit.GetDetectorID()
363 top = ROOT.TVector3()
364 bot = ROOT.TVector3()
365 global_variables.modules["MufluxSpectrometer"].TubeEndPoints(detID, bot, top)
366 # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop)
367 #distance to wire, and smear it.
368 dw = ahit.dist2Wire()
369 smear = dw
370 if not no_amb:
371 smear = abs(self.random.Gaus(dw,self.sigma_spatial))
372
373 SmearedHits.append( {'digiHit':key,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(),'dist':smear, 'detID':detID} )
374 # Note: top.z()==bot.z() unless misaligned, so only add key 'z' to smearedHit
375
376 if abs(top.y())==abs(bot.y()): h['disty'].Fill(dw)
377 if abs(top.y())>abs(bot.y()): h['distu'].Fill(dw)
378 if abs(top.y())<abs(bot.y()): h['distv'].Fill(dw)
379
380 return SmearedHits
381
382
383 def stationInfo(self, hit):
384 # Taken from charmdet/drifttubeMonitoring.py
385 detid = hit.GetDetectorID()
386 statnb = detid/10000000;
387 vnb = (detid - statnb*10000000)/1000000
388 pnb = (detid - statnb*10000000 - vnb*1000000)/100000
389 lnb = (detid - statnb*10000000 - vnb*1000000 - pnb*100000)/10000
390 view = "_x"
391 if vnb==0 and statnb==2: view = "_v"
392 if vnb==1 and statnb==1: view = "_u"
393 if pnb>1:
394 print("something wrong with detector id",detid)
395 pnb = 0
396 return statnb,vnb,pnb,lnb,view
397
398 def correctAlignment(self, hit):
399 # Taken from charmdet/drifttubeMonitoring.py
400 sqrt2 = ROOT.TMath.Sqrt(2.)
401 cos30 = ROOT.TMath.Cos(30./180.*ROOT.TMath.Pi())
402 sin30 = ROOT.TMath.Sin(30./180.*ROOT.TMath.Pi())
403 #delX=[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
404 delX=[[-1.2,1.4,-0.9,0.7],[-1.69-1.31,0.3+2.82,0.6-0.4,-0.8-0.34],[-0.9,1.14,-0.80,1.1],[-0.8,1.29,0.15,2.0]]
405 #delUV = [[1.8,-2.,0.8,-1.],[1.2,-1.1,-0.5,0.3]]
406 delUV = [[0,0,0,0],[0,0,0,0]]
407 vtop = ROOT.TVector3()
408 vbot = ROOT.TVector3()
409 rc = hit.MufluxSpectrometerEndPoints(vbot,vtop)
410 s,v,p,l,view = self.stationInfo(hit)
411 if view=='_x':
412 cor = delX[s-1][2*l+p]
413 vbot[0]=vbot[0]+cor
414 vtop[0]=vtop[0]+cor
415 else:
416 if view=='_u': cor = delUV[0][2*l+p]
417 elif view=='_v': cor = delUV[1][2*l+p]
418 vbot[0]=vbot[0]+cor*cos30
419 vtop[0]=vtop[0]+cor*cos30
420 vbot[1]=vbot[1]+cor*sin30
421 vtop[1]=vtop[1]+cor*sin30
422 #tmp = vbot[0]*ROOT.TMath.Cos(rotUV)-vbot[1]*ROOT.TMath.Sin(rotUV)
423 #vbot[1]=vbot[1]*ROOT.TMath.Cos(rotUV)+vbot[0]*ROOT.TMath.Sin(rotUV)
424 #vbot[0]=tmp
425 #tmp = vtop[0]*ROOT.TMath.Cos(rotUV)-vtop[1]*ROOT.TMath.Sin(rotUV)
426 #vtop[1]=vtop[1]*ROOT.TMath.Cos(rotUV)+vtop[0]*ROOT.TMath.Sin(rotUV)
427 #vbot[0]=tmp
428 return vbot,vtop
429
430
431 def RT(self, s,t,function='parabola'):
432 # Taken from charmdet/drifttubeMonitoring.py
433 # rt relation, drift time to distance, drift time?
434 tMinAndTmax = {1:[587,1860],2:[587,1860],3:[610,2300],4:[610,2100]}
435 R = global_variables.ShipGeo.MufluxSpectrometer.InnerTubeDiameter/2. # = 3.63*u.cm
436 # parabola
437 if function == 'parabola' or 'rtTDC'+str(s)+'000_x' not in h:
438 p1p2 = {1:[688.,7.01],2:[688.,7.01],3:[923.,4.41],4:[819.,0.995]}
439 t0_corr = max(0,t-tMinAndTmax[s][0])
440 tmp1 = ROOT.TMath.Sqrt(p1p2[s][0]**2+4.*p1p2[s][1]*t0_corr)
441 r = (2*tmp1-2*p1p2[s][0]+p1p2[s][0]*ROOT.TMath.Log(-p1p2[s][0]+2*tmp1)-p1p2[s][0]*ROOT.TMath.Log(p1p2[s][0]) )/(8*p1p2[s][1])
442 else:
443 if t > tMinAndTmax[s][1]: r = R
444 elif t< tMinAndTmax[s][0]: r = 0
445 else: r = h['rtTDC'+str(s)+'000_x'].Eval(t)
446 # h['TDC2R'].Fill(t,r)
447 return r
448
449
450 def smearHits_realData(self):
451
452 # smear strawtube points
453 SmearedHits = []
454 key = -1
455 for ahit in self.sTree.Digi_MufluxSpectrometerHits:
456 key+=1
457 detID = ahit.GetDetectorID()
458 top = ROOT.TVector3()
459 bot = ROOT.TVector3()
460 bot, top = self.correctAlignment(ahit)
461 # global_variables.modules["MufluxSpectrometer"].TubeEndPoints(detID, bot, top)
462 # ahit.MufluxSpectrometerEndPoints(bot,top)
463 # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop)
464 # distance to wire.
465 # dist = ahit.GetDigi() * 3.7 / (2000. * 2.)
466 tdc = ahit.GetDigi()
467 dist = self.RT(ahit.GetDetectorID()/10000000, tdc)
468
469 SmearedHits.append( {'digiHit':key,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(),'dist':dist, 'detID':detID} )
470
471 if abs(top.y())==abs(bot.y()): h['disty'].Fill(dist)
472 if abs(top.y())>abs(bot.y()): h['distu'].Fill(dist)
473 if abs(top.y())<abs(bot.y()): h['distv'].Fill(dist)
474
475 return SmearedHits
476
477
478 def muonTaggerHits_mc(self):
479
480 # smear strawtube points
481 TaggerHits = []
482 key = -1
483 for ahit in self.sTree.MuonTaggerPoint:
484 key+=1
485 detID = ahit.GetDetectorID()
486
487 TaggerHits.append( {'digiHit':key,'xtop':ahit.GetX(),'ytop':ahit.GetY(),'z':ahit.GetZ(),
488 'xbot':ahit.GetX(),'ybot':ahit.GetY(), 'detID':detID} )
489
490 return TaggerHits
491
492 def muonTaggerHits_realData(self):
493
494 # smear strawtube points
495 TaggerHits = []
496 key = -1
497 for ahit in self.sTree.Digi_MuonTaggerHits:
498 key+=1
499 detID = ahit.GetDetectorID()
500 top = ROOT.TVector3()
501 bot = ROOT.TVector3()
502 ahit.EndPoints(bot,top)
503
504 TaggerHits.append( {'digiHit':key,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(), 'detID':detID} )
505
506 return TaggerHits
507
508
509 def getPtruthFirst(self,mcPartKey):
510 Ptruth,Ptruthx,Ptruthy,Ptruthz = -1.,-1.,-1.,-1.
511 for ahit in self.sTree.MufluxSpectrometerPoint:
512 if ahit.GetTrackID() == mcPartKey:
513 Ptruthx,Ptruthy,Ptruthz = ahit.GetPx(),ahit.GetPy(),ahit.GetPz()
514 Ptruth = ROOT.TMath.Sqrt(Ptruthx**2+Ptruthy**2+Ptruthz**2)
515 break
516 return Ptruth,Ptruthx,Ptruthy,Ptruthz
517
518 def getPtruthAtOrigin(self,mcPartKey):
519 Ptruth,Ptruthx,Ptruthy,Ptruthz = -1.,-1.,-1.,-1.
520 atrack=self.sTree.MCTrack[mcPartKey]
521 Ptruthx= atrack.GetPx()
522 Ptruthy= atrack.GetPy()
523 Ptruthz= atrack.GetPz()
524 Ptruth = ROOT.TMath.Sqrt(Ptruthx**2+Ptruthy**2+Ptruthz**2)
525 return Ptruth,Ptruthx,Ptruthy,Ptruthz
526
527
528 def fracMCsame(self, trackids):
529
530 track = {}
531 nh = len(trackids)
532
533 for tid in trackids:
534 if tid in track:
535 track[tid] += 1
536 else:
537 track[tid] = 1
538
539 # now get track with largest number of hits
540 if track != {}:
541 tmax = max(track, key=track.get)
542 else:
543 track = {-999:0}
544 tmax = -999
545
546 frac = 0.
547 if nh > 0:
548 frac = float(track[tmax]) / float(nh)
549
550 return frac,tmax
551
552
553 def recognition_metrics(self, track_hit_ids, mode='all'):
554
555 n_tracks_y12 = 0
556 n_recognized_y12 = 0
557 n_clones_y12 = 0
558 n_ghosts_y12 = 0
559 n_others_y12 = 0
560
561 n_tracks_stereo12 = 0
562 n_recognized_stereo12 = 0
563 n_clones_stereo12 = 0
564 n_ghosts_stereo12 = 0
565 n_others_stereo12 = 0
566
567 n_tracks_34 = 0
568 n_recognized_34 = 0
569 n_clones_34 = 0
570 n_ghosts_34 = 0
571 n_others_34 = 0
572
573 true_track_ids_y12 = []
574 true_track_ids_stereo12 = []
575 true_track_ids_34 = []
576
577 n_true_track_hits_y12 = {}
578 n_true_track_hits_stereo12 = {}
579 n_true_track_hits_34 = {}
580
581 n_true_track_hits_1 = {}
582 n_true_track_hits_2 = {}
583 n_true_track_hits_3 = {}
584 n_true_track_hits_4 = {}
585
586 for ahit in self.sTree.MufluxSpectrometerPoint:
587
588 track_id = ahit.GetTrackID()
589
590 detID = ahit.GetDetectorID()
591 statnb = detID // 10000000
592 vnb = (detID - statnb * 10000000) // 1000000
593
594 is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1)
595 is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0)
596 is_34 = (statnb == 3) + (statnb == 4)
597
598 if statnb == 1:
599 if track_id in n_true_track_hits_1:
600 n_true_track_hits_1[track_id] += 1
601 else:
602 n_true_track_hits_1[track_id] = 1
603
604 if statnb == 2:
605 if track_id in n_true_track_hits_2:
606 n_true_track_hits_2[track_id] += 1
607 else:
608 n_true_track_hits_2[track_id] = 1
609
610 if statnb == 3:
611 if track_id in n_true_track_hits_3:
612 n_true_track_hits_3[track_id] += 1
613 else:
614 n_true_track_hits_3[track_id] = 1
615
616 if statnb == 4:
617 if track_id in n_true_track_hits_4:
618 n_true_track_hits_4[track_id] += 1
619 else:
620 n_true_track_hits_4[track_id] = 1
621
622 if is_y12:
623 if track_id in n_true_track_hits_y12:
624 n_true_track_hits_y12[track_id] += 1
625 else:
626 n_true_track_hits_y12[track_id] = 1
627
628 if is_stereo12:
629 if track_id in n_true_track_hits_stereo12:
630 n_true_track_hits_stereo12[track_id] += 1
631 else:
632 n_true_track_hits_stereo12[track_id] = 1
633
634 if is_34:
635 if track_id in n_true_track_hits_34:
636 n_true_track_hits_34[track_id] += 1
637 else:
638 n_true_track_hits_34[track_id] = 1
639
640
641 if mode == 'all':
642 min_hits = 1
643 elif mode == '3hits':
644 min_hits = 3
645
646 if mode == 'all' or mode == '3hits':
647 for key in n_true_track_hits_y12.keys():
648 if n_true_track_hits_y12[key] >= min_hits:
649 true_track_ids_y12.append(key)
650
651 for key in n_true_track_hits_stereo12.keys():
652 if n_true_track_hits_stereo12[key] >= min_hits:
653 true_track_ids_stereo12.append(key)
654
655 for key in n_true_track_hits_34.keys():
656 if n_true_track_hits_34[key] >= min_hits:
657 true_track_ids_34.append(key)
658
659 if mode == 'Tr4':
660 min_hits = 1
661 for key in n_true_track_hits_y12.keys():
662 if key in n_true_track_hits_1 and key in n_true_track_hits_2:
663 if key in n_true_track_hits_3 and key in n_true_track_hits_4:
664 if n_true_track_hits_y12[key] >= min_hits:
665 true_track_ids_y12.append(key)
666
667 for key in n_true_track_hits_stereo12.keys():
668 if key in n_true_track_hits_1 and key in n_true_track_hits_2:
669 if key in n_true_track_hits_3 and key in n_true_track_hits_4:
670 if n_true_track_hits_stereo12[key] >= min_hits:
671 true_track_ids_stereo12.append(key)
672
673 for key in n_true_track_hits_34.keys():
674 if key in n_true_track_hits_1 and key in n_true_track_hits_2:
675 if key in n_true_track_hits_3 and key in n_true_track_hits_4:
676 if n_true_track_hits_34[key] >= min_hits:
677 true_track_ids_34.append(key)
678
679
680
681 n_tracks_y12 = len(true_track_ids_y12)
682 n_tracks_stereo12 = len(true_track_ids_stereo12)
683 n_tracks_34 = len(true_track_ids_34)
684
685 found_track_ids_y12 = []
686 found_track_ids_stereo12 = []
687 found_track_ids_34 = []
688
689 for i_track in track_hit_ids.keys():
690
691 atrack = track_hit_ids[i_track]
692 atrack_y12 = atrack['y12']
693 atrack_stereo12 = atrack['stereo12']
694 atrack_34 = atrack['34']
695 atrack_y_tagger = atrack['y_tagger']
696
697 if len(atrack_y12) > 0:
698 reco_hit_ids_y12 = []
699 for i_hit in atrack_y12:
700 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
701 reco_hit_ids_y12.append(ahit.GetTrackID())
702 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
703 if frac_y12 >= 0.7:
704 if tmax_y12 in true_track_ids_y12 and tmax_y12 not in found_track_ids_y12:
705 n_recognized_y12 += 1
706 found_track_ids_y12.append(tmax_y12)
707 elif tmax_y12 in true_track_ids_y12 and tmax_y12 in found_track_ids_y12:
708 n_clones_y12 += 1
709 elif tmax_y12 not in true_track_ids_y12:
710 n_others_y12 += 1
711 else:
712 n_ghosts_y12 += 1
713
714
715 if len(atrack_stereo12) > 0:
716 reco_hit_ids_stereo12 = []
717 for i_hit in atrack_stereo12:
718 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
719 reco_hit_ids_stereo12.append(ahit.GetTrackID())
720 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
721 if frac_stereo12 >= 0.7:
722 if tmax_stereo12 in true_track_ids_stereo12 and tmax_stereo12 not in found_track_ids_stereo12:
723 n_recognized_stereo12 += 1
724 found_track_ids_stereo12.append(tmax_stereo12)
725 elif tmax_stereo12 in true_track_ids_stereo12 and tmax_stereo12 in found_track_ids_stereo12:
726 n_clones_stereo12 += 1
727 elif tmax_stereo12 not in true_track_ids_stereo12:
728 n_others_stereo12 += 1
729 else:
730 n_ghosts_stereo12 += 1
731
732
733 if len(atrack_34) > 0:
734 reco_hit_ids_34 = []
735 for i_hit in atrack_34:
736 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
737 reco_hit_ids_34.append(ahit.GetTrackID())
738 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
739 if frac_34 >= 0.7:
740 if tmax_34 in true_track_ids_34 and tmax_34 not in found_track_ids_34:
741 n_recognized_34 += 1
742 found_track_ids_34.append(tmax_34)
743 elif tmax_34 in true_track_ids_34 and tmax_34 in found_track_ids_34:
744 n_clones_34 += 1
745 elif tmax_34 not in true_track_ids_34:
746 n_others_34 += 1
747 else:
748 n_ghosts_34 += 1
749
750 output = (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12,
751 n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12,
752 n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34)
753
754 return output
755
756
757 def target_metrics(self, track_hit_ids):
758
759 n_tracks = 0
760 n_recognized = 0
761 n_clones = 0
762 n_ghosts = 0
763 n_others = 0
764
765 true_track_ids = []
766
767 n_true_track_hits_y12 = {}
768 n_true_track_hits_stereo12 = {}
769 n_true_track_hits_34 = {}
770
771 for ahit in self.sTree.MufluxSpectrometerPoint:
772
773 track_id = ahit.GetTrackID()
774
775 detID = ahit.GetDetectorID()
776 statnb = detID // 10000000
777 vnb = (detID - statnb * 10000000) // 1000000
778
779 is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1)
780 is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0)
781 is_34 = (statnb == 3) + (statnb == 4)
782
783 if is_y12:
784 if track_id in n_true_track_hits_y12:
785 n_true_track_hits_y12[track_id] += 1
786 else:
787 n_true_track_hits_y12[track_id] = 1
788
789 if is_stereo12:
790 if track_id in n_true_track_hits_stereo12:
791 n_true_track_hits_stereo12[track_id] += 1
792 else:
793 n_true_track_hits_stereo12[track_id] = 1
794
795 if is_34:
796 if track_id in n_true_track_hits_34:
797 n_true_track_hits_34[track_id] += 1
798 else:
799 n_true_track_hits_34[track_id] = 1
800
801 min_hits = 3
802 for key in n_true_track_hits_y12.keys():
803 if n_true_track_hits_y12[key] >= min_hits:
804 if key in n_true_track_hits_stereo12:
805 if n_true_track_hits_stereo12[key] >= min_hits:
806 if key in n_true_track_hits_34:
807 if n_true_track_hits_34[key] >= min_hits:
808 true_track_ids.append(key)
809
810 n_tracks = len(true_track_ids)
811
812 found_track_ids = []
813
814 for i_track in track_hit_ids.keys():
815
816 atrack = track_hit_ids[i_track]
817 atrack_y12 = atrack['y12']
818 atrack_stereo12 = atrack['stereo12']
819 atrack_34 = atrack['34']
820 atrack_y_tagger = atrack['y_tagger']
821
822 if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits:
823
824 reco_hit_ids_y12 = []
825 for i_hit in atrack_y12:
826 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
827 reco_hit_ids_y12.append(ahit.GetTrackID())
828 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
829
830 reco_hit_ids_stereo12 = []
831 for i_hit in atrack_stereo12:
832 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
833 reco_hit_ids_stereo12.append(ahit.GetTrackID())
834 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
835
836 reco_hit_ids_34 = []
837 for i_hit in atrack_34:
838 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
839 reco_hit_ids_34.append(ahit.GetTrackID())
840 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
841
842
843 min_eff = 0.7
844 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34:
845 if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff:
846
847 if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids:
848 n_recognized += 1
849 found_track_ids.append(tmax_y12)
850 elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids:
851 n_clones += 1
852 elif tmax_y12 not in true_track_ids:
853 n_others += 1
854
855 else:
856 n_ghosts += 1
857 else:
858 n_ghosts += 1
859
860 else:
861 n_ghosts += 1
862
863 output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others)
864
865 return output
866
867
868 def all_tracks_metrics(self, track_hit_ids):
869
870 n_tracks = 0
871 n_recognized = 0
872 n_clones = 0
873 n_ghosts = 0
874 n_others = 0
875
876 true_track_ids = []
877 for ahit in self.sTree.MufluxSpectrometerPoint:
878 track_id = ahit.GetTrackID()
879 if track_id not in true_track_ids:
880 true_track_ids.append(track_id)
881
882 n_tracks = len(true_track_ids)
883 found_track_ids = []
884 min_hits = 3
885
886 for i_track in track_hit_ids.keys():
887
888 atrack = track_hit_ids[i_track]
889 atrack_y12 = atrack['y12']
890 atrack_stereo12 = atrack['stereo12']
891 atrack_34 = atrack['34']
892 atrack_y_tagger = atrack['y_tagger']
893
894 if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits:
895
896 reco_hit_ids_y12 = []
897 for i_hit in atrack_y12:
898 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
899 reco_hit_ids_y12.append(ahit.GetTrackID())
900 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
901
902 reco_hit_ids_stereo12 = []
903 for i_hit in atrack_stereo12:
904 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
905 reco_hit_ids_stereo12.append(ahit.GetTrackID())
906 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
907
908 reco_hit_ids_34 = []
909 for i_hit in atrack_34:
910 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
911 reco_hit_ids_34.append(ahit.GetTrackID())
912 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
913
914 min_eff = 0.7
915 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34:
916 if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff:
917
918 if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids:
919 n_recognized += 1
920 found_track_ids.append(tmax_y12)
921 elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids:
922 n_clones += 1
923 elif tmax_y12 not in true_track_ids:
924 n_others += 1
925
926 else:
927 n_ghosts += 1
928 else:
929 n_ghosts += 1
930
931 else:
932 n_ghosts += 1
933
934 output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others)
935
936 return output
937
938
939 def muon_metrics(self, track_hit_ids):
940
941 n_tracks = 0
942 n_recognized = 0
943 n_clones = 0
944 n_ghosts = 0
945 n_others = 0
946
947 true_track_ids = []
948 for ahit in self.sTree.MufluxSpectrometerPoint:
949 track_id = ahit.GetTrackID()
950 pdg = ahit.PdgCode()
951 if track_id not in true_track_ids:
952 if abs(pdg)==13:
953 true_track_ids.append(track_id)
954
955 n_tracks = len(true_track_ids)
956 found_track_ids = []
957 min_hits = 3
958
959 for i_track in track_hit_ids.keys():
960
961 atrack = track_hit_ids[i_track]
962 atrack_y12 = atrack['y12']
963 atrack_stereo12 = atrack['stereo12']
964 atrack_34 = atrack['34']
965 atrack_y_tagger = atrack['y_tagger']
966
967 if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits:
968
969 reco_hit_ids_y12 = []
970 for i_hit in atrack_y12:
971 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
972 reco_hit_ids_y12.append(ahit.GetTrackID())
973 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
974
975 reco_hit_ids_stereo12 = []
976 for i_hit in atrack_stereo12:
977 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
978 reco_hit_ids_stereo12.append(ahit.GetTrackID())
979 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
980
981 reco_hit_ids_34 = []
982 for i_hit in atrack_34:
983 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
984 reco_hit_ids_34.append(ahit.GetTrackID())
985 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
986
987 min_eff = 0.7
988 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34:
989 if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff:
990
991 if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids:
992 n_recognized += 1
993 found_track_ids.append(tmax_y12)
994 elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids:
995 n_clones += 1
996 elif tmax_y12 not in true_track_ids:
997 n_others += 1
998
999 else:
1000 n_ghosts += 1
1001 else:
1002 n_ghosts += 1
1003
1004 else:
1005 n_ghosts += 1
1006
1007 output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others)
1008
1009 return output
1010
1011
1012 def muon_with_tagger_metrics(self, track_hit_ids):
1013
1014 n_tracks = 0
1015 n_recognized = 0
1016 n_clones = 0
1017 n_ghosts = 0
1018 n_others = 0
1019
1020 true_tagger_track_ids = []
1021 for ahit in self.sTree.MuonTaggerPoint:
1022 track_id = ahit.GetTrackID()
1023 pdg = ahit.PdgCode()
1024 if track_id not in true_tagger_track_ids:
1025 if abs(pdg)==13:
1026 true_tagger_track_ids.append(track_id)
1027
1028 true_track_ids = []
1029 for ahit in self.sTree.MufluxSpectrometerPoint:
1030 track_id = ahit.GetTrackID()
1031 pdg = ahit.PdgCode()
1032 if track_id not in true_track_ids:
1033 if abs(pdg)==13:
1034 if track_id in true_tagger_track_ids:
1035 true_track_ids.append(track_id)
1036
1037 n_tracks = len(true_track_ids)
1038 found_track_ids = []
1039 min_hits = 3
1040
1041 for i_track in track_hit_ids.keys():
1042
1043 atrack = track_hit_ids[i_track]
1044 atrack_y12 = atrack['y12']
1045 atrack_stereo12 = atrack['stereo12']
1046 atrack_34 = atrack['34']
1047 atrack_y_tagger = atrack['y_tagger']
1048
1049 if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits:
1050
1051 reco_hit_ids_y12 = []
1052 for i_hit in atrack_y12:
1053 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1054 reco_hit_ids_y12.append(ahit.GetTrackID())
1055 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
1056
1057 reco_hit_ids_stereo12 = []
1058 for i_hit in atrack_stereo12:
1059 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1060 reco_hit_ids_stereo12.append(ahit.GetTrackID())
1061 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
1062
1063 reco_hit_ids_34 = []
1064 for i_hit in atrack_34:
1065 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1066 reco_hit_ids_34.append(ahit.GetTrackID())
1067 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
1068
1069 min_eff = 0.7
1070 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34:
1071 if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff:
1072
1073 if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids:
1074 n_recognized += 1
1075 found_track_ids.append(tmax_y12)
1076 elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids:
1077 n_clones += 1
1078 elif tmax_y12 not in true_track_ids:
1079 n_others += 1
1080
1081 else:
1082 n_ghosts += 1
1083 else:
1084 n_ghosts += 1
1085
1086 else:
1087 n_ghosts += 1
1088
1089 output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others)
1090
1091 return output
1092
1093
1094 def muon_metrics_vs_p(self, track_hit_ids):
1095
1096 n_tracks = 0
1097 n_recognized = 0
1098 n_clones = 0
1099 n_ghosts = 0
1100 n_others = 0
1101
1102 true_track_ids = []
1103 true_track_p = {}
1104 for ahit in self.sTree.MufluxSpectrometerPoint:
1105 track_id = ahit.GetTrackID()
1106 pdg = ahit.PdgCode()
1107 if track_id not in true_track_ids:
1108 if abs(pdg)==13:
1109 true_track_ids.append(track_id)
1110 if track_id not in true_track_p:
1111 Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(track_id)
1112 true_track_p[track_id] = Ptruth
1113 h['True_all_tracks_vs_p_true'].Fill(Ptruth)
1114 if abs(pdg)==13:
1115 h['True_muons_vs_p_true'].Fill(Ptruth)
1116
1117 n_tracks = len(true_track_ids)
1118 found_track_ids = []
1119 min_hits = 3
1120
1121 for i_track in track_hit_ids.keys():
1122
1123 atrack = track_hit_ids[i_track]
1124 atrack_y12 = atrack['y12']
1125 atrack_stereo12 = atrack['stereo12']
1126 atrack_34 = atrack['34']
1127 atrack_y_tagger = atrack['y_tagger']
1128
1129 reco_hit_ids_y12 = []
1130 for i_hit in atrack_y12:
1131 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1132 reco_hit_ids_y12.append(ahit.GetTrackID())
1133 frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12)
1134
1135 reco_hit_ids_stereo12 = []
1136 for i_hit in atrack_stereo12:
1137 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1138 reco_hit_ids_stereo12.append(ahit.GetTrackID())
1139 frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12)
1140
1141 reco_hit_ids_34 = []
1142 for i_hit in atrack_34:
1143 ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']]
1144 reco_hit_ids_34.append(ahit.GetTrackID())
1145 frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34)
1146
1147 if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits:
1148
1149 min_eff = 0.7
1150 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34:
1151 if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff:
1152
1153 if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids:
1154 n_recognized += 1
1155 found_track_ids.append(tmax_y12)
1156 h['Reco_muons_vs_p_true'].Fill(true_track_p[tmax_y12])
1157 elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids:
1158 n_clones += 1
1159 elif tmax_y12 not in true_track_ids:
1160 n_others += 1
1161
1162 else:
1163 n_ghosts += 1
1164 h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12])
1165 else:
1166 n_ghosts += 1
1167 h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12])
1168
1169 else:
1170 n_ghosts += 1
1171 h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12])
1172
1173 output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others)
1174
1175 return output
1176
1177
1178
1179 def findTracks(self):
1180
1181 hitPosLists = {}
1182 hitPosLists_noT4 = {}
1183 stationCrossed = {}
1184 stationCrossed_noT4 = {}
1185 trackDigiHits = {}
1186 trackDigiHits_noT4 = {}
1187 trackMomentums = {}
1188 trackMomentums_noT4 = {}
1189 trackCandidates = []
1190 trackCandidates_noT4 = []
1191 trackCandidates_all = []
1192 nTrack = -1
1193
1194 self.fGenFitArray.Delete()
1195 self.fitTrack2MC.clear()
1196
1197 # hit smearing
1198 if self.sTree.GetBranch("MufluxSpectrometerPoint"):
1199 if global_variables.withT0:
1200 self.SmearedHits = self.withT0Estimate()
1201 self.TaggerHits = self.muonTaggerHits_mc()
1202 else:
1203 self.SmearedHits = self.smearHits(global_variables.withNoStrawSmearing)
1204 self.TaggerHits = self.muonTaggerHits_mc()
1205 elif self.sTree.GetBranch("Digi_MufluxSpectrometerHits"):
1206 self.SmearedHits = self.smearHits_realData()
1207 self.TaggerHits = self.muonTaggerHits_realData()
1208 else:
1209 print('No branches "MufluxSpectrometer" or "Digi_MufluxSpectrometerHits".')
1210 return nTrack
1211
1212 if global_variables.realPR:
1213
1214 # Do real PatRec
1215 track_hits = MufluxPatRec.execute(self.SmearedHits, self.TaggerHits, withNTaggerHits, withDist2Wire)
1216
1217 # Create hitPosLists for track fit
1218 for i_track in track_hits.keys():
1219
1220 atrack = track_hits[i_track]
1221 atrack_y12 = atrack['y12']
1222 atrack_stereo12 = atrack['stereo12']
1223 atrack_34 = atrack['34']
1224 atrack_smeared_hits = list(atrack_y12) + list(atrack_stereo12) + list(atrack_34)
1225 atrack_p = np.abs(atrack['p'])
1226 atrack_y_in_magnet = atrack['y_in_magnet']
1227 atrack_x_in_magnet = atrack['x_in_magnet']
1228
1229 for sm in atrack_smeared_hits:
1230
1231 # detID = self.digiMufluxSpectrometer[sm['digiHit']].GetDetectorID()
1232 detID = sm['detID']
1233 station = int(detID/10000000)
1234 trID = i_track
1235
1236 # T1-4 for track fit
1237 if trID not in hitPosLists:
1238 hitPosLists[trID] = ROOT.std.vector('TVectorD')()
1239 stationCrossed[trID] = {}
1240 trackDigiHits[trID] = []
1241 trackMomentums[trID] = atrack_p
1242 m = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']])
1243 hitPosLists[trID].push_back(ROOT.TVectorD(7,m))
1244 if station not in stationCrossed[trID]:
1245 stationCrossed[trID][station]=0
1246 stationCrossed[trID][station]+=1
1247 trackDigiHits[trID].append(sm['digiHit'])
1248
1249 # T1-3 for track fit
1250 if (int(detID/1000000)!=40):
1251 if trID not in hitPosLists_noT4:
1252 hitPosLists_noT4[trID] = ROOT.std.vector('TVectorD')()
1253 stationCrossed_noT4[trID] = {}
1254 trackDigiHits_noT4[trID] = []
1255 trackMomentums_noT4[trID] = atrack_p
1256 m_noT4 = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']])
1257 hitPosLists_noT4[trID].push_back(ROOT.TVectorD(7,m_noT4))
1258 if station not in stationCrossed_noT4[trID]:
1259 stationCrossed_noT4[trID][station]=0
1260 stationCrossed_noT4[trID][station]+=1
1261 trackDigiHits_noT4[trID].append(sm['digiHit'])
1262
1263 # Do fake PatRec
1264 else:
1265
1266 if not self.sTree.GetBranch("MufluxSpectrometerPoint"):
1267 print('Fake PatRec is impossible. No "MufluxSpectrometerPoint" branch.')
1268 return nTrack
1269
1270 track_hits = {}
1271 for sm in self.SmearedHits:
1272
1273 # detID = self.digiMufluxSpectrometer[sm['digiHit']].GetDetectorID()
1274 detID = sm['detID']
1275 station = int(detID/10000000)
1276 vnb = (detID - station * 10000000) // 1000000
1277 is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1)
1278 is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0)
1279 is_34 = (station == 3) + (station == 4)
1280
1281 trID = self.sTree.MufluxSpectrometerPoint[sm['digiHit']].GetTrackID()
1282
1283 # PatRec
1284 if trID not in track_hits:
1285 atrack = {'y12': [], 'stereo12': [], '34': []}
1286 track_hits[trID] = atrack
1287 if is_y12:
1288 track_hits[trID]['y12'].append(sm)
1289 if is_stereo12:
1290 track_hits[trID]['stereo12'].append(sm)
1291 if is_34:
1292 track_hits[trID]['34'].append(sm)
1293
1294 # T1-4 for track fit
1295 if trID not in hitPosLists:
1296 hitPosLists[trID] = ROOT.std.vector('TVectorD')()
1297 stationCrossed[trID] = {}
1298 trackDigiHits[trID] = []
1299 trackMomentums[trID] = 3.
1300 m = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']])
1301 hitPosLists[trID].push_back(ROOT.TVectorD(7,m))
1302 if station not in stationCrossed[trID]:
1303 stationCrossed[trID][station]=0
1304 stationCrossed[trID][station]+=1
1305 trackDigiHits[trID].append(sm['digiHit'])
1306
1307 # T1-3 for track fit
1308 if (int(detID/1000000)!=40):
1309 if trID not in hitPosLists_noT4:
1310 hitPosLists_noT4[trID] = ROOT.std.vector('TVectorD')()
1311 stationCrossed_noT4[trID] = {}
1312 trackDigiHits_noT4[trID] = []
1313 trackMomentums_noT4[trID] = 3.
1314 m_noT4 = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']])
1315 hitPosLists_noT4[trID].push_back(ROOT.TVectorD(7,m_noT4))
1316 if station not in stationCrossed_noT4[trID]:
1317 stationCrossed_noT4[trID][station]=0
1318 stationCrossed_noT4[trID][station]+=1
1319 trackDigiHits_noT4[trID].append(sm['digiHit'])
1320
1321
1322 # All tracks fit
1323 for atrack in hitPosLists:
1324
1325 if atrack < 0:
1326 continue # these are hits not assigned to MC track because low E cut
1327 pdg = 13 # only muons
1328 if not abs(pdg)==13:
1329 continue # only keep muons
1330
1331 meas = hitPosLists[atrack]
1332 mom_init = trackMomentums[atrack]
1333 nM = meas.size()
1334
1335 #if nM < 12 : continue # not enough hits to make a good trackfit
1336 #comment for hits in t1-3
1337 #if len(stationCrossed[atrack]) < 4 :
1338 # continue # not enough stations crossed to make a good trackfit
1339
1340 charge = self.PDG.GetParticle(pdg).Charge()/(3.)
1341 posM = ROOT.TVector3(0, 0, 0)
1342 momM = ROOT.TVector3(0,0,mom_init * u.GeV)
1343 # approximate covariance
1344 covM = ROOT.TMatrixDSym(6)
1345 resolution = self.sigma_spatial
1346 if global_variables.withT0:
1347 resolution = resolution*1.4 # worse resolution due to t0 estimate
1348 for i in range(3):
1349 covM[i][i] = resolution*resolution
1350 covM[0][0]=resolution*resolution*100.
1351 for i in range(3,6):
1352 covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2)
1353 # trackrep
1354 rep = ROOT.genfit.RKTrackRep(pdg)
1355 # smeared start state
1356 stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep)
1357 rep.setPosMomCov(stateSmeared, posM, momM, covM)
1358 # create track
1359 seedState = ROOT.TVectorD(6)
1360 seedCov = ROOT.TMatrixDSym(6)
1361 rep.get6DStateCov(stateSmeared, seedState, seedCov)
1362 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
1363 hitCov = ROOT.TMatrixDSym(7)
1364 hitCov[6][6] = resolution*resolution
1365 for m in meas:
1366 tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to
1367 measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to
1368 measurement.setMaxDistance(1.85*u.cm)
1369 # measurement.setLeftRightResolution(-1)
1370 tp.addRawMeasurement(measurement) # package measurement in the TrackPoint
1371 theTrack.insertPoint(tp) # add point to Track
1372 trackCandidates_all.append([theTrack,atrack])
1373
1374
1375 # T1-4 track fit
1376 for atrack in hitPosLists:
1377
1378 if atrack < 0:
1379 continue # these are hits not assigned to MC track because low E cut
1380 pdg = 13 # only muons
1381 if not abs(pdg)==13:
1382 continue # only keep muons
1383
1384 meas = hitPosLists[atrack]
1385 mom_init = trackMomentums[atrack]
1386 nM = meas.size()
1387
1388 #if nM < 12 : continue # not enough hits to make a good trackfit
1389 #comment for hits in t1-3
1390 if len(stationCrossed[atrack]) < 4 :
1391 continue # not enough stations crossed to make a good trackfit
1392
1393 charge = self.PDG.GetParticle(pdg).Charge()/(3.)
1394 posM = ROOT.TVector3(0, 0, 0)
1395 momM = ROOT.TVector3(0,0,mom_init * u.GeV)
1396 # approximate covariance
1397 covM = ROOT.TMatrixDSym(6)
1398 resolution = self.sigma_spatial
1399 if global_variables.withT0:
1400 resolution = resolution*1.4 # worse resolution due to t0 estimate
1401 for i in range(3):
1402 covM[i][i] = resolution*resolution
1403 covM[0][0]=resolution*resolution*100.
1404 for i in range(3,6):
1405 covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2)
1406 # trackrep
1407 rep = ROOT.genfit.RKTrackRep(pdg)
1408 # smeared start state
1409 stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep)
1410 rep.setPosMomCov(stateSmeared, posM, momM, covM)
1411 # create track
1412 seedState = ROOT.TVectorD(6)
1413 seedCov = ROOT.TMatrixDSym(6)
1414 rep.get6DStateCov(stateSmeared, seedState, seedCov)
1415 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
1416 hitCov = ROOT.TMatrixDSym(7)
1417 hitCov[6][6] = resolution*resolution
1418 for m in meas:
1419 tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to
1420 measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to
1421 measurement.setMaxDistance(1.85*u.cm)
1422 # measurement.setLeftRightResolution(-1)
1423 tp.addRawMeasurement(measurement) # package measurement in the TrackPoint
1424 theTrack.insertPoint(tp) # add point to Track
1425 trackCandidates.append([theTrack,atrack])
1426
1427 # T1-3 track fit
1428 for atrack in hitPosLists_noT4:
1429
1430 if atrack < 0:
1431 continue # these are hits not assigned to MC track because low E cut
1432 pdg = 13 # only muons
1433 if not abs(pdg)==13:
1434 continue # only keep muons
1435 meas = hitPosLists_noT4[atrack]
1436 mom_init = trackMomentums_noT4[atrack]
1437 nM = meas.size()
1438
1439 #if nM < 6 : continue # not enough hits to make a good trackfit
1440 #comment for hits in t1-3
1441 if len(stationCrossed[atrack]) < 3 :
1442 continue # not enough stations crossed to make a good trackfit
1443
1444 charge = self.PDG.GetParticle(pdg).Charge()/(3.)
1445 posM = ROOT.TVector3(0, 0, 0)
1446 momM = ROOT.TVector3(0,0,mom_init * u.GeV)
1447 # approximate covariance
1448 covM = ROOT.TMatrixDSym(6)
1449 resolution = self.sigma_spatial
1450 if global_variables.withT0:
1451 resolution = resolution*1.4 # worse resolution due to t0 estimate
1452 for i in range(3):
1453 covM[i][i] = resolution*resolution
1454 covM[0][0]=resolution*resolution*100.
1455 for i in range(3,6):
1456 covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2)
1457 # trackrep
1458 rep = ROOT.genfit.RKTrackRep(pdg)
1459 # smeared start state
1460 stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep)
1461 rep.setPosMomCov(stateSmeared, posM, momM, covM)
1462 # create track
1463 seedState = ROOT.TVectorD(6)
1464 seedCov = ROOT.TMatrixDSym(6)
1465 rep.get6DStateCov(stateSmeared, seedState, seedCov)
1466 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
1467 hitCov = ROOT.TMatrixDSym(7)
1468 hitCov[6][6] = resolution*resolution
1469 for m in meas:
1470 tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to
1471 measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to
1472 measurement.setMaxDistance(1.85*u.cm)
1473 # measurement.setLeftRightResolution(-1)
1474 tp.addRawMeasurement(measurement) # package measurement in the TrackPoint
1475 theTrack.insertPoint(tp) # add point to Track
1476 trackCandidates_noT4.append([theTrack,atrack])
1477
1478 if global_variables.debug:
1479
1480 z_hits_y = []
1481 x_hits_y = []
1482 for ahit in self.SmearedHits:
1483 detID = ahit['detID']
1484 station = int(detID/10000000)
1485 vnb = (detID - station * 10000000) // 1000000
1486 is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1)
1487 is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0)
1488 is_34 = (station == 3) + (station == 4)
1489 if is_y12 or is_34:
1490 z_hits_y.append(ahit['z'])
1491 x_hits_y.append(ahit['xtop'])
1492 plt.figure(figsize=(22, 6))
1493 plt.subplot(1, 2, 1)
1494 plt.scatter(z_hits_y, x_hits_y, color='b')
1495
1496 z_tagger_hits_y = []
1497 x_tagger_hits_y = []
1498 for ahit in self.TaggerHits:
1499 detID = ahit['detID']
1500 station = detID // 10000
1501 direction = (detID - station * 10000) // 1000
1502 if direction == 1 or detID < 10:
1503 z_tagger_hits_y.append(ahit['z'])
1504 x_tagger_hits_y.append((ahit['xtop'] + ahit['xbot']) / 2.)
1505 plt.scatter(z_tagger_hits_y, x_tagger_hits_y, color='b')
1506
1507 for i_track in track_hits:
1508 atrack = track_hits[i_track]
1509 atrack_y12 = atrack['y12']
1510 atrack_stereo12 = atrack['stereo12']
1511 atrack_34 = atrack['34']
1512 atrack_tagger = atrack['y_tagger']
1513 atrack_y = list(atrack_y12) + list(atrack_34) + list(atrack_tagger)
1514 z_atrack_y = [ahit['z'] for ahit in atrack_y]
1515 x_atrack_y = [ahit['xtop'] for ahit in atrack_y]
1516 plt.scatter(z_atrack_y, x_atrack_y, color=np.random.rand(3,))
1517
1518 plt.xlabel('z')
1519 plt.ylabel('x')
1520 plt.xlim(0, 1200)
1521 plt.ylim(-150, 150)
1522
1523
1524 plt.subplot(1, 2, 2)
1525 for ahit in self.SmearedHits:
1526 detID = ahit['detID']
1527 station = int(detID/10000000)
1528 vnb = (detID - station * 10000000) // 1000000
1529 is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1)
1530 is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0)
1531 is_34 = (station == 3) + (station == 4)
1532 if is_y12 or is_stereo12:
1533 xbot, xtop = ahit['xbot'], ahit['xtop']
1534 ybot, ytop = ahit['ybot'], ahit['ytop']
1535 plt.plot([ybot, ytop], [xbot, xtop], color='b')
1536
1537 for i_track in track_hits:
1538 atrack = track_hits[i_track]
1539 atrack_y12 = atrack['y12']
1540 atrack_stereo12 = atrack['stereo12']
1541 atrack_34 = atrack['34']
1542 atrack_tagger = atrack['y_tagger']
1543 atrack_12 = list(atrack_y12) + list(atrack_stereo12)
1544 color = np.random.rand(3,)
1545 for ahit in atrack_12:
1546 xbot, xtop = ahit['xbot'], ahit['xtop']
1547 ybot, ytop = ahit['ybot'], ahit['ytop']
1548 plt.plot([ybot, ytop], [xbot, xtop], color=color)
1549
1550 plt.xlabel('y')
1551 plt.ylabel('x')
1552 plt.xlim(-50, 50)
1553 plt.ylim(-50, 50)
1554
1555 plt.savefig('pics/event_'+str(self.iEvent)+'.pdf', format='pdf', bbox_inches='tight')
1556 plt.clf()
1557 plt.close()
1558
1559
1560
1561
1562
1563 # Metrics
1564 if self.sTree.GetBranch("MufluxSpectrometerPoint"):
1565
1566
1567 (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.target_metrics(track_hits)
1568 h['Reco_target'].Fill("N total", n_tracks)
1569 h['Reco_target'].Fill("N recognized tracks", n_recognized)
1570 h['Reco_target'].Fill("N clones", n_clones)
1571 h['Reco_target'].Fill("N ghosts", n_ghosts)
1572 h['Reco_target'].Fill("N others", n_others)
1573
1574 (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_metrics(track_hits)
1575 h['Reco_muon'].Fill("N total", n_tracks)
1576 h['Reco_muon'].Fill("N recognized tracks", n_recognized)
1577 h['Reco_muon'].Fill("N clones", n_clones)
1578 h['Reco_muon'].Fill("N ghosts", n_ghosts)
1579 h['Reco_muon'].Fill("N others", n_others)
1580
1581 (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_with_tagger_metrics(track_hits)
1582 h['Reco_muon_with_tagger'].Fill("N total", n_tracks)
1583 h['Reco_muon_with_tagger'].Fill("N recognized tracks", n_recognized)
1584 h['Reco_muon_with_tagger'].Fill("N clones", n_clones)
1585 h['Reco_muon_with_tagger'].Fill("N ghosts", n_ghosts)
1586 h['Reco_muon_with_tagger'].Fill("N others", n_others)
1587
1588 (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_metrics_vs_p(track_hits)
1589
1590 (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.all_tracks_metrics(track_hits)
1591 h['Reco_all_tracks'].Fill("N total", n_tracks)
1592 h['Reco_all_tracks'].Fill("N recognized tracks", n_recognized)
1593 h['Reco_all_tracks'].Fill("N clones", n_clones)
1594 h['Reco_all_tracks'].Fill("N ghosts", n_ghosts)
1595 h['Reco_all_tracks'].Fill("N others", n_others)
1596
1597
1598 (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12,
1599 n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12,
1600 n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits)
1601
1602 h['NTrueTracks'].Fill("Stations 1&2, Y views", n_tracks_y12)
1603 h['NTrueTracks'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12)
1604 h['NTrueTracks'].Fill("Stations 3&4", n_tracks_34)
1605
1606 h['Reco_y12'].Fill("N total", n_tracks_y12)
1607 h['Reco_y12'].Fill("N recognized tracks", n_recognized_y12)
1608 h['Reco_y12'].Fill("N clones", n_clones_y12)
1609 h['Reco_y12'].Fill("N ghosts", n_ghosts_y12)
1610 h['Reco_y12'].Fill("N others", n_others_y12)
1611
1612 h['Reco_stereo12'].Fill("N total", n_tracks_stereo12)
1613 h['Reco_stereo12'].Fill("N recognized tracks", n_recognized_stereo12)
1614 h['Reco_stereo12'].Fill("N clones", n_clones_stereo12)
1615 h['Reco_stereo12'].Fill("N ghosts", n_ghosts_stereo12)
1616 h['Reco_stereo12'].Fill("N others", n_others_stereo12)
1617
1618 h['Reco_34'].Fill("N total", n_tracks_34)
1619 h['Reco_34'].Fill("N recognized tracks", n_recognized_34)
1620 h['Reco_34'].Fill("N clones", n_clones_34)
1621 h['Reco_34'].Fill("N ghosts", n_ghosts_34)
1622 h['Reco_34'].Fill("N others", n_others_34)
1623
1624
1625 (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12,
1626 n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12,
1627 n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits, mode='3hits')
1628
1629 h['NTrueTracks_3hits'].Fill("Stations 1&2, Y views", n_tracks_y12)
1630 h['NTrueTracks_3hits'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12)
1631 h['NTrueTracks_3hits'].Fill("Stations 3&4", n_tracks_34)
1632
1633 h['Reco_y12_3hits'].Fill("N total", n_tracks_y12)
1634 h['Reco_y12_3hits'].Fill("N recognized tracks", n_recognized_y12)
1635 h['Reco_y12_3hits'].Fill("N clones", n_clones_y12)
1636 h['Reco_y12_3hits'].Fill("N ghosts", n_ghosts_y12)
1637 h['Reco_y12_3hits'].Fill("N others", n_others_y12)
1638
1639 h['Reco_stereo12_3hits'].Fill("N total", n_tracks_stereo12)
1640 h['Reco_stereo12_3hits'].Fill("N recognized tracks", n_recognized_stereo12)
1641 h['Reco_stereo12_3hits'].Fill("N clones", n_clones_stereo12)
1642 h['Reco_stereo12_3hits'].Fill("N ghosts", n_ghosts_stereo12)
1643 h['Reco_stereo12_3hits'].Fill("N others", n_others_stereo12)
1644
1645 h['Reco_34_3hits'].Fill("N total", n_tracks_34)
1646 h['Reco_34_3hits'].Fill("N recognized tracks", n_recognized_34)
1647 h['Reco_34_3hits'].Fill("N clones", n_clones_34)
1648 h['Reco_34_3hits'].Fill("N ghosts", n_ghosts_34)
1649 h['Reco_34_3hits'].Fill("N others", n_others_34)
1650
1651
1652 (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12,
1653 n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12,
1654 n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits, mode='Tr4')
1655
1656 h['NTrueTracks_Tr4'].Fill("Stations 1&2, Y views", n_tracks_y12)
1657 h['NTrueTracks_Tr4'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12)
1658 h['NTrueTracks_Tr4'].Fill("Stations 3&4", n_tracks_34)
1659
1660 h['Reco_y12_Tr4'].Fill("N total", n_tracks_y12)
1661 h['Reco_y12_Tr4'].Fill("N recognized tracks", n_recognized_y12)
1662 h['Reco_y12_Tr4'].Fill("N clones", n_clones_y12)
1663 h['Reco_y12_Tr4'].Fill("N ghosts", n_ghosts_y12)
1664 h['Reco_y12_Tr4'].Fill("N others", n_others_y12)
1665
1666 h['Reco_stereo12_Tr4'].Fill("N total", n_tracks_stereo12)
1667 h['Reco_stereo12_Tr4'].Fill("N recognized tracks", n_recognized_stereo12)
1668 h['Reco_stereo12_Tr4'].Fill("N clones", n_clones_stereo12)
1669 h['Reco_stereo12_Tr4'].Fill("N ghosts", n_ghosts_stereo12)
1670 h['Reco_stereo12_Tr4'].Fill("N others", n_others_stereo12)
1671
1672 h['Reco_34_Tr4'].Fill("N total", n_tracks_34)
1673 h['Reco_34_Tr4'].Fill("N recognized tracks", n_recognized_34)
1674 h['Reco_34_Tr4'].Fill("N clones", n_clones_34)
1675 h['Reco_34_Tr4'].Fill("N ghosts", n_ghosts_34)
1676 h['Reco_34_Tr4'].Fill("N others", n_others_34)
1677
1678
1679
1680 n_true_track_hits_y12 = {}
1681 n_true_track_hits_stereo12 = {}
1682 n_true_track_hits_34 = {}
1683 for ahit in self.sTree.MufluxSpectrometerPoint:
1684
1685 track_id = ahit.GetTrackID()
1686
1687 detID = ahit.GetDetectorID()
1688 statnb = detID // 10000000
1689 vnb = (detID - statnb * 10000000) // 1000000
1690
1691 is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1)
1692 is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0)
1693 is_34 = (statnb == 3) + (statnb == 4)
1694
1695 if is_y12:
1696 if track_id in n_true_track_hits_y12:
1697 n_true_track_hits_y12[track_id] += 1
1698 else:
1699 n_true_track_hits_y12[track_id] = 1
1700 if is_stereo12:
1701 if track_id in n_true_track_hits_stereo12:
1702 n_true_track_hits_stereo12[track_id] += 1
1703 else:
1704 n_true_track_hits_stereo12[track_id] = 1
1705 if is_34:
1706 if track_id in n_true_track_hits_34:
1707 n_true_track_hits_34[track_id] += 1
1708 else:
1709 n_true_track_hits_34[track_id] = 1
1710
1711 for key in n_true_track_hits_y12.keys():
1712 n = n_true_track_hits_y12[key]
1713 h['NHits_true_y12'].Fill(n)
1714 for key in n_true_track_hits_stereo12.keys():
1715 n = n_true_track_hits_stereo12[key]
1716 h['NHits_true_stereo12'].Fill(n)
1717 for key in n_true_track_hits_34.keys():
1718 n = n_true_track_hits_34[key]
1719 h['NHits_true_34'].Fill(n)
1720
1721 for i_track in track_hits.keys():
1722 atrack = track_hits[i_track]
1723 atrack_y12 = atrack['y12']
1724 atrack_stereo12 = atrack['stereo12']
1725 atrack_34 = atrack['34']
1726
1727 if len(atrack_y12) > 0:
1728 h['NHits_reco_y12'].Fill(len(atrack_y12))
1729 if len(atrack_stereo12) > 0:
1730 h['NHits_reco_stereo12'].Fill(len(atrack_stereo12))
1731 if len(atrack_34) > 0:
1732 h['NHits_reco_34'].Fill(len(atrack_34))
1733
1734
1735
1736
1737 for entry in trackCandidates:
1738 #check
1739 #print "fitting with stereo"
1740 atrack = entry[1]
1741 theTrack = entry[0]
1742 if not theTrack.checkConsistency():
1743 print('Problem with track before fit, not consistent',atrack,theTrack)
1744 continue
1745 # do the fit
1746 try: self.fitter.processTrack(theTrack) # processTrackWithRep(theTrack,rep,True)
1747 except:
1748 print("genfit failed to fit track")
1749 continue
1750 #check
1751 if not theTrack.checkConsistency():
1752 #print 'Problem with track after fit, not consistent',atrack,theTrack
1753 continue
1754 fitStatus = theTrack.getFitStatus()
1755 nmeas = fitStatus.getNdf()
1756 chi2 = fitStatus.getChi2()/nmeas
1757 pvalue = fitStatus.getPVal()
1758 #if pvalue < 0.05:
1759 # print "P value too low. Rejecting track."
1760 # continue
1761 h['nmeas'].Fill(nmeas)
1762 h['chi2'].Fill(chi2)
1763 h['p-value'].Fill(pvalue)
1764 try:
1765 fittedState = theTrack.getFittedState()
1766 fittedMom = fittedState.getMomMag()
1767 h['p-fittedtracks'].Fill(fittedMom)
1768 h['1/p-fittedtracks'].Fill(1./fittedMom)
1769 Px,Py,Pz = fittedState.getMom().x(),fittedState.getMom().y(),fittedState.getMom().z()
1770 P = fittedMom
1771 Pt = ROOT.TMath.Sqrt(Px*Px+Py*Py)
1772 h['p/pt'].Fill(P,Pt)
1773 h['pt-fittedtracks'].Fill(Pt)
1774 h['1/pt-fittedtracks'].Fill(1./Pt)
1775
1776 if self.sTree.GetBranch("MufluxSpectrometerPoint"):
1777
1778 atrack_ids = []
1779 for digi_hit in trackDigiHits[atrack]:
1780 ahit = self.sTree.MufluxSpectrometerPoint[digi_hit]
1781 atrack_ids.append(ahit.GetTrackID())
1782 frac, tmax = self.fracMCsame(atrack_ids)
1783 Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(tmax) # MC Truth
1784 Pgun,Pgunx,Pguny,Pgunz = self.getPtruthAtOrigin(tmax) # MC Truth
1785 Pttruth = ROOT.TMath.Sqrt(Ptruthx*Ptruthx+Ptruthy*Ptruthy)
1786 h['p/pt_truth'].Fill(Ptruth,Pttruth)
1787 perr = (P - Ptruth) / Ptruth
1788 pterr = (Pt - Pttruth) / Pttruth
1789 h['p_rel_error'].Fill(perr)
1790 h['pt_rel_error'].Fill(pterr)
1791
1792 mom_init = trackMomentums[atrack]
1793 print("P_precalc, P_truth: ", mom_init, Ptruth)
1794
1795 if Pz !=0:
1796 pxpzfitted = Px/Pz
1797 pypzfitted = Py/Pz
1798 if Ptruthz !=0:
1799 pxpztrue = Ptruthx/Ptruthz
1800 pypztrue = Ptruthy/Ptruthz
1801 h['Px/Pzfitted'].Fill(pxpzfitted)
1802 h['Py/Pzfitted'].Fill(pypzfitted)
1803 h['Px/Pztrue'].Fill(pxpztrue)
1804 h['Py/Pztrue'].Fill(pypztrue)
1805 h['Px/Pzfitted-Px/Pztruth'].Fill(Ptruth,pxpzfitted-pxpztrue)
1806 h['Py/Pzfitted-Py/Pztruth'].Fill(Ptruth,pypzfitted-pypztrue)
1807 h['ptruth'].Fill(Ptruth)
1808 delPOverP = (P/Ptruth)-1
1809 invdelPOverP = (Ptruth/P)-1
1810 if 1==0:
1811 if invdelPOverP < -0.8:
1812 print("invdelPOverP = ",invdelPOverP)
1813 print("Ptruth =",Ptruth," Pfitted =",P)
1814 for n in range(hitPosLists[atrack].size()):
1815 print("hit=",n," x(top) ",hitPosLists[atrack][n][0]," y(top) ",hitPosLists[atrack][n][1]," z ",hitPosLists[atrack][n][2]," x(bot) ",hitPosLists[atrack][n][3]," y(bot) ", hitPosLists[atrack][n][4], " dist ", hitPosLists[atrack][n][6])
1816 nMufluxHits = self.sTree.MufluxSpectrometerPoint.GetEntriesFast()
1817 for i in range(nMufluxHits):
1818 MufluxHit = self.sTree.MufluxSpectrometerPoint[i]
1819 if ((hitPosLists[atrack][n][0]+1.8 > MufluxHit.GetX()) or(hitPosLists[atrack][n][3]+1.8 > MufluxHit.GetX())) and ((hitPosLists[atrack][n][0]-1.8<MufluxHit.GetX()) or (hitPosLists[atrack][n][3]-1.8<MufluxHit.GetX())) and (hitPosLists[atrack][n][2]+1.>MufluxHit.GetZ()) and (hitPosLists[atrack][n][2]-1.<MufluxHit.GetZ()):
1820 print("hit x=",MufluxHit.GetX()," y=",MufluxHit.GetY()," z=",MufluxHit.GetZ())
1821
1822
1823 h['delPOverP'].Fill(Ptruth,delPOverP)
1824 h['invdelPOverP'].Fill(Ptruth,invdelPOverP)
1825 h['deltaPOverP'].Fill(Ptruth,delPOverP)
1826 h['Pfitted-Pgun'].Fill(Pgun,P)
1827 #print "end fitting with stereo"
1828
1829 except:
1830 print("problem with fittedstate")
1831 continue
1832
1833 #if 1==0:
1834 for entry in trackCandidates_noT4:
1835 #check
1836 #print "fitting without stereo hits"
1837 atrack = entry[1]
1838 theTrack = entry[0]
1839 if not theTrack.checkConsistency():
1840 print('Problem with track before fit, not consistent',atrack,theTrack)
1841 continue
1842 # do the fit
1843 try: self.fitter.processTrack(theTrack) # processTrackWithRep(theTrack,rep,True)
1844 except:
1845 print("genfit failed to fit track")
1846 continue
1847 #check
1848 if not theTrack.checkConsistency():
1849 print('Problem with track after fit, not consistent',atrack,theTrack)
1850 continue
1851 fitStatus = theTrack.getFitStatus()
1852 nmeas = fitStatus.getNdf()
1853 chi2 = fitStatus.getChi2()/nmeas
1854 pvalue = fitStatus.getPVal()
1855 #if pvalue < 0.05:
1856 # print "P value too low. Rejecting track."
1857 # continue
1858 h['nmeas-noT4'].Fill(nmeas)
1859 h['chi2-noT4'].Fill(chi2)
1860 h['p-value-noT4'].Fill(pvalue)
1861 try:
1862
1863 fittedState = theTrack.getFittedState()
1864 fittedMom = fittedState.getMomMag()
1865 h['p-fittedtracks-noT4'].Fill(fittedMom)
1866 h['1/p-fittedtracks-noT4'].Fill(1./fittedMom)
1867 Px,Py,Pz = fittedState.getMom().x(),fittedState.getMom().y(),fittedState.getMom().z()
1868 P = fittedMom
1869 Pt = ROOT.TMath.Sqrt(Px*Px+Py*Py)
1870 h['p/pt_noT4'].Fill(P,Pt)
1871 h['pt-fittedtracks-noT4'].Fill(Pt)
1872 h['1/pt-fittedtracks-noT4'].Fill(1./Pt)
1873
1874 if self.sTree.GetBranch("MufluxSpectrometerPoint"):
1875
1876 atrack_ids = []
1877 for digi_hit in trackDigiHits_noT4[atrack]:
1878 ahit = self.sTree.MufluxSpectrometerPoint[digi_hit]
1879 atrack_ids.append(ahit.GetTrackID())
1880 frac, tmax = self.fracMCsame(atrack_ids)
1881 Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(tmax)
1882 Pgun,Pgunx,Pguny,Pgunz = self.getPtruthAtOrigin(tmax)
1883 Pttruth = ROOT.TMath.Sqrt(Ptruthx*Ptruthx+Ptruthy*Ptruthy)
1884 h['p/pt_truth_noT4'].Fill(Ptruth,Pttruth)
1885 perr = (P - Ptruth) / Ptruth
1886 pterr = (Pt - Pttruth) / Pttruth
1887 h['p_rel_error_noT4'].Fill(perr)
1888 h['pt_rel_error_noT4'].Fill(pterr)
1889
1890 if Pz !=0:
1891 pxpzfitted = Px/Pz
1892 pypzfitted = Py/Pz
1893 if Ptruthz !=0:
1894 pxpztrue = Ptruthx/Ptruthz
1895 pypztrue = Ptruthy/Ptruthz
1896 h['Px/Pzfitted-Px/Pztruth-noT4'].Fill(Ptruth,pxpzfitted-pxpztrue)
1897 h['Py/Pzfitted-Py/Pztruth-noT4'].Fill(Ptruth,pypzfitted-pypztrue)
1898 h['Px/Pzfitted-noT4'].Fill(pxpzfitted)
1899 h['Py/Pzfitted-noT4'].Fill(pypzfitted)
1900 h['Px/Pztrue-noT4'].Fill(pxpztrue)
1901 h['Py/Pztrue-noT4'].Fill(pypztrue)
1902
1903 h['ptruth-noT4'].Fill(Ptruth)
1904 delPOverP = (P/Ptruth)-1
1905 invdelPOverP = (Ptruth/P)-1
1906 h['delPOverP-noT4'].Fill(Ptruth,delPOverP)
1907 h['invdelPOverP-noT4'].Fill(Ptruth,invdelPOverP)
1908 h['deltaPOverP-noT4'].Fill(Ptruth,delPOverP)
1909 h['Pfitted-Pgun-noT4'].Fill(Pgun,P)
1910 #print "end fitting without stereo hits"
1911 except:
1912 print("noT4 track: problem with fittedstate")
1913 continue
1914
1915
1916 for entry in trackCandidates_all:
1917 #check
1918 #print "fitting without stereo hits"
1919 atrack = entry[1]
1920 theTrack = entry[0]
1921 if not theTrack.checkConsistency():
1922 print('Problem with track before fit, not consistent',atrack,theTrack)
1923 continue
1924 # do the fit
1925 try: self.fitter.processTrack(theTrack) # processTrackWithRep(theTrack,rep,True)
1926 except:
1927 print("genfit failed to fit track")
1928 continue
1929 #check
1930 if not theTrack.checkConsistency():
1931 print('Problem with track after fit, not consistent',atrack,theTrack)
1932 continue
1933 fitStatus = theTrack.getFitStatus()
1934 nmeas = fitStatus.getNdf()
1935 chi2 = fitStatus.getChi2()/nmeas
1936 pvalue = fitStatus.getPVal()
1937 #if pvalue < 0.05:
1938 # print "P value too low. Rejecting track."
1939 # continue
1940 h['nmeas-all'].Fill(nmeas)
1941 h['chi2-all'].Fill(chi2)
1942 h['p-value-all'].Fill(pvalue)
1943 try:
1944
1945 fittedState = theTrack.getFittedState()
1946 fittedMom = fittedState.getMomMag()
1947 h['p-fittedtracks-all'].Fill(fittedMom)
1948 h['1/p-fittedtracks-all'].Fill(1./fittedMom)
1949 Px,Py,Pz = fittedState.getMom().x(),fittedState.getMom().y(),fittedState.getMom().z()
1950 P = fittedMom
1951 Pt = ROOT.TMath.Sqrt(Px*Px+Py*Py)
1952 h['p/pt_all'].Fill(P,Pt)
1953 h['pt-fittedtracks-all'].Fill(Pt)
1954 h['1/pt-fittedtracks-all'].Fill(1./Pt)
1955
1956 if self.sTree.GetBranch("MufluxSpectrometerPoint"):
1957
1958 atrack_ids = []
1959 for digi_hit in trackDigiHits_noT4[atrack]:
1960 ahit = self.sTree.MufluxSpectrometerPoint[digi_hit]
1961 atrack_ids.append(ahit.GetTrackID())
1962 frac, tmax = self.fracMCsame(atrack_ids)
1963 Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(tmax)
1964 Pgun,Pgunx,Pguny,Pgunz = self.getPtruthAtOrigin(tmax)
1965 Pttruth = ROOT.TMath.Sqrt(Ptruthx*Ptruthx+Ptruthy*Ptruthy)
1966 h['p/pt_truth_all'].Fill(Ptruth,Pttruth)
1967 perr = (P - Ptruth) / Ptruth
1968 pterr = (Pt - Pttruth) / Pttruth
1969 h['p_rel_error_all'].Fill(perr)
1970 h['pt_rel_error_all'].Fill(pterr)
1971
1972 if Pz !=0:
1973 pxpzfitted = Px/Pz
1974 pypzfitted = Py/Pz
1975 if Ptruthz !=0:
1976 pxpztrue = Ptruthx/Ptruthz
1977 pypztrue = Ptruthy/Ptruthz
1978 h['Px/Pzfitted-Px/Pztruth-all'].Fill(Ptruth,pxpzfitted-pxpztrue)
1979 h['Py/Pzfitted-Py/Pztruth-all'].Fill(Ptruth,pypzfitted-pypztrue)
1980 h['Px/Pzfitted-all'].Fill(pxpzfitted)
1981 h['Py/Pzfitted-all'].Fill(pypzfitted)
1982 h['Px/Pztrue-all'].Fill(pxpztrue)
1983 h['Py/Pztrue-all'].Fill(pypztrue)
1984
1985 h['ptruth-all'].Fill(Ptruth)
1986 delPOverP = (P/Ptruth)-1
1987 invdelPOverP = (Ptruth/P)-1
1988 h['delPOverP-all'].Fill(Ptruth,delPOverP)
1989 h['invdelPOverP-all'].Fill(Ptruth,invdelPOverP)
1990 h['deltaPOverP-all'].Fill(Ptruth,delPOverP)
1991 h['Pfitted-Pgun-all'].Fill(Pgun,P)
1992 #print "end fitting without stereo hits"
1993 except:
1994 print("All tracks: problem with fittedstate")
1995 continue
1996
1997 # make track persistent
1998 nTrack = self.fGenFitArray.GetEntries()
1999 if not global_variables.debug:
2000 theTrack.prune("CFL") # http://sourceforge.net/p/genfit/code/HEAD/tree/trunk/core/include/Track.h#l280
2001 self.fGenFitArray[nTrack] = theTrack
2002 self.fitTrack2MC.push_back(atrack)
2003 if global_variables.debug:
2004 print('save track',theTrack,chi2,nM,fitStatus.isFitConverged())
2005
2006 self.fitTracks.Fill()
2007 self.mcLink.Fill()
2008 return nTrack+1
2009
2010 def finish(self):
2011 del self.fitter
2012 print('finished writing tree')
2013 self.sTree.Write()
2014 ut.errorSummary()
2015 ut.writeHists(h,"recohists.root")
2016 # if global_variables.realPR:
2017 # ut.writeHists(shipPatRec.h,"recohists_patrec.root")
set
set(INCLUDE_DIRECTORIES ${SYSTEM_INCLUDE_DIRECTORIES} ${VMC_INCLUDE_DIRS} ${CMAKE_SOURCE_DIR}/shipdata ${CMAKE_SOURCE_DIR}/shipLHC ${CMAKE_SOURCE_DIR}/analysis/cuts ${CMAKE_SOURCE_DIR}/analysis/tools ${FMT_INCLUDE_DIR}) include_directories($
Definition CMakeLists.txt:1
MufluxDigiReco.MufluxDigiReco.RT
RT(self, s, t, function='parabola')
Definition MufluxDigiReco.py:431
MufluxDigiReco.MufluxDigiReco.muon_metrics_vs_p
muon_metrics_vs_p(self, track_hit_ids)
Definition MufluxDigiReco.py:1094
MufluxDigiReco.MufluxDigiReco.target_metrics
target_metrics(self, track_hit_ids)
Definition MufluxDigiReco.py:757
MufluxDigiReco.MufluxDigiReco.muon_metrics
muon_metrics(self, track_hit_ids)
Definition MufluxDigiReco.py:939
MufluxDigiReco.MufluxDigiReco.recognition_metrics
recognition_metrics(self, track_hit_ids, mode='all')
Definition MufluxDigiReco.py:553
MufluxDigiReco.MufluxDigiReco.digitizeMuonTagger
digitizeMuonTagger(self, fake_clustering=True)
Definition MufluxDigiReco.py:162
MufluxDigiReco.MufluxDigiReco.muon_with_tagger_metrics
muon_with_tagger_metrics(self, track_hit_ids)
Definition MufluxDigiReco.py:1012
MufluxDigiReco.MufluxDigiReco.all_tracks_metrics
all_tracks_metrics(self, track_hit_ids)
Definition MufluxDigiReco.py:868
MufluxPatRec.execute
execute(SmearedHits, TaggerHits, withNTaggerHits, withDist2Wire, debug=0)
Definition MufluxPatRec.py:15