SND@LHC Software
Loading...
Searching...
No Matches
mufiHitMaps.py
Go to the documentation of this file.
1#python -i mufiMIP.py -r 46 -p /eos/experiment/sndlhc/testbeam/MuFilter/TB_data_commissioning/sndsw/ -g geofile_sndlhc_H6.root
2
3import ROOT,os
4import rootUtils as ut
5import shipunit as u
6h={}
7from argparse import ArgumentParser
8parser = ArgumentParser()
9parser.add_argument("-r", "--runNumber", dest="runNumber", help="run number", type=int,required=True)
10parser.add_argument("-p", "--path", dest="path", help="run number",required=False,default="")
11parser.add_argument("-f", "--inputFile", dest="fname", help="file name for MC", type=str,default=None,required=False)
12parser.add_argument("-g", "--geoFile", dest="geoFile", help="geofile", required=True)
13options = parser.parse_args()
14
15path = options.path
16if path.find('eos')>0: path = os.environ['EOSSHIP']+options.path
17
18import SndlhcGeo
19geo = SndlhcGeo.GeoInterface(path+options.geoFile)
20MuFilter = geo.modules['MuFilter']
21
22A,B = ROOT.TVector3(),ROOT.TVector3()
23zPos={}
24systemAndPlanes = {1:2,2:5,3:7}
25systemAndBars = {1:7,2:10,3:60}
26
27freq = 160.316E6
28
29for s in systemAndPlanes:
30 for plane in range(systemAndPlanes[s]):
31 bar = 4
32 p = plane
33 if s==3 and plane%2==0:
34 bar = 90
35 p = plane//2
36 if s==3 and plane%2==1:
37 bar = 30
38 p = plane//2
39 MuFilter.GetPosition(s*10000+p*1000+bar,A,B)
40 zPos[s*10+plane] = (A.Z()+B.Z())/2.
41
42def smallSiPMchannel(i):
43 if i==2 or i==5 or i==10 or i==13: return True
44 else: return False
45
46if options.runNumber>0:
47 f=ROOT.TFile.Open(path+'sndsw_raw_'+str(options.runNumber).zfill(6)+'.root')
48 eventTree = f.rawConv
49else:
50 f=ROOT.TFile.Open(options.fname)
51 eventTree = f.cbmsim
52
53# backward compatbility for early converted events
54eventTree.GetEvent(0)
55if eventTree.GetBranch('Digi_MuFilterHit'): eventTree.Digi_MuFilterHits = eventTree.Digi_MuFilterHit
56
57import SndlhcTracking
58trackTask = SndlhcTracking.Tracking()
59trackTask.InitTask(eventTree)
60
61def hitMaps(Nev = -1):
62 # veto system 2 layers with 7 bars and 8 sipm channels on both ends
63 # US system 5 layers with 10 bars and 8 sipm channels on both ends
64 # DS system horizontal(3) planes, 60 bars, readout on both sides, single channel
65 # vertical(4) planes, 60 bar, readout on top, single channel
66 for s in systemAndPlanes:
67 for l in range(systemAndPlanes[s]):
68 ut.bookHist(h,'hit_'+str(s*10+l),'channel map / plane '+str(s*10+l),160,-0.5,159.5)
69 if s==3: ut.bookHist(h,'bar_'+str(s*10+l),'hit map / plane '+str(s*10+l),60,-0.5,59.5)
70 else: ut.bookHist(h,'bar_'+str(s*10+l),'hit map / plane '+str(s*10+l),10,-0.5,9.5)
71 ut.bookHist(h,'sig_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
72 if s==2: ut.bookHist(h,'sigS_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
73 ut.bookHist(h,'sigL_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
74 ut.bookHist(h,'sigR_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
75 ut.bookHist(h,'Tsig_'+str(s*10+l),'signal / plane '+str(s*10+l),200,0.0,200.)
76 ut.bookHist(h,'occ_'+str(s*10+l),'channel occupancy '+str(s*10+l),100,0.0,200.)
77 ut.bookHist(h,'occTag_'+str(s*10+l),'channel occupancy '+str(s*10+l),100,0.0,200.)
78
79 ut.bookHist(h,'leftvsright_1','Veto hits in left / right',10,-0.5,9.5,10,-0.5,9.5)
80 ut.bookHist(h,'leftvsright_2','US hits in left / right',10,-0.5,9.5,10,-0.5,9.5)
81 ut.bookHist(h,'leftvsright_3','DS hits in left / right',2,-0.5,1.5,2,-0.5,1.5)
82 ut.bookHist(h,'leftvsright_signal_1','Veto signal in left / right',100,-0.5,200.,100,-0.5,200.)
83 ut.bookHist(h,'leftvsright_signal_2','US signal in left / right',100,-0.5,200.,100,-0.5,200.)
84 ut.bookHist(h,'leftvsright_signal_3','DS signal in left / right',100,-0.5,200.,100,-0.5,200.)
85
86 ut.bookHist(h,'dtime','delta event time; dt [ns]',100,0.0,1000.)
87 ut.bookHist(h,'dtimeu','delta event time; dt [us]',100,0.0,1000.)
88 ut.bookHist(h,'dtimem','delta event time; dt [ms]',100,0.0,1000.)
89
90 ut.bookHist(h,'bs','beam spot',100,-100.,10.,100,0.,80.)
91 ut.bookHist(h,'bsDS','beam spot',60,-0.5,59.5,60,-0.5,59.5)
92 ut.bookHist(h,'slopes','track slopes',100,-0.1,0.1,100,-0.1,0.1)
93
94 N=-1
95 Tprev = 0
96 if Nev < 0 : Nev = eventTree.GetEntries()
97 eventTree.GetEvent(0)
98 t0 = eventTree.EventHeader.GetEventTime()/freq
99
100 for event in eventTree:
101 N+=1
102 if N>Nev: break
103 withX = False
104 planes = {}
105 for aHit in event.Digi_MuFilterHits:
106 if not aHit.isValid(): continue
107 detID = aHit.GetDetectorID()
108 if aHit.isVertical(): withX = True
109 s = detID//10000
110 l = (detID%10000)//1000 # plane number
111 bar = (detID%1000)
112 key = s*100+l
113 if s>2:
114 l=2*l
115 if bar>59:
116 bar=bar-60
117 l+=1
118 if not key in planes: planes[key] = {}
119 sumSignal = aHit.SumOfSignals()
120 planes[key][bar] = [sumSignal['SumL'],sumSignal['SumR']]
121 nSiPMs = aHit.GetnSiPMs()
122 nSides = aHit.GetnSides()
123
124# check left/right
125 allChannels = aHit.GetAllSignals(False) # masking not yet correctly done in the raw conversion
126 if nSides==2:
127 Nleft = 0
128 Nright = 0
129 Sleft = 0
130 Sright = 0
131 for c in allChannels:
132 if nSiPMs > c[0]: # left side
133 Nleft+=1
134 Sleft+=c[1]
135 else:
136 Nright+=1
137 Sright+=c[1]
138 rc = h['leftvsright_'+str(s)].Fill(Nleft,Nright)
139 rc = h['leftvsright_signal_'+str(s)].Fill(Sleft,Sright)
140 for c in allChannels:
141 channel = bar*nSiPMs*nSides + c[0]
142 rc = h['hit_'+str(s)+str(l)].Fill( int(channel))
143 rc = h['bar_'+str(s)+str(l)].Fill(bar)
144 if s==2 and smallSiPMchannel(c[0]) : rc = h['sigS_'+str(s)+str(l)].Fill(c[1])
145 elif c[0]<nSiPMs: rc = h['sigL_'+str(s)+str(l)].Fill(c[1])
146 else : rc = h['sigR_'+str(s)+str(l)].Fill(c[1])
147 rc = h['sig_'+str(s)+str(l)].Fill(c[1])
148 maxOneBar = True
149 for key in planes:
150 if len(planes[key]) > 2: maxOneBar = False
151 if withX and maxOneBar: beamSpot()
152
153 installed_stations = {}
154 for s in range(1,4):
155 for l in range(systemAndPlanes[s]):
156 if h['hit_'+str(s)+str(l)].GetEntries()>0:
157 if not s in installed_stations: installed_stations[s]=0
158 installed_stations[s]+=1
159 x = 0
160 y = 0
161 for s in installed_stations:
162 if installed_stations[s]>0: x+=1
163 if installed_stations[s]>y: y = installed_stations[s]
164 ut.bookCanvas(h,'hitmaps',' ',1200,1600,x,y)
165 ut.bookCanvas(h,'barmaps',' ',1200,1600,x,y)
166 ut.bookCanvas(h,'signal',' ',1200,1600,x,y)
167 ut.bookCanvas(h,'Tsignal',' ',1200,1600,x,y)
168
169 for S in installed_stations:
170 for l in range(installed_stations[S]):
171 n = S-1 + l*x
172 tc = h['hitmaps'].cd(n)
173 h['hit_'+str(S)+str(l)].Draw()
174 tc = h['barmaps'].cd(n)
175 h['bar_'+str(S)+str(l)].Draw()
176 tc = h['signal'].cd(n)
177 h['sig_'+str(S)+str(l)].Draw()
178 tc = h['Tsignal'].cd(n)
179 h['Tsig_'+str(S)+str(l)].Draw()
180
181 ut.bookCanvas(h,'USBars',' ',1200,900,1,1)
182 colours = {0:ROOT.kOrange,1:ROOT.kRed,2:ROOT.kGreen,3:ROOT.kBlue,4:ROOT.kMagenta}
183 for i in range(5):
184 h['bar_2'+str(i)].SetLineColor(colours[i])
185 h['bar_2'+str(i)].SetLineWidth(2)
186 h['bar_2'+str(i)].SetStats(0)
187 h['bar_20'].Draw()
188 h['bar_21'].Draw('same')
189 h['bar_22'].Draw('same')
190 h['bar_23'].Draw('same')
191 h['bar_24'].Draw('same')
192 h['lbar2']=ROOT.TLegend(0.6,0.6,0.99,0.99)
193 for i in range(5):
194 h['lbar2'].AddEntry(h['bar_2'+str(i)],'plane '+str(i+1),"f")
195 h['lbar2'].Draw()
196 for i in range(4):
197 h['hit_3'+str(i)].SetLineColor(colours[i])
198 h['hit_3'+str(i)].SetLineWidth(2)
199 h['hit_3'+str(i)].SetStats(0)
200 h['hit_30'].Draw()
201 h['hit_31'].Draw('same')
202 h['hit_32'].Draw('same')
203 h['hit_33'].Draw('same')
204 h['lbar3']=ROOT.TLegend(0.6,0.6,0.99,0.99)
205 for i in range(4):
206 h['lbar3'].AddEntry(h['hit_3'+str(i)],'plane '+str(i+1),"f")
207 h['lbar3'].Draw()
208
209 ut.bookCanvas(h,'LR',' ',1800,900,3,2)
210 h['LR'].cd(1)
211 h['leftvsright_'+str(1)].Draw('textBox')
212 h['LR'].cd(2)
213 h['leftvsright_'+str(2)].Draw('textBox')
214 h['LR'].cd(3)
215 h['leftvsright_'+str(3)].Draw('textBox')
216 h['LR'].cd(4)
217 h['leftvsright_signal_1'].SetMaximum(h['leftvsright_signal_1'].GetBinContent(10,10))
218 h['leftvsright_signal_2'].SetMaximum(h['leftvsright_signal_2'].GetBinContent(10,10))
219 h['leftvsright_signal_3'].SetMaximum(h['leftvsright_signal_3'].GetBinContent(10,10))
220 h['leftvsright_signal_'+str(1)].Draw('colz')
221 h['LR'].cd(5)
222 h['leftvsright_signal_'+str(2)].Draw('colz')
223 h['LR'].cd(6)
224 h['leftvsright_signal_'+str(3)].Draw('colz')
225
226 ut.bookCanvas(h,'LRinEff',' ',1800,450,3,1)
227 for s in range(1,4):
228 h['lLRinEff'+str(s)]=ROOT.TLegend(0.6,0.54,0.99,0.93)
229 name = 'leftvsright_signal_'+str(s)
230 h[name+'0Y'] = h[name].ProjectionY(name+'0Y',1,1)
231 h[name+'0X'] = h[name].ProjectionX(name+'0X',1,1)
232 h[name+'1X'] = h[name].ProjectionY(name+'1Y')
233 h[name+'1Y'] = h[name].ProjectionX(name+'1X')
234 tc = h['LRinEff'].cd(s)
235 tc.SetLogy()
236 h[name+'0X'].SetStats(0)
237 h[name+'0Y'].SetStats(0)
238 h[name+'1X'].SetStats(0)
239 h[name+'1Y'].SetStats(0)
240 h[name+'0X'].SetLineColor(ROOT.kRed)
241 h[name+'0Y'].SetLineColor(ROOT.kGreen)
242 h[name+'1X'].SetLineColor(ROOT.kMagenta)
243 h[name+'1Y'].SetLineColor(ROOT.kCyan)
244 h[name+'0X'].SetMaximum(max(h[name+'1X'].GetMaximum(),h[name+'1Y'].GetMaximum()))
245 h[name+'0X'].Draw()
246 h[name+'0Y'].Draw('same')
247 h[name+'1X'].Draw('same')
248 h[name+'1Y'].Draw('same')
249 # Fill(Sleft,Sright)
250 h['lLRinEff'+str(s)].AddEntry(h[name+'0X'],'left with no signal right',"f")
251 h['lLRinEff'+str(s)].AddEntry(h[name+'0Y'],'right with no signal left',"f")
252 h['lLRinEff'+str(s)].AddEntry(h[name+'1X'],'left all',"f")
253 h['lLRinEff'+str(s)].AddEntry(h[name+'1Y'],'right all',"f")
254 h['lLRinEff'+str(s)].Draw()
255
256
257
258 ut.bookCanvas(h,'signalUSVeto',' ',1200,1600,3,7)
259 s = 1
260 l = 1
261 for plane in range(2):
262 for side in ['L','R','S']:
263 tc = h['signalUSVeto'].cd(l)
264 l+=1
265 if side=='S': continue
266 h['sig'+side+'_'+str( s*10+plane)].Draw()
267 s=2
268 for plane in range(5):
269 for side in ['L','R','S']:
270 tc = h['signalUSVeto'].cd(l)
271 l+=1
272 h['sig'+side+'_'+str( s*10+plane)].Draw()
273 ut.bookCanvas(h,'signalDS',' ',900,1600,2,7)
274 s = 3
275 l = 1
276 for plane in range(7):
277 for side in ['L','R']:
278 tc = h['signalDS'].cd(l)
279 l+=1
280 h['sig'+side+'_'+str( s*10+plane)].Draw()
281
282
283 for canvas in ['signalUSVeto','LR','hitmaps','barmaps','signal','Tsignal','USBars']:
284 h[canvas].Update()
285 h[canvas].Print(canvas+'-run'+str(options.runNumber)+'.png')
286
287def eventTime(Nev=-1):
288 Tprev = -1
289 if Nev < 0 : Nev = eventTree.GetEntries()
290 ut.bookHist(h,'Etime','delta event time; dt [s]',100,0.0,1.)
291 ut.bookHist(h,'EtimeZ','delta event time; dt [ns]',1000,0.0,10000.)
292 ut.bookCanvas(h,'T',' ',1024,2*768,1,2)
293 eventTree.GetEvent(0)
294 t0 = eventTree.EventHeader.GetEventTime()/160.E6
295 eventTree.GetEvent(Nev-1)
296 tmax = eventTree.EventHeader.GetEventTime()/160.E6
297 ut.bookHist(h,'time','elapsed time; t [s]',1000,0,tmax-t0)
298
299 N=-1
300 for event in eventTree:
301 N+=1
302 if N>Nev: break
303 T = event.EventHeader.GetEventTime()
304 dT = 0
305 if Tprev >0: dT = T-Tprev
306 Tprev = T
307 rc = h['Etime'].Fill(dT/freq)
308 rc = h['EtimeZ'].Fill(dT*1E9/160.316E6)
309 rc = h['time'].Fill( (T/freq-t0))
310 tc = h['T'].cd(1)
311 tc.SetLogy(1)
312 h['EtimeZ'].Draw()
313 tc.Update()
314 tc = h['T'].cd(2)
315 h['time'].Draw()
316 h['T'].Update()
317 h['T'].Print('time-run'+str(options.runNumber)+'.png')
318
319def beamSpot():
320 trackTask.ExecuteTask()
321 Xbar = -10
322 Ybar = -10
323 for aTrack in eventTree.fittedTracks:
324 state = aTrack.getFittedState()
325 pos = state.getPos()
326 rc = h['bs'].Fill(pos.x(),pos.y())
327 points = aTrack.getPoints()
328 for p in points:
329 m = p.getRawMeasurement()
330 detID = m.getDetId()
331 key = m.getHitId()//1000 # for mufi
332 aHit = eventTree.Digi_MuFilterHits[key]
333 if aHit.GetDetectorID() != detID: continue # not a Mufi hit
334 s = detID//10000
335 l = (detID%10000)//1000 # plane number
336 bar = (detID%1000)
337 if s>2:
338 l=2*l
339 if bar>59:
340 bar=bar-60
341 l+=1
342 if s==3 and l%2==0: Ybar=bar
343 if s==3 and l%2==1: Xbar=bar
344 nSiPMs = aHit.GetnSiPMs()
345 nSides = aHit.GetnSides()
346 for p in range(nSides):
347 c=bar*nSiPMs*nSides + p*nSiPMs
348 for i in range(nSiPMs):
349 signal = aHit.GetSignal(i+p*nSiPMs)
350 if signal > 0:
351 rc = h['Tsig_'+str(s)+str(l)].Fill(signal)
352 mom = state.getMom()
353 slopeY= mom.X()/mom.Z()
354 slopeX= mom.Y()/mom.Z()
355 h['slopes'].Fill(slopeX,slopeY)
356 if not Ybar<0 and not Xbar<0 and abs(slopeY)<0.01: rc = h['bsDS'].Fill(Xbar,Ybar)
357
358 aTrack.Delete()
359
360def DS_track():
361# check for low occupancy and enough hits in DS stations
362 stations = {}
363 for s in systemAndPlanes:
364 for plane in range(systemAndPlanes[s]):
365 stations[s*10+plane] = {}
366 k=-1
367 for aHit in eventTree.Digi_MuFilterHits:
368 k+=1
369 if not aHit.isValid(): continue
370 s = aHit.GetDetectorID()//10000
371 p = (aHit.GetDetectorID()//1000)%10
372 bar = aHit.GetDetectorID()%1000
373 plane = s*10+p
374 if s==3:
375 if bar<60: plane = s*10+2*p
376 else: plane = s*10+2*p+1
377 stations[plane][k] = aHit
378 if not len(stations[30])*len(stations[31])*len(stations[32])*len(stations[33]) == 1: return -1
379# build trackCandidate
380 hitlist = {}
381 for p in range(30,34):
382 k = list(stations[p].keys())[0]
383 hitlist[k] = stations[p][k]
384 theTrack = trackTask.fitTrack(hitlist)
385 return theTrack
386
387def USEfficiency(Nev=-1):
388 name = {1:'Veto',2:'US',3:'DS'}
389 for s in systemAndPlanes:
390 for l in range(systemAndPlanes[s]):
391 ut.bookHist(h,'dtLRvsX_'+name[s]+str(s*10+l),'dt vs x track '+str(s*10+l)+";X [cm]; dt [ns]",80,-70.,10.,100,-10.,10.)
392 if s!=2 : continue
393 for l in range(systemAndPlanes[s]):
394 ut.bookHist(h,'resY_'+name[s]+str(s*10+l),'residual Y '+str(s*10+l),100,-20.,20.)
395 ut.bookHist(h,'resY_'+name[s]+'L'+str(s*10+l),'residual Y '+str(s*10+l),100,-20.,20.)
396 ut.bookHist(h,'resY_'+name[s]+'R'+str(s*10+l),'residual Y '+str(s*10+l),100,-20.,20.)
397 ut.bookHist(h,'resY_'+name[s]+'S'+str(s*10+l),'residual Y '+str(s*10+l),100,-20.,20.)
398 ut.bookHist(h,'track_'+name[s]+str(s*10+l),'track x/y '+str(s*10+l),80,-70.,10.,80,0.,80.)
399 for bar in range(10):
400 ut.bookHist(h,'nSiPMs_'+name[s]+str(s*10+l)+'_'+str(bar),'#sipms',16,-0.5,15.5)
401 ut.bookHist(h,'signalS_'+name[s]+str(s*10+l)+'_'+str(bar),'signal',100,0.,200.)
402 ut.bookHist(h,'signalL_'+name[s]+str(s*10+l)+'_'+str(bar),'signal',100,0.,200.)
403
404 if Nev < 0 : Nev = eventTree.GetEntries()
405 N=0
406 for event in eventTree:
407 N+=1
408 if N>Nev: break
409 theTrack = DS_track()
410 if not hasattr(theTrack,"getFittedState"): continue
411# now extrapolate to US and check for hits.
412 state = theTrack.getFittedState()
413 pos = state.getPos()
414 mom = state.getMom()
415
416 dsHits = {}
417 for s in range(2,4):
418 for p in range(5): dsHits[s*10+p]=[]
419 for aHit in eventTree.Digi_MuFilterHits:
420 if not aHit.isValid(): continue
421 plane = (aHit.GetDetectorID()//1000)
422 dsHits[plane].append(aHit)
423 s = 2
424 for plane in range(5):
425 z = zPos[s*10+plane]
426 lam = (z-pos.z())/mom.z()
427 xEx,yEx = pos.x()+lam*mom.x(),pos.y()+lam*mom.y()
428 # tag with station close by
429 if plane ==0: tag = 1
430 else: tag = plane -1
431 tagged = False
432 for aHit in dsHits[s*10+tag]:
433 detID = aHit.GetDetectorID()
434 MuFilter.GetPosition(detID,A,B)
435 dy = (A[1]+B[1])/2. - yEx
436 if abs(dy)<5: tagged = True
437 if not tagged: continue
438 rc = h['track_US'+str(s*10+plane)].Fill(xEx,yEx)
439 for aHit in dsHits[s*10+plane]:
440 detID = aHit.GetDetectorID()
441 bar = detID%1000
442 MuFilter.GetPosition(detID,A,B)
443 dy = (A[1]+B[1])/2. - yEx
444 rc = h['resY_US'+str(s*10+plane)].Fill(dy)
445 S = aHit.GetAllSignals()
446 # check for signal in left / right or small sipm
447 left,right,small = False,False,False
448 for x in S:
449 if x.first<8: left = True
450 else: right = True
451 if smallSiPMchannel(x.first): small = True
452 if left: rc = h['resY_USL'+str(s*10+plane)].Fill(dy)
453 if right: rc = h['resY_USR'+str(s*10+plane)].Fill(dy)
454 if small: rc = h['resY_USS'+str(s*10+plane)].Fill(dy)
455 if abs(dy)<2.0: # check channels
456 rc = h['nSiPMs_US'+str(s*10+plane)+'_'+str(bar)].Fill(S.size())
457 for x in S:
458 if smallSiPMchannel(x.first): rc = h['signalS_US'+str(s*10+plane)+'_'+str(bar)].Fill(x.second)
459 else: rc = h['signalL_US'+str(s*10+plane)+'_'+str(bar)].Fill(x.second)
460# look at delta time vs track X
461 dt = aHit.GetDeltaT(False)
462 h['dtLRvsX_US'+str(s*10+plane)].Fill(xEx,dt*1E9/freq)
463# check DS
464 s = 3
465 for plane in range(4):
466 z = zPos[s*10+plane]
467 lam = (z-pos.z())/mom.z()
468 xEx,yEx = pos.x()+lam*mom.x(),pos.y()+lam*mom.y()
469 for aHit in dsHits[s*10+plane]:
470 # only horizontal layers have two sides
471 if not aHit.GetnSides()==2: continue
472 detID = aHit.GetDetectorID()
473 MuFilter.GetPosition(detID,A,B)
474 dt = aHit.GetDeltaT(False)
475 h['dtLRvsX_DS'+str(s*10+plane)].Fill(xEx,dt*1E9/freq)
476
477 theTrack.Delete()
478
479latex = ROOT.TLatex()
480def analyze_USefficiency():
481 ut.bookCanvas(h,'E','',1200,2000,4,5)
482 ut.bookCanvas(h,'T','',800,1600,1,5)
483 ut.bookCanvas(h,'TDS','',800,900,1,2)
484 s=2
485 for plane in range(5):
486 tc = h['E'].cd(4*plane+1)
487 h['track_US'+str(s*10+plane)].Draw('colz')
488 l = 0
489 for side in ['L','R','S']:
490 tc = h['E'].cd(4*plane+2+l)
491 l+=1
492 res = h['resY_US'+side+str(s*10+plane)]
493 res.Draw()
494 binw = res.GetBinWidth(1)
495 myGauss = ROOT.TF1('gauss','abs([0])*'+str(binw)+'/(abs([2])*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+abs([3])',4)
496 myGauss.SetParameter(0,res.GetEntries())
497 myGauss.SetParameter(1,0)
498 myGauss.SetParameter(2,2.)
499 rc = res.Fit(myGauss,'SL','',-15.,15.)
500 fitResult = rc.Get()
501 tc.Update()
502 stats = res.FindObject('stats')
503 stats.SetOptFit(111111)
504 stats.SetX1NDC(0.63)
505 stats.SetY1NDC(0.25)
506 stats.SetX2NDC(0.98)
507 stats.SetY2NDC(0.94)
508 tracks = h['track_US'+str(s*10+plane)].GetEntries()
509 eff = fitResult.Parameter(0)/tracks
510 effErr = fitResult.ParError(0)/tracks
511 latex.DrawLatexNDC(0.2,0.8,'eff=%5.2F+/-%5.2F%%'%(eff,effErr))
512 h['E'].Print('Eff'+'-run'+str(options.runNumber)+'.png')
513
514 latex.SetTextColor(ROOT.kRed)
515 for plane in range(5):
516 tc = h['T'].cd(plane+1)
517 hist = h['dtLRvsX_US'+str(s*10+plane)]
518 hist.SetStats(0)
519 hist.Draw('colz')
520# get time x correlation, X = m*dt + b
521 h['gdtLRvsX_US'+str(s*10+plane)] = ROOT.TGraph()
522 g = h['gdtLRvsX_US'+str(s*10+plane)]
523 xproj = hist.ProjectionX('tmpx')
524 for nx in range(1,hist.GetNbinsX()+1):
525 tmp = hist.ProjectionY('tmp',nx,nx)
526 X = xproj.GetBinCenter(nx)
527 dt = tmp.GetMean()
528 g.SetPoint(nx-1,X,dt)
529 g.SetLineColor(ROOT.kBlue)
530 g.SetLineWidth(2)
531 g.Draw('same')
532 rc = g.Fit('pol1','S','',-50.,-10.)
533 result = rc.Get()
534 m = 1./result.Parameter(1)
535 b = -result.Parameter(0) * m
536 txt = 'dt X relation: X = #frac{dt}{%5.2F} +%5.2F '%(1./m,b)
537 latex.DrawLatexNDC(0.2,0.8,txt)
538 h['T'].Print('dTvsX'+'-run'+str(options.runNumber)+'.png')
539
540 # for DS
541 s = 3
542 for plane in range(2):
543 tc = h['TDS'].cd(plane+1)
544 hist = h['dtLRvsX_DS'+str(s*10+plane)]
545 hist.SetStats(0)
546 hist.Draw('colz')
547# get time x correlation, X = m*dt + b
548 h['gdtLRvsX_DS'+str(s*10+plane)] = ROOT.TGraph()
549 g = h['gdtLRvsX_DS'+str(s*10+plane)]
550 xproj = hist.ProjectionX('tmpx')
551 for nx in range(1,hist.GetNbinsX()+1):
552 tmp = hist.ProjectionY('tmp',nx,nx)
553 X = xproj.GetBinCenter(nx)
554 dt = tmp.GetMean()
555 g.SetPoint(nx-1,X,dt)
556 g.SetLineColor(ROOT.kBlue)
557 g.SetLineWidth(2)
558 g.Draw('same')
559 rc = g.Fit('pol1','S','',-50.,-10.)
560 result = rc.Get()
561 m = 1./result.Parameter(1)
562 b = -result.Parameter(0) * m
563 txt = 'dt X relation: X = #frac{dt}{%5.2F} +%5.2F '%(1./m,b)
564 latex.DrawLatexNDC(0.2,0.8,txt)
565 h['TDS'].Print('dTvsXDS'+'-run'+str(options.runNumber)+'.png')
566
567def timing(Nev = -1):
568 binning = {1:50,2:50,3:50}
569 for s in systemAndPlanes:
570 for l in range(systemAndPlanes[s]):
571 ut.bookHist(h,'timeDiffsL_'+str(s*10+l),' deviation from mean'+str(s*10+l)+'; [ns]',100,-binning[s],binning[s])
572 ut.bookHist(h,'timeDiffsR_'+str(s*10+l),' deviation from mean'+str(s*10+l)+'; [ns]',100,-binning[s],binning[s])
573 ut.bookHist(h,'timeDiffsLR_'+str(s*10+l),' mean left - mean right'+str(s*10+l)+'; [ns]',100,-binning[s],binning[s])
574
575 N=-1
576 if Nev < 0 : Nev = eventTree.GetEntries()
577 for event in eventTree:
578 N+=1
579 if N>Nev: break
580 for aHit in event.Digi_MuFilterHits:
581 if not aHit.isValid(): continue
582 detID = aHit.GetDetectorID()
583 s = detID//10000
584 l = (detID%10000)//1000 # plane number
585 bar = (detID%1000)
586 if s>2:
587 l=2*l
588 if bar>59:
589 bar=bar-60
590 l+=1
591 nSides = aHit.GetnSides()
592 nSiPMs = aHit.GetnSiPMs()
593
594 allChannels = aHit.GetAllTimes()
595 meanL, meanR = 0,0
596 nL,nR = 0,0
597 for c in allChannels:
598 if nSiPMs > c[0]: # left side
599 nL+=1
600 meanL+=c[1]
601 else: # right side
602 nR+=1
603 meanR+=c[1]
604 if nL >0: meanL = meanL/nL * 1E9/freq
605 if nR >0: meanR = meanR/nR * 1E9/freq
606 if nL>0 and nR>0: rc = h['timeDiffsLR_'+str(s*10+l)].Fill(meanL-meanR)
607 for c in allChannels:
608 if nSiPMs > c[0]: # left side
609 dt = c[1]* 1E9/freq - meanL
610 rc = h['timeDiffsL_'+str(s*10+l)].Fill(dt)
611 else:
612 dt = c[1]* 1E9/freq - meanR
613 rc = h['timeDiffsR_'+str(s*10+l)].Fill(dt)
614
615 installed_stations = {0:0,1:5,2:4}
616 x = 0
617 y = 0
618 for s in installed_stations:
619 if installed_stations[s]>0: x+=1
620 if installed_stations[s]>y: y = installed_stations[s]
621 ut.bookCanvas(h,'dt',' ',800,1600,1,5)
622 ut.bookCanvas(h,'dtLR',' ',1200,1600,x,y)
623
624 for S in installed_stations:
625 for l in range(installed_stations[S]):
626 if S==1:
627 n = 1+ l
628 tc = h['dt'].cd(n)
629 tc.SetLogy(1)
630 h['timeDiffsL_'+str(S+1)+str(l)].SetLineColor(ROOT.kRed)
631 h['timeDiffsL_'+str(S+1)+str(l)].Draw()
632 h['timeDiffsR_'+str(S+1)+str(l)].SetLineColor(ROOT.kGreen)
633 h['timeDiffsR_'+str(S+1)+str(l)].Draw('same')
634 n = S + l*x
635 tc = h['dtLR'].cd(n)
636 h['timeDiffsLR_'+str(S+1)+str(l)].Draw()
637 h['dt'].Print('dt'+'-run'+str(options.runNumber)+'.png')
638 h['dtLR'].Print('dtLR'+'-run'+str(options.runNumber)+'.png')
639
640
641
MuFilter::GetPosition
void GetPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition MuFilter.cxx:639
mufiHitMaps.smallSiPMchannel
smallSiPMchannel(i)
Definition mufiHitMaps.py:42
mufiHitMaps.hitMaps
hitMaps(Nev=-1)
Definition mufiHitMaps.py:61
mufiHitMaps.timing
timing(Nev=-1)
Definition mufiHitMaps.py:567
mufiHitMaps.USEfficiency
USEfficiency(Nev=-1)
Definition mufiHitMaps.py:387
mufiHitMaps.eventTime
eventTime(Nev=-1)
Definition mufiHitMaps.py:287