SND@LHC Software
Loading...
Searching...
No Matches
MufiCTR.py
Go to the documentation of this file.
1import ROOT,os,sys
2import rootUtils as ut
3
4from argparse import ArgumentParser
5parser = ArgumentParser()
6parser.add_argument("-f", dest="file", help="input offline monitoring file",default='/eos/experiment/sndlhc/www/offline/run004705.root')
7options = parser.parse_args()
8
9f = ROOT.TFile(options.file)
10ROOT.gROOT.cd()
11h={}
12S = {1:[1800,800,2,1],2:[1800,1500,2,3],3:[1800,1800,2,4]}
13channelsPerPlane = {1:7*16,3:60*2}
14Planes = {1:2,3:7}
15
16runN=options.file[options.file.find(".root")-6:options.file.find(".root")]
17
18def start():
19 for x in ['dTScifiDS','dTcorScifiDS']:
20 for d in ['','B2noB1','B1only']:
21 dd=''
22 if d!='': dd = d+'/'
23 # safety net if B2noB1 and B1Only plots don't exist
24 try:
25 h['mufi-'+x+d] = f.mufilter.Get(dd+'mufi-'+x+d).Clone('mufi-'+x+d)
26 except:
27 try:
28 #2024: changed structure of offline monitoring file with shifter and expert subdirs
29 h['mufi-'+x+d] = f.mufilter.Get('expert/'+dd+'mufi-'+x+d).Clone('mufi-'+x+d)
30 except:
31 continue
32 h['mufi-'+x+d].Draw()
33 h['mufi-'+x+d].Print('mufi-'+x+d+'.png')
34 for pad in h['mufi-'+x+d].GetListOfPrimitives():
35 for y in pad.GetListOfPrimitives():
36 if not y.ClassName().find('TH')<0:
37 hname = y.GetName()
38 h[hname] = y.Clone(hname)
39
40def start0(plist=['tmp4705p9'],s=3):
41 detector = 'mufi-'
42 wanted = []
43 system = 'DS'
44 if s==1: system = 'Veto'
45 for x in ['','cor']:
46 for d in ['','B2noB1','B1only']:
47 ut.bookCanvas(h,detector+'dT'+x+'Scifi'+system+d,'dt rel to scifi '+system,S[s][0],S[s][1],S[s][2],S[s][3])
48 for l in range(7):
49 tag = str(10*s+l)+d
50 wanted.append(detector+'dT'+x+'_'+tag)
51 for p in plist:
52 ut.readHists(h,p,wanted=wanted)
53 for x in ['','cor']:
54 for d in ['','B2noB1','B1only']:
55 for l in range(S[s][2]*S[s][3]):
56 if s==3 and l>6: continue
57 n = l+1
58 tag = str(10*s+l)+d
59 if s==3 and n==7: n=8
60 # safety net if B2noB1 and B1Only plots don't exist
61 try:
62 tc = h[detector+'dT'+x+'Scifi'+system+d].cd(n)
63 h[detector+'dT'+x+'_'+tag].Draw('colz')
64 except:
65 continue
66 h[detector+'dT'+x+'Scifi'+system+d].Print('mufi-dT'+x+'Scifi'+system+d+'.png')
67
68def xCheck():
69 ut.bookCanvas(h,'xCheckH',' ',1200,600,1,1)
70 ut.bookCanvas(h,'xCheckV',' ',1200,600,1,1)
71 tcolors = [0,ROOT.kRed,ROOT.kRed,ROOT.kBlue,ROOT.kBlue,ROOT.kGreen,ROOT.kGreen,ROOT.kCyan,ROOT.kCyan]
72 h['meanAndSig'] = {}
73 for b in ['','B2noB1']:
74 j=0
75 try:
76 tc = h['mufi-dTcorScifiDS'+b]
77 except:
78 continue
79 h['meanAndSig'][b] = {}
80 xmin,xmax = -5.,5.
81 if b=='B2noB1': xmin,xmax = -20.,-8.
82 for pad in tc.GetListOfPrimitives():
83 j+=1
84 if j==7: j==8
85 tpad = tc.cd(j)
86 for x in pad.GetListOfPrimitives():
87 if not x.ClassName().find('TH')<0:
88 if j%2==1: tcX = h['xCheckH'].cd()
89 else: tcX = h['xCheckV'].cd()
90 hname = x.GetName()
91 h[hname] = x.Clone(hname)
92 h[hname+'B1'] = h[hname].ProjectionX(hname+'B1')
93 h[hname+'B1'].SetStats(0)
94 h[hname+'B1'].SetLineColor(tcolors[j])
95 h[hname+'B1'].SetMarkerColor(tcolors[j])
96 h[hname+'B1'].SetMarkerStyle(30+j)
97 rc = h[hname+'B1'].Fit('gaus','SQ','',xmin,xmax)
98 fitRes = rc.Get()
99 h[hname+'B1'].GetFunction('gaus').SetLineColor(tcolors[j])
100 h['meanAndSig'][b][j] = [fitRes.Parameter(1),fitRes.Parameter(2)]
101 plots = {'H':[0,2,4],'V':[1,3,5,6]}
102 b=''
103 txt = ROOT.TLatex()
104 txt.SetTextFont(42)
105 txt.SetTextSize(0.04)
106 for c in plots:
107 tcX = h['xCheck'+c].cd()
108 tcX.SetLogy(1)
109 j=0
110 for x in plots[c]:
111 if x<2: h['mufi-dTcor_3'+str(x)+'B1'].Draw()
112 else: h['mufi-dTcor_3'+str(x)+'B1'].Draw('same')
113 if x==6: m,s = h['meanAndSig'][b][x+2]
114 else: m,s = h['meanAndSig'][b][x+1]
115 h['txt'+str(x)] = txt.DrawLatexNDC(0.15,0.85-j*0.05,"DS%i%s m=%5.2Fns #sigma=%5.2Fns"%(x//2,c,m,s))
116 h['txt'+str(x)].SetTextColor(h['mufi-dTcor_3'+str(x)+'B1'].GetLineColor())
117 j+=1
118 j+=1
119 try:
120 if c=='H':
121 h['txtB2H'] = txt.DrawLatexNDC(0.15,0.85-j*0.05,"BW tracks DS0H: %5.2Fns DS1H: %5.2Fns DS2H: %5.2Fns "%(
122 h['meanAndSig']['B2noB1'][1][0],h['meanAndSig']['B2noB1'][3][0],h['meanAndSig']['B2noB1'][5][0]))
123 h['txtB2H'].SetTextColor(ROOT.kMagenta)
124 else:
125 h['txtB2V'] = txt.DrawLatexNDC(0.15,0.85-j*0.05,"BW tracks DS0V: %5.2Fns DS1V: %5.2Fns DS2V: %5.2Fns DS3V: %5.2Fns "%(
126 h['meanAndSig']['B2noB1'][2][0],h['meanAndSig']['B2noB1'][4][0],h['meanAndSig']['B2noB1'][6][0],h['meanAndSig']['B2noB1'][8][0]))
127 h['txtB2V'].SetTextColor(ROOT.kMagenta)
128 except: pass
129 h['xCheck'+c].Update()
130 h['xCheck'+c].Print('mufi-dTScifiDS-xCheck'+c+'.png')
131
132
133par0={}
134alignTpar = {}
135alignTparB2 = {}
136
137def execute(s=3):
138 detector = 'mufi-'
139 system = 'DS'
140 if s==1: system = 'Veto'
141 ut.bookCanvas(h,'c1','',640,480,1,1)
142 h['c1'].cd()
143 fitLimits ={1:[-15,-5],11:[-15,-5],3:[-10,-2],13:[-25,-15]}
144 for l in range(S[s][2]*S[s][3]):
145 if s==3 and l>6: continue
146 tag = str(s*10+l)
147 hist = h[detector+'dT_'+tag]
148 alignTpar[s*10+l] = {}
149 alignTparB2[s*10+l] = {}
150 for i in range(hist.GetNbinsY()):
151 if s==1 and i > 8*16: continue
152 tagi = str(i*1000+s*10+l)
153 alignTpar[s*10+l][i] = [100,-100]
154 h['tmp'+tagi] = hist.ProjectionX('tmp'+tagi,i+1,i+1)
155 tmp = h['tmp'+tagi]
156 rc = tmp.Fit('gaus','SQ','',fitLimits[s][0],fitLimits[s][1])
157 fitres = rc.Get()
158 if fitres:
159 alignTpar[s*10+l][i]=[fitres.Parameter(1),fitres.Parameter(2)]
160 alignTparB2[s*10+l][i] = [100,-100]
161 try:
162 tmp = h[detector+'dT_'+tag+'B2noB1'].ProjectionX('tmp',i+1,i+1)
163 rc = tmp.Fit('gaus','SQ','',fitLimits[s+10][0],fitLimits[s+10][1])
164 fitres = rc.Get()
165 if fitres:
166 alignTparB2[s*10+l][i]=[fitres.Parameter(1),fitres.Parameter(2)]
167 except: pass
168 nbins = channelsPerPlane[s] * Planes[s]
169 ut.bookHist(h,'talign','time alignment B1;channel;offset [ns]',nbins,-0.5,nbins-0.5)
170 ut.bookHist(h,'talign2','time alignment B2;channel;offset [ns]',nbins,-0.5,nbins-0.5)
171 ut.bookHist(h,'tres','time diff res;channel;sigma [ns]',nbins,-0.5,nbins-0.5)
172 ut.bookHist(h,'t12','time diff B1 B2;channel;sigma [ns]',nbins,-0.5,nbins-0.5)
173 for l in range(S[s][2]*S[s][3]):
174 if s==3 and l>6: continue
175 for i in range(hist.GetNbinsY()):
176 k = l*channelsPerPlane[s]
177 rc = h['talign'].Fill(k+i,alignTpar[s*10+l][i][0])
178 rc = h['talign2'].Fill(k+i,alignTparB2[s*10+l][i][0])
179 rc = h['tres'].Fill(k+i,abs(alignTpar[s*10+l][i][1]))
180 if alignTpar[s*10+l][i][0]<99 and alignTparB2[s*10+l][i][0]<99:
181 rc = h['t12'].Fill(k+i,alignTparB2[s*10+l][i][0]-alignTpar[s*10+l][i][0])
182 else: rc = h['t12'].Fill(k+i,100)
183 if s==3:
184 h['talign'].SetMinimum(-10)
185 h['talign'].SetMaximum(1)
186 else:
187 h['talign'].SetMinimum(-15)
188 h['talign'].SetMaximum(-5)
189 h['talign2'].SetMinimum(-25)
190 h['talign2'].SetMaximum(1)
191 h['tres'].SetMinimum(0)
192 h['tres'].SetMaximum(2.5)
193 h['talign'].SetMarkerStyle(30)
194 h['talign'].SetMarkerColor(ROOT.kRed)
195 h['talign2'].SetMarkerStyle(32)
196 h['talign2'].SetMarkerColor(ROOT.kBlue)
197 h['tres'].SetMarkerStyle(31)
198 h['tres'].SetMarkerColor(ROOT.kGreen)
199 h['tres'].SetStats(0)
200 h['talign'].SetStats(0)
201 h['talign2'].SetStats(0)
202 ut.bookCanvas(h,'results',system+' time alignment',2400,1200,1,2)
203 tc = h['results'].cd(1)
204 h['talign'].Draw('phist')
205 tc = h['results'].cd(2)
206 h['tres'].Draw('phist')
207 h['results'].Print(system+'timeAlign-run'+runN+'.png')
208 h['t12'].SetMinimum(-30)
209 h['t12'].SetMaximum(0)
210 h['t12'].SetStats(0)
211 h['t12'].SetMarkerColor(ROOT.kMagenta)
212 h['t12'].SetMarkerStyle(53)
213
214 ut.bookCanvas(h,'results2','time diff between FW and BW tracks',1200,600,1,1)
215 h ['results2'].cd()
216 h['t12'].Draw('phist')
217 h['results2'].Print(system+'timeAlignFWBW-run'+runN+'.png')
218
219 for k in range(0,4):
220 h['talignReorder'+str(k)] = h['talign'].Clone('talignReorder'+str(k))
221 h['talignReorder'+str(k)].Reset()
222 if s==3:
223 slope = 0.0833
224 for i in range(h['talign'].GetNbinsX()):
225 if i<120 and i%2==0: p=0
226 elif i<120 and i%2==1: p=1
227 elif i<240 and i%2==0: p=2
228 elif i<360 and i%2==0: p=3
229 elif i<360 and i%2==1: p=4
230 elif i<480 and i%2==0: p=5
231 elif i<600 and i%2==0: p=6
232 elif i<600 and i%2==1: p=7
233 elif i<720 and i%2==0: p=8
234 elif i<840 and i%2==0: p=9
235 else: continue
236 if i%2==0: j=(i//2)%60
237 if i%2==1: j=((i-1)//2)%60
238 X = h['talign'].GetBinContent(i+1)
239 h['talignReorder0'].SetBinContent(j+1+p*60,X)
240 h['talignReorder1'].SetBinContent(j+1+p*60,X)
241 channel = (j+1+p*60)%60
242 if channel<30: X+= channel*slope - 15*slope
243 else: X-=(channel-30)*slope - 15*slope
244 h['talignReorder2'].SetBinContent(j+1+p*60,X)
245
246 for l in range(0,3):
247 c = 'talignReorder'+str(l)
248 ut.bookCanvas(h,'results'+c,system+' time alignment '+str(l),1200,600,1,1)
249 h['results'+c].SetGridy(1)
250 h['results'+c].cd()
251 par0[l]={}
252 for x in h[c].GetListOfFunctions(): x.Delete()
253 if l==2:
254 par1 = 0
255 planes = 20
256 iv = 30
257 h[c].SetMinimum(-7)
258 h[c].SetMaximum(-4)
259 else:
260 par1 = slope
261 planes = 20
262 iv = 30
263 h[c].SetMinimum(-20)
264 h[c].SetMaximum(3)
265 for p in range(planes):
266 if p==10:
267 imin = 0.+p*iv+7
268 imax = 0.+imin+iv-3
269 elif p==3:
270 imin = 0.+p*iv
271 imax = 0.+imin+iv-10
272 elif p==6:
273 imin = 0.+p*iv+5
274 imax = 0.+imin+iv-5
275 else:
276 imin = 0.+p*iv
277 imax = 0.+imin+iv
278 theFun = 'fun'+str(l*100+p)
279 h[theFun]=ROOT.TF1(theFun,'[0]+x*[1]')
280 sign = 1
281 if p%2==0: sign = -1
282 h[theFun].SetParameter(0, -6)
283 if l>0: h[theFun].FixParameter(1, sign*par1)
284 rc = h[c].Fit(h[theFun],'SQw+','',imin+2,imax-2)
285 res = rc.Get()
286 ix = 2
287 while res.Chi2() > 1.5:
288 ix+=2
289 if imax-ix < imin+2: 1/0
290 if h[c].GetFunction(theFun): h[c].GetFunction(theFun).Delete()
291 if p==3: rc = h[c].Fit(h[theFun],'SQw+','',imin+ix,imax-ix)
292 else : rc = h[c].Fit(h[theFun],'SQw+','',imin+2,imax-ix) # vertical no statistics
293 res = rc.Get()
294 par0[l][p]=[res.Parameter(0),res.Parameter(1)]
295 print("%i slope %5.3F b = %5.2F"%(p,res.Parameter(1),res.Parameter(0)))
296 h[c].GetXaxis().SetRange(1,600)
297 h[c].Draw()
298 h['results'+c].Print(system+'timeAlignFit'+str(l)+'-run'+runN+'.png')
299
300 c = 'talignReorder3'
301 ut.bookHist(h,'resT','t residuals',100,-3.,3.)
302 for i in range(h[c].GetNbinsX()):
303 p = i//30
304 if p<20:
305 res = h[c.replace('3','2')].GetBinContent(i+1)-par0[2][p][0]
306 h[c].SetBinContent(i+1,res)
307 rc = h['resT'].Fill(res)
308 ut.bookCanvas(h,'results'+c,system+' time alignment 3',1200,600,1,1)
309 h['results'+c].SetGridy(1)
310 h['results'+c].cd()
311 h[c].GetXaxis().SetRange(1,600)
312 h[c].SetMinimum(-1)
313 h[c].SetMaximum(1)
314 h[c].Draw()
315 h['results'+c].Print(system+'timeAlignFit3'+'-run'+runN+'.png')
316 h['c1'].cd()
317 rc = h['resT'].Fit('gaus','S','')
318 h['c1'].Update()
319 stats = h['resT'].FindObject('stats')
320 stats.SetOptFit(1)
321 stats.SetOptStat(1000000000)
322 stats.SetX1NDC(0.595611)
323 stats.SetY1NDC(0.614035)
324 stats.SetX2NDC(0.979624)
325 stats.SetY2NDC(0.936404)
326 h['c1'].Print(system+'timeResiudals'+'-run'+runN+'.png')
327
328 h['c1'].cd()
329 for I in [ [10.5,89.5],[249.5,328.5],[492.5,571.5] ]:
330 p+=1
331 rc = h['t12'].Fit('pol1','SQ','',I[0],I[1])
332 res = rc.Get()
333 print("%i delta t = %5.2F"%(p,res.Parameter(0)))
334 m=0
335 for i in par0[0]:
336 m+=abs(par0[0][i][1])
337 print('average slope = ',m/len(par0[0]))
338
339# golden nu event
340# $EOSSHIP/eos/experiment/sndlhc/convertedData/physics/2022/run_005056/sndsw_raw-0101.root -g geofile_sndlhc_TI18_V7_22November2022.root
341# 849600
342def timingOfEvent(makeCluster=False,debug=False):
343 firstScifi_z = 300*u.cm
344 TDC2ns = 1E9/160.316E6
345 ut.bookHist(h,'evTimeDS','cor time of hits;[ns]',100,-20.,20)
346 ut.bookHist(h,'evTimeScifi','cor time of hits blue DS red Scifi;[ns]',100,-20.,20)
347 ut.bookCanvas(h,'tevTime','cor time of hits',1024,768,1,1)
348 h['evTimeScifi'].SetLineColor(ROOT.kRed)
349 h['evTimeDS'].SetLineColor(ROOT.kBlue)
350 h['evTimeScifi'].SetStats(0)
351 h['evTimeDS'].SetStats(0)
352 if makeCluster: trackTask.scifiCluster()
353 meanXY = {}
354 for siCl in trackTask.clusScifi:
355 detID = siCl.GetFirst()
356 s = detID//1000000
357 isVertical = detID%1000000//100000
358 siCl.GetPosition(A,B)
359 z=(A[2]+B[2])/2.
360 pos = (A[1]+B[1])/2.
361 L = abs(A[0]-B[0])/2.
362 if isVertical:
363 pos = (A[0]+B[0])/2.
364 L = abs(A[1]-B[1])/2.
365 corTime = geo.modules['Scifi'].GetCorrectedTime(detID, siCl.GetTime(), L) - (z-firstScifi_z)/u.speedOfLight
366 h['evTimeScifi'].Fill(corTime)
367 if debug: print(detID,corTime,pos)
368 for aHit in eventTree.Digi_MuFilterHits:
369 detID = aHit.GetDetectorID()
370 if not detID//10000==3: continue
371 if aHit.isVertical(): nmax = 1
372 else: nmax=2
373 geo.modules['MuFilter'].GetPosition(detID,A,B)
374 z=(A[2]+B[2])/2.
375 pos = (A[1]+B[1])/2.
376 L = abs(A[0]-B[0])/2.
377 if isVertical:
378 pos = (A[0]+B[0])/2.
379 L = abs(A[1]-B[1])/2.
380 for i in range(nmax):
381 corTime = geo.modules['MuFilter'].GetCorrectedTime(detID, i, aHit.GetTime(i)*TDC2ns, L)- (z-firstScifi_z)/u.speedOfLight
382 h['evTimeDS'].Fill(corTime)
383 if debug: print(detID,i,corTime,pos)
384 tc=h['tevTime'].cd()
385 h['evTimeScifi'].Draw()
386 h['evTimeDS'].Draw('same')
387
388
389
MufiCTR.start0
start0(plist=['tmp4705p9'], s=3)
Definition MufiCTR.py:40
MufiCTR.execute
execute(s=3)
Definition MufiCTR.py:137
MufiCTR.start
start()
Definition MufiCTR.py:18
MufiCTR.timingOfEvent
timingOfEvent(makeCluster=False, debug=False)
Definition MufiCTR.py:342
MufiCTR.xCheck
xCheck()
Definition MufiCTR.py:68