SND@LHC Software
Loading...
Searching...
No Matches
Mufi_monitoring.py
Go to the documentation of this file.
1#!/usr/bin/env python
2import ROOT,os,sys
3import rootUtils as ut
4import shipunit as u
5import numpy as np
6
7A,B = ROOT.TVector3(),ROOT.TVector3()
8detector = "mufi-"
9# a dictionary to keep track of US sides with gel/no gel
10# L first, next R as usual
11US_gel_dict = {20:["no GEL", "with GEL"],
12 21:["with GEL", "no GEL"],
13 22:["no GEL", "with GEL"],
14 23:["with GEL", "no GEL"],
15 24:["with GEL", "no GEL"]}
16
17class Mufi_hitMaps(ROOT.FairTask):
18 " produce hitmaps for MuFilter, Veto/US/DS"
19 """
20 veto system 2 layers with 7 bars and 8 sipm channels on both ends
21 1 layer with 7 bars and 8 sipms on the top
22 US system 5 layers with 10 bars and 8 sipm channels on both ends
23 DS system horizontal(3) planes, 60 bars, readout on both sides, single channel
24 vertical(4) planes, 60 bar, readout on top, single channel
25 """
26 def Init(self,options,monitor):
27 self.M = monitor
28 sdict = self.M.sdict
29 h = self.M.h
30 run = ROOT.FairRunAna.Instance()
31 self.trackTask = self.M.trackTask
32 if not self.trackTask: self.trackTask = run.GetTask('houghTransform')
33 ioman = ROOT.FairRootManager.Instance()
34 self.OT = ioman.GetSink().GetOutTree()
35 self.mufi_vsignal = self.M.Scifi.GetConfParF("Scifi/signalSpeed")
36
37 channelsPerSystem = {1:self.M.MuFilter.GetConfParI("MuFilter/VetonSiPMs"),
38 2:self.M.MuFilter.GetConfParI("MuFilter/UpstreamnSiPMs"),
39 3:self.M.MuFilter.GetConfParI("MuFilter/DownstreamnSiPMs")}
40
41# type of crossing, check for b1only,b2nob1,nobeam
42 if self.M.fsdict or self.M.hasBunchInfo: self.xing = {'':True,'B1only':False,'B2noB1':False,'noBeam':False}
43 else: self.xing = {'':True}
44 for xi in self.xing:
45 ut.bookHist(h,detector+'Noise'+xi,'events with hits in single plane; s*10+l;',40,0.5,39.5)
46 for s in monitor.systemAndPlanes:
47 ut.bookHist(h,sdict[s]+'Mult'+xi,'QDCs vs nr hits; #hits; QDC [a.u.]',200,0.,800.,200,0.,300.)
48 for l in range(monitor.systemAndPlanes[s]):
49 ut.bookHist(h,detector+'hitmult_'+str(s*10+l)+xi,'hit mult / plane '+sdict[s]+str(l)+'; #hits',61,-0.5,60.5)
50 ut.bookHist(h,detector+'hit_'+str(s*10+l)+xi,'channel map / plane '+sdict[s]+str(l)+'; #channel',160,-0.5,159.5)
51 ut.bookHist(h,detector+'Xhit_'+str(s*10+l)+xi,'Xchannel map / plane '+sdict[s]+str(l)+'; #channel',160,-0.5,159.5)
52
53 # Only large SiPMs will be monitored.
54 # Small SiPMs suffer from large time offsets extending to other events!
55 note = ''
56 NSmallSiPMs = 0
57 if s==2:
58 note = ', large only'
59 NSmallSiPMs = 2
60 ut.bookHist(h,detector+'chanActiveRight_'+str(s*10+l),
61 sdict[s]+' '+str(l)+'R channel hit multiplicity; #channel; bar',
62 channelsPerSystem[s],-0.5,channelsPerSystem[s]-0.5,
63 monitor.systemAndBars[s],-0.5,monitor.systemAndBars[s]-0.5)
64 ut.bookHist(h,detector+'chanNfiredRight_'+str(s*10+l),
65 sdict[s]+' '+str(l)+'R number of fired channels '+note+'; N fired channels; bar',
66 channelsPerSystem[s]+1-NSmallSiPMs,-0.5,channelsPerSystem[s]+0.5-NSmallSiPMs,
67 monitor.systemAndBars[s],-0.5,monitor.systemAndBars[s]-0.5)
68 side = 'L'
69 if (s==1 and l==2) or (s==3 and (l%2==1 or l==6)):
70 side = 'T'
71 ut.bookHist(h,detector+'chanActiveLeft_'+str(s*10+l),
72 sdict[s]+' '+str(l)+side+' channel hit multiplicity; #channel; bar',
73 channelsPerSystem[s],-0.5,channelsPerSystem[s]-0.5,
74 monitor.systemAndBars[s],-0.5,monitor.systemAndBars[s]-0.5)
75 ut.bookHist(h,detector+'chanNfiredLeft_'+str(s*10+l),
76 sdict[s]+' '+str(l)+side+' number of fired channels '+note+'; N fired channels; bar',
77 channelsPerSystem[s]+1-NSmallSiPMs,-0.5,channelsPerSystem[s]+0.5-NSmallSiPMs,
78 monitor.systemAndBars[s],-0.5,monitor.systemAndBars[s]-0.5)
79
80 if s==3:
81 ut.bookHist(h,detector+'bar_'+str(s*10+l)+xi,'bar map / plane '+sdict[s]+str(l)+'; bar',60,-0.5,59.5)
82 ut.bookHist(h,detector+'dT_'+str(s*10+l)+xi,'dT with respect to first scifi '+sdict[s]+str(l)+'; dt [ns] ;# bar + channel', 100,-25.,5.,120,-0.5,2*60-0.5)
83 ut.bookHist(h,detector+'dTcor_'+str(s*10+l)+xi,'dTcor with respect to first scifi '+sdict[s]+str(l)+'; dt [ns] ;# bar + channel',100,-25.,5.,120,-0.5,2*60-0.5)
84 if l == 4:
85 for ss in range(1,6):
86 ut.bookHist(h,'deltaTScifiMufiHit_'+str(ss)+xi,'deltaT scifi earliest hit versus DS hit 2H',200,-25.,25.)
87 else:
88 ut.bookHist(h,detector+'bar_'+str(s*10+l)+xi,'bar map / plane '+sdict[s]+str(l)+'; bar',10,-0.5,9.5)
89 if s==1:
90 ut.bookHist(h,detector+'dT_'+str(s*10+l)+xi,'dT with respect to first scifi '+sdict[s]+str(l)+'; dt [ns] ;# bar + channel', 100,-25.,5.,120,-0.5,2*8*7-0.5)
91 ut.bookHist(h,detector+'dTcor_'+str(s*10+l)+xi,'dTcor with respect to first scifi '+sdict[s]+str(l)+'; dt [ns] ;# bar + channel',100,-25.,5.,120,-0.5,2*8*7-0.5)
92 ut.bookHist(h,detector+'sig_'+str(s*10+l)+xi,'signal / plane '+sdict[s]+str(l)+'; QDC [a.u.]',200,0.0,200.)
93 if s==2:
94 ut.bookHist(h,detector+'sigS_'+str(s*10+l)+xi,'signal / plane '+sdict[s]+str(l)+'; QDC [a.u.]',200,0.0,200.)
95 ut.bookHist(h,detector+'TsigS_'+str(s*10+l)+xi,'signal / plane '+sdict[s]+str(l)+'; QDC [a.u.]',200,0.0,200.)
96 tagL = US_gel_dict[s*10+l][0]+": "
97 tagR = US_gel_dict[s*10+l][1]+": "
98 else:
99 tagL = ""
100 tagR = ""
101 histo_title_helper = 'signal / plane '+sdict[s]+str(l)
102 ut.bookHist(h,detector+'sigL_'+str(s*10+l)+xi,tagL+histo_title_helper+'; QDC [a.u.]',200,0.0,200.)
103 ut.bookHist(h,detector+'sigR_'+str(s*10+l)+xi,tagR+histo_title_helper+'; QDC [a.u.]',200,0.0,200.)
104 ut.bookHist(h,detector+'Tsig_'+str(s*10+l)+xi,histo_title_helper+'; QDC [a.u.]',200,0.0,200.)
105 ut.bookHist(h,detector+'TsigL_'+str(s*10+l)+xi,tagL+histo_title_helper+'; QDC [a.u.]',200,0.0,200.)
106 ut.bookHist(h,detector+'TsigR_'+str(s*10+l)+xi,tagR+histo_title_helper+'; QDC [a.u.]',200,0.0,200.)
107 # not used currently?
108 ut.bookHist(h,detector+'occ_'+str(s*10+l)+xi,'channel occupancy '+sdict[s]+str(l),100,0.0,200.)
109 ut.bookHist(h,detector+'occTag_'+str(s*10+l)+xi,'channel occupancy '+sdict[s]+str(l),100,0.0,200.)
110
111 ut.bookHist(h,detector+'leftvsright_1'+xi,'Veto hits in left / right; Left: # hits; Right: # hits',10,-0.5,9.5,10,-0.5,9.5)
112 ut.bookHist(h,detector+'leftvsright_2'+xi,'US hits in left / right; L: # hits; R: # hits',10,-0.5,9.5,10,-0.5,9.5)
113 ut.bookHist(h,detector+'leftvsright_3'+xi,'DS hits in left / right; L: # hits; R: # hits',2,-0.5,1.5,2,-0.5,1.5)
114 ut.bookHist(h,detector+'leftvsright_signal_1'+xi,'Veto signal in left / right; Left: QDC [a.u.]; Right: QDC [a.u.]',100,-0.5,200.,100,-0.5,200.)
115 ut.bookHist(h,detector+'leftvsright_signal_2'+xi,'US signal in left / right; L: QDC [a.u.]; R: QDC [a.u.]',100,-0.5,200.,100,-0.5,200.)
116 ut.bookHist(h,detector+'leftvsright_signal_3'+xi,'DS signal in left / right; L: QDC [a.u.]; R: QDC [a.u.]',100,-0.5,200.,100,-0.5,200.)
117
118 ut.bookHist(h,detector+'dtime'+xi,'delta event time; dt [ns]',100,0.0,1000.)
119 ut.bookHist(h,detector+'dtimeu'+xi,'delta event time; dt [us]',100,0.0,1000.)
120 ut.bookHist(h,detector+'dtimem'+xi,'delta event time; dt [ms]',100,0.0,1000.)
121
122 ut.bookHist(h,detector+'bs'+xi,'beam spot; x[cm]; y[cm]',100,-100.,10.,100,0.,80.)
123 ut.bookHist(h,detector+'bsDS'+xi,'beam spot,#bar X, #bar Y',60,-0.5,59.5,60,-0.5,59.5)
124 ut.bookHist(h,detector+'slopes'+xi,'muon DS track slopes; slope X [rad]; slope Y [rad]',150,-1.5,1.5,150,-1.5,1.5)
125 ut.bookHist(h,detector+'trackPos'+xi,'muon DS track pos; x [cm]; y [cm]',100,-90,10.,80,0.,80.)
126 ut.bookHist(h,detector+'trackPosBeam'+xi,'beam track pos slopes<0.1rad; x [cm]; y [cm]',100,-90,10.,80,0.,80.)
127
128 for bar in range(monitor.systemAndBars[s]):
129 ut.bookHist(h,detector+'chanmult_'+str(s*1000+100*l+bar)+xi,'channels firing per bar '+sdict[s]+str(l)+" bar "+str(bar)+'; fired channels',20,-0.5,19.5)
130#
131 xmin = options.Mufixmin
132 xmax = -xmin
133 ut.bookHist(h,detector+'resX_'+sdict[s]+str(s*10+l)+xi,'residual X'+str(s*10+l)+'; [#cm]',
134 100,xmin,xmax,60,-60.,0.)
135 ut.bookHist(h,detector+'resY_'+sdict[s]+str(s*10+l)+xi,'residual Y'+str(s*10+l)+'; [#cm]',
136 100,xmin,xmax,70,2.,68.)
137
138 for x in h:
139 if isinstance(h[x], ROOT.TH2):
140 h[x].SetStats(0)
141
142 self.listOfHits = {1:[],2:[],3:[]}
143 def ExecuteEvent(self,event):
144 systemAndPlanes =self.M.systemAndPlanes
145 sdict = self.M.sdict
146 h = self.M.h
147 mult = {}
148 planes = {}
149 for i in self.listOfHits: self.listOfHits[i].clear()
150 for s in systemAndPlanes:
151 for l in range(systemAndPlanes[s]): mult[s*10+l]=0
152
153 self.beamSpot(event)
154 withDSTrack = False
155 for aTrack in self.M.Reco_MuonTracks:
156 if aTrack.GetUniqueID()==3: withDSTrack = True
157
158 for aHit in event.Digi_MuFilterHits:
159 Minfo = self.M.MuFilter_PlaneBars(aHit.GetDetectorID())
160 s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
161 nSiPMs = aHit.GetnSiPMs()
162 nSides = aHit.GetnSides()
163 for c in aHit.GetAllSignals(False,False):
164 if aHit.isMasked(c.first):
165 channel = bar*nSiPMs*nSides + c.first
166 self.M.fillHist1(detector+'Xhit_'+str(s)+str(l),channel)
167
168 if not aHit.isValid(): continue
169 mult[s*10+l]+=1
170 key = s*100+l
171 if not key in planes: planes[key] = {}
172 sumSignal = self.M.map2Dict(aHit,'SumOfSignals')
173 planes[key][bar] = [sumSignal['SumL'],sumSignal['SumR']]
174# check left/right
175 allChannels = self.M.map2Dict(aHit,'GetAllSignals')
176 for c in allChannels:
177 self.listOfHits[s].append(allChannels[c])
178 Nleft,Nright,Sleft,Sright = 0,0,0,0
179 # count the small SiPMs seperately
180 NSmallLeft, NSmallRight = 0,0
181 for c in allChannels:
182 if nSiPMs > c: # left side
183 Nleft+=1
184 if s==2 and (c==2 or c==5): NSmallLeft+=1
185 Sleft+=allChannels[c]
186 h[detector+'chanActiveLeft_'+str(s*10+l)].Fill(c, bar)
187 else:
188 Nright+=1
189 if s==2 and (c==10 or c==13): NSmallRight+=1
190 Sright+=allChannels[c]
191 h[detector+'chanActiveRight_'+str(s*10+l)].Fill(c-nSiPMs, bar)
192 self.M.fillHist1(detector+'chanmult_'+str(s*1000+100*l+bar),Nleft)
193 self.M.fillHist1(detector+'chanmult_'+str(s*1000+100*l+bar),10+Nright)
194 h[detector+'chanNfiredLeft_'+str(s*10+l)].Fill(Nleft-NSmallLeft, bar)
195 h[detector+'chanNfiredRight_'+str(s*10+l)].Fill(Nright-NSmallRight, bar)
196 if not aHit.isVertical(): # vertical DS plane is read out only on one side
197 self.M.fillHist2(detector+'leftvsright_'+str(s),Nleft,Nright)
198 self.M.fillHist2(detector+'leftvsright_signal_'+str(s),Sleft,Sright)
199#
200 for c in allChannels:
201 channel = bar*nSiPMs*nSides + c
202 self.M.fillHist1(detector+'hit_'+str(s)+str(l),int(channel))
203 self.M.fillHist1(detector+'bar_'+str(s)+str(l),bar)
204 if s==2 and self.M.smallSiPMchannel(c) :
205 self.M.fillHist1(detector+'sigS_'+str(s)+str(l),allChannels[c])
206 if withDSTrack: self.M.fillHist1(detector+'TsigS_'+str(s)+str(l),allChannels[c])
207 elif c<nSiPMs:
208 self.M.fillHist1(detector+'sigL_'+str(s)+str(l),allChannels[c])
209 if withDSTrack: self.M.fillHist1(detector+'TsigL_'+str(s)+str(l),allChannels[c])
210 else :
211 self.M.fillHist1(detector+'sigR_'+str(s)+str(l),allChannels[c])
212 if withDSTrack: self.M.fillHist1(detector+'sigR_'+str(s)+str(l),allChannels[c])
213 self.M.fillHist1(detector+'sig_'+str(s)+str(l),allChannels[c])
214 if withDSTrack: self.M.fillHist1(detector+'sig_'+str(s)+str(l),allChannels[c])
215 allChannels.clear()
216#
217 # noise event with many hits in one plane
218 onePlane = []
219 for x in mult:
220 if mult[x]>3: onePlane.append(x)
221 if len(onePlane)==1:
222 self.M.fillHist1(detector+'Noise',onePlane[0])
223
224#
225 for s in self.listOfHits:
226 nhits = len(self.listOfHits[s])
227 qcdsum = 0
228 for i in range(nhits):
229 self.M.fillHist2(sdict[s]+'Mult',nhits, self.listOfHits[s][i])
230 for s in systemAndPlanes:
231 for l in range(systemAndPlanes[s]):
232 self.M.fillHist1(detector+'hitmult_'+str(s*10+l),mult[s*10+l])
233# mufi residuals with scifi tracks
234 for aTrack in self.M.Reco_MuonTracks:
235 if not aTrack.GetUniqueID()==1: continue
236 fitStatus = aTrack.getFitStatus()
237 if not fitStatus.isFitConverged(): continue
238 posMom = {}
239 fstate = aTrack.getFittedState()
240 posMom['first'] = [fstate.getPos(),fstate.getMom()]
241 # fstate = aTrack.getFittedState(aTrack.getNumPointsWithMeasurement()-1) does not make a difference
242 posMom['last'] = [fstate.getPos(),fstate.getMom()]
243 rc = self.trackTask.trackDir(aTrack)
244 scifi_time0 = rc[2]
245 pos,mom = posMom['first']
246 lam = (self.trackTask.firstScifi_z-pos.z())/mom.z()
247 # nominal first position
248 pos1 = ROOT.TVector3(pos.x()+lam*mom.x(),pos.y()+lam*mom.y(),self.trackTask.firstScifi_z)
249 dsHitTimes = []
250 for aHit in event.Digi_MuFilterHits:
251 if not aHit.isValid(): continue
252 detID = aHit.GetDetectorID()
253 Minfo = self.M.MuFilter_PlaneBars(detID)
254 s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
255 self.M.MuFilter.GetPosition(detID,A,B)
256# calculate DOCA
257 if s==1: pos,mom = posMom['first']
258 else: pos,mom = posMom['last']
259 zEx = self.M.zPos['MuFilter'][s*10+l]
260 lam = (zEx-pos.z())/mom.z()
261 xEx,yEx = pos.x()+lam*mom.x(),pos.y()+lam*mom.y()
262 pq = A-pos
263 uCrossv= (B-A).Cross(mom)
264 doca = pq.Dot(uCrossv)/uCrossv.Mag()
265 self.M.fillHist2(detector+'resX_'+sdict[s]+str(s*10+l),doca/u.cm,xEx)
266 self.M.fillHist2(detector+'resY_'+sdict[s]+str(s*10+l),doca/u.cm,yEx)
267# calculate time difference for DS
268 if (s==3 and abs(doca)<2.5*u.cm) or (s==1 and abs(doca)<6*u.cm):
269 # horizontal layers have left and right sipms
270 if aHit.isVertical(): nmax = 1
271 else: nmax = 2
272 barMult = 2
273 if s==1:
274 nmax = 16
275 barMult = 16
276 for i in range(nmax):
277 if aHit.GetTime(i) < 0: continue # not valid time
278 posM = ROOT.TVector3(xEx,yEx,zEx)
279 # correct for flight length
280 trajLength = (posM-pos1).Mag()
281 # correct for signal speed, need to know left or right
282 if s==3:
283 if i==1: X = B-posM # B is right only horizontal planes have a second readout
284 else: X = A-posM # A is on the left, or top for vertical planes
285 if s==1:
286 if i<8: X = A-posM
287 else: X = B-posM
288 L = X.Mag()/self.mufi_vsignal
289 tM = aHit.GetTime(i)*self.M.TDC2ns - L - trajLength/u.speedOfLight
290 self.M.fillHist2(detector+'dT_'+str(s*10+l),tM-scifi_time0,bar*barMult+i)
291 # use corrected time
292 if s==3:
293 corTime = self.M.MuFilter.GetCorrectedTime(detID, i, aHit.GetTime(i)*self.M.TDC2ns, X.Mag())
294 tM = corTime - trajLength/u.speedOfLight
295 self.M.fillHist2(detector+'dTcor_'+str(s*10+l),tM-scifi_time0,bar*barMult+i)
296 if s==3 and l==2:
297 timeLeft = aHit.GetTime(0)
298 timeRight = aHit.GetTime(1)
299 if timeLeft>0 and timeRight>0:
300 dL = abs(A[0]-B[0])
301 avTime = self.M.MuFilter.GetCorrectedTime(detID, 0, timeLeft*self.M.TDC2ns,0) + \
302 self.M.MuFilter.GetCorrectedTime(detID, 1, timeRight*self.M.TDC2ns,0)
303 dsHitTimes.append( (avTime-abs(A[0]-B[0])/15)/2) # L/2 / 15cm/ns
304# fill histograms with time difference of earliest scifi hit in station i and last horizontal DS time, average left and right
305 if len(dsHitTimes)>0:
306 dsHitTimes.sort()
307 scifiHitTimes = {1:[],2:[],3:[],4:[],5:[]}
308 for scifiHit in event.Digi_ScifiHits:
309 detID = scifiHit.GetDetectorID()
310 s = int(scifiHit.GetDetectorID()/1000000)
311 scifiHitTimes[s].append(self.M.Scifi.GetCorrectedTime(detID,scifiHit.GetTime()*self.M.TDC2ns,0))
312 for s in scifiHitTimes:
313 if len(scifiHitTimes[s])<1: continue
314 scifiHitTimes[s].sort()
315 deltaT = dsHitTimes[0] - scifiHitTimes[s][0] - (self.M.zPos['MuFilter'][34]-self.M.zPos['Scifi'][s*10])/u.speedOfLight
316 self.M.fillHist1('deltaTScifiMufiHit_'+str(s),deltaT)
317 def beamSpot(self,event):
318 if not self.trackTask: return
319 h = self.M.h
320 W = self.M.Weight
321 Xbar = -10
322 Ybar = -10
323 for aTrack in self.M.Reco_MuonTracks:
324 if not aTrack.GetUniqueID()==3: continue
325 state = aTrack.getFittedState()
326 pos = state.getPos()
327 rc = h[detector+'bs'].Fill(pos.x(),pos.y(),W)
328 mom = state.getMom()
329 slopeX= mom.X()/mom.Z()
330 slopeY= mom.Y()/mom.Z()
331 pos = state.getPos()
332
333 self.M.fillHist2(detector+'slopes',slopeX,slopeY)
334 self.M.fillHist2(detector+'trackPos',pos.X(),pos.Y())
335 if abs(slopeX)<0.1 and abs(slopeY)<0.1: self.M.fillHist2(detector+'trackPosBeam',pos.X(),pos.Y())
336 if not Ybar<0 and not Xbar<0 and abs(slopeY)<0.01: self.M.fillHist2(detector+'bsDS',Xbar,Ybar)
337
338 def Plot(self):
339 h = self.M.h
340 sdict = self.M.sdict
341 systemAndPlanes =self.M.systemAndPlanes
342 S = {1:[1800,800,systemAndPlanes[1],1],2:[1800,1500,2,3],3:[1800,1800,2,4]}
343 for xi in self.xing:
344 if not self.M.fsdict and not self.M.hasBunchInfo and xi!='': continue
345
346 for s in S:
347 ut.bookCanvas(h,detector+'hitmaps' +sdict[s]+xi,'hitmaps' +sdict[s],S[s][0],S[s][1],S[s][2],S[s][3])
348 ut.bookCanvas(h,detector+'Xhitmaps' +sdict[s]+xi,'Xhitmaps' +sdict[s],S[s][0],S[s][1],S[s][2],S[s][3])
349 ut.bookCanvas(h,detector+'barmaps'+sdict[s]+xi,'barmaps'+sdict[s],S[s][0],S[s][1],S[s][2],S[s][3])
350 if s==3 or s==1:
351 ut.bookCanvas(h,detector+'dTScifi'+sdict[s]+xi,'dt rel to scifi'+sdict[s],S[s][0],S[s][1],S[s][2],S[s][3])
352 ut.bookCanvas(h,detector+'dTcorScifi'+sdict[s]+xi,'dtcor rel to scifi'+sdict[s],S[s][0],S[s][1],S[s][2],S[s][3])
353
354 for l in range(systemAndPlanes[s]):
355 n = l+1
356 if s==3 and n==7: n=8
357 tc = h[detector+'hitmaps'+sdict[s]+xi].cd(n)
358 tag = str(s)+str(l)+xi
359 h[detector+'hit_'+tag].Draw()
360 tc = h[detector+'Xhitmaps'+sdict[s]+xi].cd(n)
361 h[detector+'Xhit_'+tag].Draw()
362
363 tc = h[detector+'barmaps'+sdict[s]+xi].cd(n)
364 h[detector+'bar_'+tag].Draw()
365 if s==3 or s==1:
366 tc = h[detector+'dTScifi'+sdict[s]+xi].cd(n)
367 h[detector+'dT_'+tag].Draw('colz')
368 tc = h[detector+'dTcorScifi'+sdict[s]+xi].cd(n)
369 h[detector+'dTcor_'+tag].Draw('colz')
370
371 ut.bookCanvas(h,detector+'hitmult'+xi,'hit multiplicities per plane',2000,1600,4,3)
372 k=1
373 for s in systemAndPlanes:
374 for l in range(systemAndPlanes[s]):
375 tc = h[detector+'hitmult'+xi].cd(k)
376 tc.SetLogy(1)
377 k+=1
378 rc = h[detector+'hitmult_'+str(s*10+l)+xi].Draw()
379 ut.bookCanvas(h,'noise'+xi,' ',1200,1800,1,1)
380 tc = h['noise'+xi].cd()
381 h[detector+'Noise'+xi].Draw()
382
383 ut.bookCanvas(h,'VETO'+xi,' ',1200,1800,1,2)
384 for l in range(2):
385 tc = h['VETO'+xi].cd(l+1)
386 hname = detector+'hit_'+str(1)+str(l)+xi
387 h[hname].SetStats(0)
388 h[hname].Draw()
389 for n in range(7):
390 x = (n+1)*16-0.5
391 y = h[detector+'hit_'+str(1)+str(l)+xi].GetMaximum()
392 lname = 'L'+str(n)+hname
393 h[lname] = ROOT.TLine(x,0,x,y)
394 h[lname].SetLineColor(ROOT.kRed)
395 h[lname].SetLineStyle(9)
396 h[lname].Draw('same')
397
398 ut.bookCanvas(h,'USBars'+xi,' ',1200,900,1,1)
399 colours = {0:ROOT.kOrange,1:ROOT.kRed,2:ROOT.kGreen,3:ROOT.kBlue,4:ROOT.kMagenta,5:ROOT.kCyan,
400 6:ROOT.kAzure,7:ROOT.kPink,8:ROOT.kSpring}
401 for i in range(5):
402 h[detector+'bar_2'+str(i)+xi].SetLineColor(colours[i])
403 h[detector+'bar_2'+str(i)+xi].SetLineWidth(2)
404 h[detector+'bar_2'+str(i)+xi].SetStats(0)
405 h[detector+'bar_20'+xi].Draw()
406 h[detector+'bar_21'+xi].Draw('same')
407 h[detector+'bar_22'+xi].Draw('same')
408 h[detector+'bar_23'+xi].Draw('same')
409 h[detector+'bar_24'+xi].Draw('same')
410 h[detector+'lbar2'+xi]=ROOT.TLegend(0.6,0.6,0.99,0.99)
411 for i in range(5):
412 h[detector+'lbar2'+xi].AddEntry(h[detector+'bar_2'+str(i)+xi],'plane '+str(i+1),"f")
413 h[detector+'lbar2'+xi].Draw()
414 for i in range(7):
415 h[detector+'hit_3'+str(i)+xi].SetLineColor(colours[i])
416 h[detector+'hit_3'+str(i)+xi].SetLineWidth(2)
417 h[detector+'hit_3'+str(i)+xi].SetStats(0)
418 h[detector+'hit_30'+xi].Draw()
419 for i in range(1,7):
420 h[detector+'hit_3'+str(i)+xi].Draw('same')
421 h[detector+'lbar3'+xi]=ROOT.TLegend(0.6,0.6,0.99,0.99)
422 for i in range(7):
423 h[detector+'lbar3'+xi].AddEntry(h[detector+'hit_3'+str(i)+xi],'plane '+str(i+1),"f")
424 h[detector+'lbar3'+xi].Draw()
425
426 ut.bookCanvas(h,detector+'LR'+xi,' ',1800,900,3,2)
427 for i in range(1,4):
428 h[detector+'LR'+xi].cd(i)
429 h[detector+'leftvsright_'+str(i)+xi].Draw('textBox')
430 h[detector+'LR'+xi].cd(i+3)
431 h[detector+'leftvsright_signal_'+str(i)+xi].SetMaximum(h[detector+'leftvsright_signal_'+str(i)+xi].GetBinContent(10,10))
432 h[detector+'leftvsright_signal_'+str(i)+xi].Draw('colz')
433
434 ut.bookCanvas(h,detector+'LRinEff'+xi,' ',1800,450,3,1)
435 for s in range(1,4):
436 h[detector+'lLRinEff'+str(s)+xi]=ROOT.TLegend(0.6,0.54,0.99,0.93)
437 name = detector+'leftvsright_signal_'+str(s)+xi
438 h[name+'0Y'] = h[name].ProjectionY(name+'0Y',1,1)
439 h[name+'0X'] = h[name].ProjectionX(name+'0X',1,1)
440 h[name+'1X'] = h[name].ProjectionY(name+'1Y')
441 h[name+'1Y'] = h[name].ProjectionX(name+'1X')
442 tc = h[detector+'LRinEff'+xi].cd(s)
443 tc.SetLogy()
444 h[name+'0X'].SetStats(0)
445 h[name+'0Y'].SetStats(0)
446 h[name+'1X'].SetStats(0)
447 h[name+'1Y'].SetStats(0)
448 h[name+'0X'].SetLineColor(ROOT.kRed)
449 h[name+'0Y'].SetLineColor(ROOT.kGreen)
450 h[name+'1X'].SetLineColor(ROOT.kMagenta)
451 h[name+'1Y'].SetLineColor(ROOT.kCyan)
452 h[name+'0X'].SetMaximum(max(h[name+'1X'].GetMaximum(),h[name+'1Y'].GetMaximum()))
453 h[name+'0X'].Draw()
454 h[name+'0Y'].Draw('same')
455 h[name+'1X'].Draw('same')
456 h[name+'1Y'].Draw('same')
457 # Fill(Sleft,Sright)
458 h[detector+'lLRinEff'+str(s)+xi].AddEntry(h[name+'0X'],'left with no signal right',"f")
459 h[detector+'lLRinEff'+str(s)+xi].AddEntry(h[name+'0Y'],'right with no signal left',"f")
460 h[detector+'lLRinEff'+str(s)+xi].AddEntry(h[name+'1X'],'left all',"f")
461 h[detector+'lLRinEff'+str(s)+xi].AddEntry(h[name+'1Y'],'right all',"f")
462 h[detector+'lLRinEff'+str(s)+xi].Draw()
463
464 listSipmTypes = ['L','R','S']
465 plane_label = {}
466 # parameters for the signal median
467 q05=np.array([0.5])
468 signal_attributes={}
469#
470 for tag in ["","T"]:
471 ut.bookCanvas(h,tag+'signalUSVeto'+xi,' ',1200,1600,3,self.M.systemAndPlanes[1]+self.M.systemAndPlanes[2])
472 s = 1
473 l = 1
474 Xaxis_bin = 1
475 for plane in range(self.M.systemAndPlanes[1]):
476 for side in listSipmTypes:
477 tc = h[tag+'signalUSVeto'+xi].cd(l)
478 l+=1
479 if side=='S' or (plane==2 and side =='R'): continue
480 rc = h[detector+tag+'sig'+side+'_'+str( s*10+plane)+xi]
481 rc.Draw()
482 if side=='S' or tag=='T': continue
483 if plane==2: plane_label[Xaxis_bin] = "Veto "+str(plane)+" T"
484 else: plane_label[Xaxis_bin] = "Veto "+str(plane)+" "+side
485 med=np.array([0.])
486 rc.GetQuantiles(1,med,q05)
487 signal_attributes[Xaxis_bin]={
488 "median":med[0],
489 "std":rc.GetStdDev(),
490 "max":rc.FindLastBinAbove(0),
491 "percent_overflow":rc.GetBinContent(rc.GetNbinsX()+1)/rc.Integral()*100. if rc.Integral()>0 else 0.}
492 Xaxis_bin += 1
493 s=2
494 for plane in range(self.M.systemAndPlanes[2]):
495 for side in listSipmTypes:
496 tc = h[tag+'signalUSVeto'+xi].cd(l)
497 l+=1
498 rc = h[detector+tag+'sig'+side+'_'+str( s*10+plane)+xi]
499 rc.Draw()
500 if side=='S' or tag=='T': continue
501 plane_label[Xaxis_bin] = "US "+str(plane)+" "+side
502 med=np.array([0.])
503 rc.GetQuantiles(1,med,q05)
504 signal_attributes[Xaxis_bin]={
505 "median":med[0],
506 "std":rc.GetStdDev(),
507 "max":rc.FindLastBinAbove(0),
508 "percent_overflow":rc.GetBinContent(rc.GetNbinsX()+1)/rc.Integral()*100. if rc.Integral()>0 else 0.}
509 Xaxis_bin += 1
510 ut.bookCanvas(h,tag+'signalDS'+xi,' ',900,1600,2,self.M.systemAndPlanes[3])
511 s = 3
512 l = 1
513 for plane in range(self.M.systemAndPlanes[3]):
514 for side in listSipmTypes:
515 if side == 'S': continue
516 tc = h[tag+'signalDS'+xi].cd(l)
517 l+=1
518 if (plane%2==1 or plane==6) and side=='R': continue
519 rc = h[detector+tag+'sig'+side+'_'+str( s*10+plane)+xi]
520 rc.Draw()
521 if () or tag=='T': continue
522 if plane%2==1 or plane==6: plane_label[Xaxis_bin] = "DS "+str(plane//2+1)+" T"
523 else:plane_label[Xaxis_bin] = "DS "+str(plane//2+1)+" "+side
524 med=np.array([0.])
525 rc.GetQuantiles(1,med,q05)
526 signal_attributes[Xaxis_bin]={
527 "median":med[0],
528 "std":rc.GetStdDev(),
529 "max":rc.FindLastBinAbove(0),
530 "percent_overflow":rc.GetBinContent(rc.GetNbinsX()+1)/rc.Integral()*100. if rc.Integral()>0 else 0.}
531 Xaxis_bin += 1
532
533# summary canvases of the median, maximum and overflow of signals
534 ut.bookCanvas(h,detector+'signalsSummary'+xi,' ',1024,768,1,3)
535 signal_medians = [ROOT.TGraphErrors(), ROOT.TGraphErrors(), ROOT.TGraphErrors()]
536 signal_maxima = [ROOT.TGraphErrors(), ROOT.TGraphErrors(), ROOT.TGraphErrors()]
537 signal_overflow = [ROOT.TGraphErrors(), ROOT.TGraphErrors(), ROOT.TGraphErrors()]
538 Area = {}
539 point_count = {}
540 for s in self.M.systemAndPlanes:
541 point_count[s-1]=-1
542 if s ==1:
543 signal_medians[s-1].SetTitle("Median of signal per plane and per side")
544 signal_medians[s-1].GetYaxis().SetTitle("median QDC [a.u.]")
545 signal_maxima[s-1].SetTitle("Maximal signal per plane and per side")
546 signal_maxima[s-1].GetYaxis().SetTitle("maximum QDC [a.u.]")
547 signal_overflow[s-1].SetTitle("Overflow/All QDC per plane and per side")
548 signal_overflow[s-1].GetYaxis().SetTitle("overflow QDC [%]")
549 for item in signal_attributes.keys():
550 if s==2 and plane_label[item].find('US')<0 : continue
551 if s==3 and plane_label[item].find('DS')<0 : continue
552 point_count[s-1] += 1
553 signal_medians[s-1].SetPoint(point_count[s-1],item,signal_attributes[item]["median"])
554 signal_medians[s-1].SetPointError(point_count[s-1],0,signal_attributes[item]["std"])
555 signal_maxima[s-1].SetPoint(point_count[s-1],item,signal_attributes[item]["max"])
556 signal_maxima[s-1].SetPointError(point_count[s-1],0,0)
557 signal_overflow[s-1].SetPoint(point_count[s-1],item,signal_attributes[item]["percent_overflow"])
558 signal_overflow[s-1].SetPointError(point_count[s-1],0,0)
559 graph_list = [signal_medians, signal_maxima, signal_overflow]
560 for counter, graph in enumerate([item[0] for item in graph_list]):
561 h[detector+'signalsSummary'+xi].cd(counter+1)
562 ROOT.gPad.SetBottomMargin(0.2)
563 ROOT.gPad.SetGrid(0)
564 graph.Draw('AP')
565 xAxis = graph.GetXaxis()
566 # get rid of the original ticks
567 xAxis.SetTickLength(0)
568 ymin = graph.GetHistogram().GetMinimum()
569 ymax = graph.GetHistogram().GetMaximum()
570 for index, item in enumerate(signal_attributes.keys()):
571 bin_index = xAxis.FindBin(item)
572 xAxis.SetBinLabel(bin_index,plane_label[item])
573 # Draw custom grid lines for the X axis
574 grid = ROOT.TLine(graph.GetPointX(index), ymin, graph.GetPointX(index), ymax)
575 grid.SetLineStyle(3)
576 grid.DrawClone()
577 # Draw ticks
578 tick = ROOT.TLine(graph.GetPointX(index), ymin, graph.GetPointX(index), ymin + 0.03*(ymax-ymin))
579 tick.DrawClone()
580 # redraw the 1st veto graph and draw the graphs for the US and DS
581 graph.Draw('P,same')
582 graph_list[counter][1].Draw('P,same')
583 graph_list[counter][2].Draw('P,same')
584 for s in self.M.systemAndPlanes:
585 if s==1: Area[counter]={}
586 Area[counter][s-1] = ROOT.TBox(ROOT.TMath.MinElement(graph_list[counter][s-1].GetN(), graph_list[counter][s-1].GetX())-0.5,
587 graph.GetHistogram().GetMinimum(),
588 ROOT.TMath.MaxElement(graph_list[counter][s-1].GetN(), graph_list[counter][s-1].GetX())+0.5,
589 graph.GetHistogram().GetMaximum())
590 Area[counter][s-1].SetLineWidth(0)
591 Area[counter][s-1].SetFillStyle(3003)
592 Area[counter][s-1].SetFillColor(s+1)
593 Area[counter][s-1].Draw("same")
594 graph_list[counter][s-1].SetMarkerStyle(21)
595 graph_list[counter][s-1].SetMarkerColor(s+1)
596 graph_list[counter][s-1].SetLineColor(s+1)
597# end of the summary QDC shifter canvases
598
599 ut.bookCanvas(h,detector+"chanbar"+xi,' ',1800,700,3,1)
600 for s in self.M.systemAndPlanes:
601 opt = ""
602 if s==3:
603 ut.bookCanvas(h,sdict[s]+"chanbar"+xi,' ',1800,1800,12,15)
604 else:
605 y = self.M.systemAndPlanes[s]
606 ut.bookCanvas(h,sdict[s]+"chanbar"+xi,' ',1800,700,y,self.M.systemAndBars[s])
607 h[sdict[s]+"chanbar"+xi].cd(1)
608 for l in range(self.M.systemAndPlanes[s]):
609 if s==3 and (l==1 or l==3 or l==5 or l==6):continue
610 maxN = 0
611 for bar in range(self.M.systemAndBars[s]):
612 hname = detector+'chanmult_'+str(s*1000+100*l+bar)+xi
613 nmax = h[hname].GetBinContent(h[hname].GetMaximumBin())
614 if nmax > maxN : maxN = nmax
615 for bar in range(self.M.systemAndBars[s]):
616 hname = detector+'chanmult_'+str(s*1000+100*l+bar)+xi
617 h[hname].SetStats(0)
618 # h[hname].SetMaximum(maxN*1.2)
619 h[detector+"chanbar"+xi].cd(s)
620 h[hname].DrawClone(opt)
621 opt="same"
622 i = l+1 + (self.M.systemAndBars[s]-bar-1)*self.M.systemAndPlanes[s]
623 if s==3:
624 ix = bar//15 + 1 + (l//2)*4
625 iy = 15 - bar%15
626 i = ix+(iy-1)*12
627 h[sdict[s]+"chanbar"+xi].cd(i)
628 h[hname].SetMaximum(h[hname].GetBinContent(h[hname].GetMaximumBin())*1.2)
629 h[hname].Draw()
630
631# shifter summary plots
632 for item in ["Active", "Nfired"]:
633 ut.bookCanvas(h,detector+item+'ChannelsPerBarVeto',' ',1024,768,2,3)
634 counter = 1
635 s = 1
636 for l in range(self.M.systemAndPlanes[s]):
637 for i in range(2):
638 h[detector+item+'ChannelsPerBarVeto'].cd(counter)
639 counter += 1
640 if l==2 and i==1: continue # vertical planes have readout on the top only
641 if i%2==0:
642 h[detector+'chan'+item+'Left_'+str(s*10+l)].Draw("colz")
643 h[detector+'chan'+item+'Left_'+str(s*10+l)].SetMinimum(0)
644 else:
645 h[detector+'chan'+item+'Right_'+str(s*10+l)].Draw("colz")
646 h[detector+'chan'+item+'Right_'+str(s*10+l)].SetMinimum(0)
647 self.M.myPrint(h[detector+item+'ChannelsPerBarVeto'],detector+item+'ChannelsPerBarUSVeto',subdir='mufilter/shifter')
648
649 ut.bookCanvas(h,detector+item+'ChannelsPerBarUS',' ',1024,768,2,5)
650 counter = 1
651 s = 2
652 for l in range(self.M.systemAndPlanes[s]):
653 for i in range(2):
654 h[detector+item+'ChannelsPerBarUS'].cd(counter)
655 if i%2==0:
656 h[detector+'chan'+item+'Left_'+str(s*10+l)].Draw("colz")
657 h[detector+'chan'+item+'Left_'+str(s*10+l)].SetMinimum(0)
658 else:
659 h[detector+'chan'+item+'Right_'+str(s*10+l)].Draw("colz")
660 h[detector+'chan'+item+'Right_'+str(s*10+l)].SetMinimum(0)
661 counter += 1
662 self.M.myPrint(h[detector+item+'ChannelsPerBarUS'],detector+item+'ChannelsPerBarUSVeto',subdir='mufilter/shifter')
663
664 s=3
665 # the DS 'Active' channels plot has the different style than the Veto and US ones
666 if item=='Active':
667 ut.bookHist(h,detector+'chanActiveDSSummaryHisto', 'DS channel hit multiplicity; ;bar',
668 10, 0, 10,
669 self.M.systemAndBars[s],-0.5,self.M.systemAndBars[s]-0.5)
670 h[detector+'chanActiveDSSummaryHisto'].SetStats(0)
671 xAxis = h[detector+'chanActiveDSSummaryHisto'].GetXaxis()
672 ut.bookCanvas(h,detector+'chanActiveDSSummary',' ',1024,768,1,1)
673 counter = 0 # counting histo bins
674 for l in range(self.M.systemAndPlanes[s]):
675 for i in range(2):# sides left <-> right
676 if (l%2==1 or l==6) and i==1: continue# vertical planes have readout on the top only
677 counter += 1
678 xAxis.SetBinLabel(counter,plane_label[15+counter])
679 # loop over bars, each DS plane has 60 bars, significant histogram bins start from 1
680 for barIndex in range(1,self.M.systemAndBars[s]+1):
681 if i%2==0:
682 h[detector+'chanActiveDSSummaryHisto'].SetBinContent(counter, barIndex,
683 h[detector+'chan'+item+'Left_'+str(s*10+l)].GetBinContent(1, barIndex))
684 else:
685 h[detector+'chanActiveDSSummaryHisto'].SetBinContent(counter, barIndex,
686 h[detector+'chan'+item+'Right_'+str(s*10+l)].GetBinContent(1, barIndex))
687 h[detector+'chanActiveDSSummaryHisto'].SetMinimum(0)
688 h[detector+'chanActiveDSSummaryHisto'].Draw('colz')
689 h[detector+'chanActiveDSSummaryHisto'].SetStats(0)
690 self.M.myPrint(h[detector+'chanActiveDSSummary'],detector+'chanActiveDSSummary',subdir='mufilter/shifter')
691 # the 'Nfired' channels plot has the same style for all MuFi systems
692 else:
693 ut.bookCanvas(h,detector+item+'ChannelsPerBarDS',' ',1024,768,3,4)
694 counter = 1
695 for l in range(self.M.systemAndPlanes[s]):
696 for i in range(2):# sides left <-> right
697 h[detector+item+'ChannelsPerBarDS'].cd(counter)
698 if (l%2==1 or l==6) and i==1: continue# vertical planes have readout on the top only
699 if i%2==0:
700 h[detector+'chan'+item+'Left_'+str(s*10+l)].Draw("colz")
701 h[detector+'chan'+item+'Left_'+str(s*10+l)].SetMinimum(0)
702 else:
703 h[detector+'chan'+item+'Right_'+str(s*10+l)].Draw("colz")
704 h[detector+'chan'+item+'Right_'+str(s*10+l)].SetMinimum(0)
705 counter += 1
706 if l==5 and i==0: counter += 2
707 self.M.myPrint(h[detector+item+'ChannelsPerBarDS'],detector+item+'ChannelsPerBarDS',subdir='mufilter/shifter')
708#expert plots
709 canvas = detector+'signalsSummary'+xi
710 self.M.h[canvas].Update()
711 if xi!='': self.M.myPrint(self.M.h[canvas],canvas,subdir='mufilter/shifter/'+xi)
712 else: self.M.myPrint(self.M.h[canvas],canvas,subdir='mufilter/shifter')
713 for canvas in ['signalUSVeto'+xi,'signalDS'+xi,detector+'LR'+xi,'USBars'+xi,
714 "Vetochanbar"+xi,"USchanbar"+xi,"DSchanbar"+xi,'noise'+xi]:
715 h[canvas].Update()
716 if x!='': self.M.myPrint(h[canvas],canvas,subdir='mufilter/expert/'+xi)
717 else: self.M.myPrint(h[canvas],canvas,subdir='mufilter/expert')
718 for canvas in [detector+'hitmaps',detector+'Xhitmaps',detector+'barmaps',detector+'dTScifi',detector+'dTcorScifi']:
719 for s in range(1,4):
720 if s<3 and canvas.find('dT')>0: continue
721 h[canvas+sdict[s]+xi].Update()
722 if x!='': self.M.myPrint(h[canvas+sdict[s]+xi],canvas+sdict[s],subdir='mufilter/expert/'+xi)
723 else: self.M.myPrint(h[canvas+sdict[s]+xi],canvas+sdict[s],subdir='mufilter/expert')
724
725# tracking
726 ut.bookCanvas(h,"muonDSTracks"+xi,' ',1200,400,3,1)
727 tc = h["muonDSTracks"+xi].cd(1)
728 h[detector+'slopes'+xi].Draw('colz')
729 tc = h["muonDSTracks"+xi].cd(2)
730 rc = h[detector+'slopes'+xi].ProjectionX("slopeX"+xi)
731 rc.Draw()
732 rc.SetTitle('track Y slope')
733 tc = h["muonDSTracks"+xi].cd(3)
734 rc = h[detector+'slopes'+xi].ProjectionY("slopeY"+xi)
735 rc.Draw()
736 rc.SetTitle('track Y slope')
737
738 ut.bookCanvas(h,detector+'TtrackPos'+xi,"track position first state",600,1200,1,2)
739 h[detector+'TtrackPos'+xi].cd(1)
740 rc = h[detector+'trackPosBeam'+xi].Draw('colz')
741 h[detector+'TtrackPos'+xi].cd(2)
742 rc = h[detector+'trackPos'+xi].Draw('colz')
743 if x!='':
744 self.M.myPrint(h["muonDSTracks"+xi],"muonDSTrackdirection"+xi,subdir='mufilter/shifter/'+xi)
745 self.M.myPrint(self.M.h[detector+'TtrackPos'+xi],detector+'trackPos'+xi,subdir='mufilter/shifter/'+xi)
746 else:
747 self.M.myPrint(h["muonDSTracks"+xi],"muonDSTrackdirection"+xi,subdir='mufilter/shifter')
748 self.M.myPrint(self.M.h[detector+'TtrackPos'+xi],detector+'trackPos'+xi,subdir='mufilter/shifter')
749
750# residuals
751 # fit all mufi planes in the canvas depending on # veto planes
752 nh, nv = self.M.systemAndPlanes[1],15//self.M.systemAndPlanes[1]
753 ut.bookCanvas(h,detector+"residualsVsX"+xi,'residualsVsX ',1200,1200, nh, nv)
754 ut.bookCanvas(h,detector+"residualsVsY"+xi,'residualsVsY ',1200,1200, nh, nv)
755 ut.bookCanvas(h,detector+"residuals"+xi,'residuals',1200,1200, nh, nv)
756# veto 2 H planes veto 1 V plane US 5 H planes DS 3 H planes DS 4 V planes = 15 in total
757 for p in ['resX_','resY_']:
758 t = detector+"residualsVs"+p.replace('res','').replace('_','')+xi
759 i = 1
760 for s in range(1,4):
761 for l in range(7):
762 if s==1 and l>self.M.systemAndPlanes[1]-1: continue
763 if s==2 and l>self.M.systemAndPlanes[2]-1: continue
764 tc = h[t].cd(i)
765 hname = detector+p+sdict[s]+str(s*10+l)+xi
766 h[hname].Draw('colz')
767 if p.find('X')>0:
768 tc = h[detector+"residuals"+xi].cd(i)
769 h[hname+'proj']=h[hname].ProjectionX(hname+'proj')
770 rc = h[hname+'proj'].Fit('gaus','SQ')
771 fitResult = rc.Get()
772 if fitResult:
773 tc.Update()
774 stats = h[hname+'proj'].FindObject('stats')
775 stats.SetOptFit(1111111)
776 stats.SetX1NDC(0.68)
777 stats.SetY1NDC(0.31)
778 stats.SetX2NDC(0.98)
779 stats.SetY2NDC(0.94)
780 h[hname+'proj'].Draw()
781 i+=1
782 if x!='':
783 self.M.myPrint(self.M.h[detector+'residualsVsX'+xi],detector+'residualsVsX',subdir='mufilter/expert/'+xi)
784 self.M.myPrint(self.M.h[detector+'residualsVsY'+xi],detector+'residualsVsY',subdir='mufilter/expert/'+xi)
785 self.M.myPrint(self.M.h[detector+'residuals'+xi],detector+'residuals',subdir='mufilter/expert/'+xi)
786 else:
787 self.M.myPrint(self.M.h[detector+'residualsVsX'+xi],detector+'residualsVsX',subdir='mufilter/expert')
788 self.M.myPrint(self.M.h[detector+'residualsVsY'+xi],detector+'residualsVsY',subdir='mufilter/expert')
789 self.M.myPrint(self.M.h[detector+'residuals'+xi],detector+'residuals',subdir='mufilter/expert')
790
791 ut.bookCanvas(self.M.h,'dt'+xi,'',1200,1200,1,2)
792 self.M.h['dt'].cd(1)
793 self.M.h['deltaTScifiMufiHit_'+str(1)].Draw('hist')
794 for s in range(1,6):
795 self.M.h['deltaTScifiMufiHit_'+str(s)].SetStats(0)
796 self.M.h['deltaTScifiMufiHit_'+str(s)].SetLineColor(s+1)
797 self.M.h['deltaTScifiMufiHit_'+str(s)].Draw('samehist')
798 self.M.h['dt'].cd(2)
799 if 'B2noB1' in self.xing:
800 self.M.h['deltaTScifiMufiHit_'+str(1)+'B2noB1'].Draw('hist')
801 for s in range(1,6):
802 self.M.h['deltaTScifiMufiHit_'+str(s)+'B2noB1'].SetStats(0)
803 self.M.h['deltaTScifiMufiHit_'+str(s)+'B2noB1'].SetLineColor(s+1)
804 self.M.h['deltaTScifiMufiHit_'+str(s)+'B2noB1'].Draw('samehist')
805 self.M.myPrint(self.M.h['dt'+xi],'scifi DS hit difference',subdir='mufilter/expert')
806
807class Mufi_largeVSsmall(ROOT.FairTask):
808 """
809 make correlation plots of small and large sipms for US and Veto"
810 """
811 def Init(self,options,monitor):
812 self.M = monitor
813 sdict = self.M.sdict
814 h = self.M.h
815 run = ROOT.FairRunAna.Instance()
816 for S in [1,2]:
817 for l in range(monitor.systemAndPlanes[S]):
818 ut.bookHist(h,'SVSl_'+str(l),'QDC large vs small sum',200,0.,200.,200,0.,200.)
819 ut.bookHist(h,'sVSl_'+str(l),'QDC large vs small average',200,0.,200.,200,0.,200.)
820 for side in ['L','R']:
821 for i1 in range(monitor.systemAndPlanes[1]+monitor.systemAndPlanes[2]):
822 for i2 in range(i1+1,8):
823 tag=''
824 if S==2 and monitor.smallSiPMchannel(i1): tag = 's'+str(i1)
825 else: tag = 'l'+str(i1)
826 if S==2 and monitor.smallSiPMchannel(i2): tag += 's'+str(i2)
827 else: tag += 'l'+str(i2)
828 ut.bookHist(h,sdict[S]+'cor'+tag+'_'+side+str(l),'QDC channel i vs channel j',200,0.,200.,200,0.,200.)
829 for bar in range(monitor.systemAndBars[S]):
830 ut.bookHist(h,sdict[S]+'cor'+tag+'_'+side+str(l)+str(bar),'QDC channel i vs channel j',200,0.,200.,200,0.,200.)
831
832 def ExecuteEvent(self,event):
833 W = self.M.Weight
834 M = self.M
835 h = self.M.h
836 sdict = self.M.sdict
837 for aHit in event.Digi_MuFilterHits:
838 if not aHit.isValid(): continue
839 detID = aHit.GetDetectorID()
840 s = detID//10000
841 bar = (detID%1000)
842 if s>2: continue
843 l = (detID%10000)//1000 # plane number
844 sumL,sumS,SumL,SumS = 0,0,0,0
845 allChannels = M.map2Dict(aHit,'GetAllSignals',mask=False)
846 nS = 0
847 nL = 0
848 for c in allChannels:
849 if s==2 and self.M.smallSiPMchannel(c) :
850 sumS+= allChannels[c]
851 nS += 1
852 else:
853 sumL+= allChannels[c]
854 nL+=1
855 if nL>0: SumL=sumL/nL
856 if nS>0: SumS=sumS/nS
857 rc = h['sVSl_'+str(l)].Fill(SumS,SumL,W)
858 rc = h['SVSl_'+str(l)].Fill(sumS/4.,sumL/12.,W)
859#
860 for side in ['L','R']:
861 offset = 0
862 if side=='R': offset = 8
863 for i1 in range(offset,offset+7):
864 if not i1 in allChannels: continue
865 qdc1 = allChannels[i1]
866 for i2 in range(i1+1,offset+8):
867 if not (i2) in allChannels: continue
868 if s==2 and self.M.smallSiPMchannel(i1): tag = 's'+str(i1-offset)
869 else: tag = 'l'+str(i1-offset)
870 if s==2 and self.M.smallSiPMchannel(i2): tag += 's'+str(i2-offset)
871 else: tag += 'l'+str(i2-offset)
872 qdc2 = allChannels[i2]
873 rc = h[sdict[s]+'cor'+tag+'_'+side+str(l)].Fill(qdc1,qdc2,W)
874 rc = h[sdict[s]+'cor'+tag+'_'+side+str(l)+str(bar)].Fill(qdc1,qdc2,W)
875 allChannels.clear()
876
877 def Plot(self):
878 h = self.M.h
879 sdict = self.M.sdict
880 systemAndPlanes = self.M.systemAndPlanes
881 ut.bookCanvas(h,'TSL','',1800,1400,3,2)
882 ut.bookCanvas(h,'STSL','',1800,1400,3,2)
883 S=2
884 for l in range(systemAndPlanes[S]):
885 tc = h['TSL'].cd(l+1)
886 tc.SetLogz(1)
887 aHist = h['sVSl_'+str(l)]
888 aHist.SetTitle(';small SiPM QCD av:large SiPM QCD av')
889 nmax = aHist.GetBinContent(aHist.GetMaximumBin())
890 aHist.SetMaximum( 0.1*nmax )
891 tc = h['sVSl_'+str(l)].Draw('colz')
892 self.M.myPrint(h['TSL'],"largeSiPMvsSmallSiPM",subdir='mufilter/expert')
893 for l in range(systemAndPlanes[S]):
894 tc = h['STSL'].cd(l+1)
895 tc.SetLogz(1)
896 aHist = h['SVSl_'+str(l)]
897 aHist.SetTitle(';small SiPM QCD sum/2:large SiPM QCD sum/6')
898 nmax = aHist.GetBinContent(aHist.GetMaximumBin())
899 aHist.SetMaximum( 0.1*nmax )
900 tc = h['SVSl_'+str(l)].Draw('colz')
901 self.M.myPrint(h['STSL'],"SumlargeSiPMvsSmallSiPM",subdir='mufilter/expert')
902 for S in [1,2]:
903 for l in range(systemAndPlanes[S]):
904 for side in ['L','R']:
905 if S == 1 and l == 2 and side == 'R': continue
906 ut.bookCanvas(h,sdict[S]+'cor'+side+str(l),'',1800,1400,7,4)
907 k=1
908 for i1 in range(7):
909 for i2 in range(i1+1,8):
910 tag=''
911 if S==2 and self.M.smallSiPMchannel(i1): tag = 's'+str(i1)
912 else: tag = 'l'+str(i1)
913 if S==2 and self.M.smallSiPMchannel(i2): tag += 's'+str(i2)
914 else: tag += 'l'+str(i2)
915 tc = h[sdict[S]+'cor'+side+str(l)].cd(k)
916 for bar in range(self.M.systemAndBars[S]):
917 if bar == 0: h[sdict[S]+'cor'+tag+'_'+side+str(l)+str(bar)].Draw('colz')
918 else: h[sdict[S]+'cor'+tag+'_'+side+str(l)+str(bar)].Draw('colzsame')
919 k+=1
920 self.M.myPrint(h[sdict[S]+'cor'+side+str(l)],'QDCcor'+side+str(l),subdir='mufilter/expert')
921
922class Veto_Efficiency(ROOT.FairTask):
923 " calculate Veto efficiency against Scifi tracks "
924 def Init(self,options,monitor):
925 self.debug = True
926 self.deadTime = 100
927 self.M = monitor
928 sdict = self.M.sdict
929 self.eventBefore={'T':-1,'N':-1,'hits':{2:0,1:0,0:0,'0L':0,'0R':0,'1L':0,'1R':0}}
930 h = self.M.h
931 run = ROOT.FairRunAna.Instance()
932 self.trackTask = run.GetTask('simpleTracking')
933 if not self.trackTask: self.trackTask = run.GetTask('houghTransform')
934 ioman = ROOT.FairRootManager.Instance()
935 self.OT = ioman.GetSink().GetOutTree()
936 s = 1
937 self.noiseCuts = [1,5,10,12]
938 self.zEx = self.M.zPos['Scifi'][10]
939 for noiseCut in self.noiseCuts:
940 ut.bookHist(h,'timeDiffPrev_'+str(noiseCut),'time diff; [clock cycles] ',100,-0.5,999.5)
941 ut.bookHist(h,'XtimeDiffPrev_'+str(noiseCut),'time diff no hits; [clock cycles] ',100,-0.5,999.5)
942 ut.bookHist(h,'timeDiffNext_'+str(noiseCut),'time diff next; [clock cycles] ',100,-0.5,999.5)
943 ut.bookHist(h,'XtimeDiffNext_'+str(noiseCut),'time diff next no hits; [clock cycles] ',100,-0.5,999.5)
944 for c in ['','NoPrev','TiNoFi']:
945 for b in ['','beam']:
946 nc = 'T'+c+str(noiseCut)+b
947 for l in range(monitor.systemAndPlanes[s]):
948 ut.bookHist(h,nc+'PosVeto_'+str(l),'track pos at veto'+str(l)+' with hit '+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
949 ut.bookHist(h,nc+'XPosVeto_'+str(l),'track pos at veto'+str(l)+' no hit'+str(l)+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
950 ut.bookHist(h,nc+'XPosVetoXL_'+str(l),'track pos at veto'+str(l)+' no hit'+str(l)+';X [cm]; Y [cm]',1100,-55.,0.,1100,10.,65.)
951 ut.bookHist(h,nc+'PosVeto_111'+str(l),'track pos at veto AND hit'+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
952 ut.bookHist(h,nc+'XPosVeto_111'+str(l),'track pos at veto no hit'+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
953 if l == 0:
954 ut.bookHist(h,nc+'PosVeto_11','track pos at veto AND hit'+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
955 ut.bookHist(h,nc+'XPosVeto_11','track pos at veto no hit'+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
956 ut.bookHist(h,nc+'PosVeto_000','track pos at veto OR hit'+';X [cm]; Y [cm]',110,-55.,0.,110,10.,65.)
957 for x in h:
958 if isinstance(h[x], ROOT.TH2) and x.find("PosVeto")>0:
959 h[x].SetStats(0)
960
961 ut.bookHist(h,'hitVeto_0','nr hits L vs R;n sipm; n sipm',25,-0.5,24.5,25,-0.5,24.5)
962 ut.bookHist(h,'hitVeto_1','nr hits L vs R;n sipm; n sipm',25,-0.5,24.5,25,-0.5,24.5)
963 ut.bookHist(h,'hitVeto_2','nr hits T ;n sipm', 25,-0.5,24.5)
964 ut.bookHist(h,'hitVeto_01','nr hits 0 vs 1;n sipm; n sipm',25,-0.5,24.5,25,-0.5,24.5)
965 ut.bookHist(h,'hitVeto_02','nr hits 0 vs 2;n sipm; n sipm',25,-0.5,24.5,25,-0.5,24.5)
966 ut.bookHist(h,'hitVeto_12','nr hits 1 vs 2;n sipm; n sipm',25,-0.5,24.5,25,-0.5,24.5)
967 ut.bookHist(h,'scaler','all no prevEvent',25,-0.5,24.5)
968 ut.bookHist(h,'deltaT','delta T DS 2 and Scifi 1',100,-20.0,20.)
969 ut.bookHist(h,'X/Y','xy matching of scifi DS',100,-20.0,20.,100,-20.0,20.)
970
971 def ExecuteEvent(self,event):
972 scifiCorTest = False
973 systemAndPlanes = self.M.systemAndPlanes
974 sdict = self.M.sdict
975 s = 1
976 h = self.M.h
977 W = self.M.Weight
978 nSiPMs = self.M.MuFilter.GetConfParI("MuFilter/VetonSiPMs")
979 hits = {2:0,1:0,0:0,'0L':0,'0R':0,'1L':0,'1R':0}
980 vetoHitsFromPrev = 0
981 if event.EventHeader.GetRunId() < 6204 and event.EventHeader.GetRunId() > 5480: vetoHitsFromPrev = 5
982 # special treatment for first 10fb-1 in 2023, wrong time alignment, again!
983 N1 = event.GetReadEntry()
984 Tcurrent = event.EventHeader.GetEventTime()
985 dT = abs(Tcurrent-self.eventBefore['T'])
986 prevAdded = False
987 for j in [0,-1]:
988 if j<0 and N1>0:
989 if dT > vetoHitsFromPrev: continue
990 rc = event.GetEvent(N1-1) # add veto hits from prev event
991 prevAdded = True
992 for aHit in event.Digi_MuFilterHits:
993 if not aHit.isValid(): continue
994 Minfo = self.M.MuFilter_PlaneBars(aHit.GetDetectorID())
995 s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
996 if s>1: continue
997 allChannels = self.M.map2Dict(aHit,'GetAllSignals')
998 hits[l]+=len(allChannels)
999 if l != 2:
1000 for c in allChannels:
1001 if nSiPMs > c: # left side
1002 hits[str(l)+'L']+=1
1003 else:
1004 hits[str(l)+'R']+=1
1005 allChannels.clear()
1006 if prevAdded and N1>1:
1007 rc = event.GetEvent(N1-2)
1008 Tprevprev = event.EventHeader.GetEventTime()
1009 dT = abs(Tcurrent-Tprevprev)
1010 if prevAdded: event.GetEvent(N1)
1011 prevEvent = False
1012 tightNoiseFilter = None
1013 otherAdvTrigger = None
1014 otherFastTrigger = None
1015 if dT < self.deadTime and dT > vetoHitsFromPrev:
1016 prevEvent = True
1017 # check type of prev event, if it would pass tight noise filter, run 6568 ++
1018 if event.EventHeader.GetRunId() > 6567 and N1>0:
1019 tightNoiseFilter, otherFastTrigger, otherAdvTrigger,Nprev,dt = self.checkOtherTriggers(event)
1020
1021 tmpT = self.eventBefore['T']
1022 tmpN = self.eventBefore['N']
1023 self.eventBefore['T'] = Tcurrent
1024 if (self.M.EventNumber - self.eventBefore['N'] > 1) and self.M.options.postScale < 2:
1025 print('what is going on?', self.M.EventNumber, self.eventBefore['N'])
1026 self.eventBefore['N'] = self.M.EventNumber
1027
1028 rc = h['scaler'].Fill(11)
1029 if self.M.Reco_MuonTracks.GetEntries()<2: return # require Scifi and DS track
1030# check that track has scifi cluster in station 1, afterthought: require measurements in all planes
1031 planes = {}
1032 for scifiTrack in self.M.Reco_MuonTracks:
1033 if not scifiTrack.GetUniqueID()==1: continue
1034 fitStatus = scifiTrack.getFitStatus()
1035 if not fitStatus.isFitConverged(): continue
1036 if fitStatus.getNdf() < 5 or fitStatus.getNdf()>12 : continue
1037 if fitStatus.getChi2()/fitStatus.getNdf() > 80: continue
1038 for nM in range(scifiTrack.getNumPointsWithMeasurement()):
1039 M = scifiTrack.getPointWithMeasurement(nM)
1040 W = M.getRawMeasurement()
1041 detID = W.getDetId()
1042 planes[detID//100000] = 1
1043 rc = h['scaler'].Fill(10)
1044 scifiOneEff = 10 in planes or 11 in planes or not scifiCorTest
1045 if not scifiOneEff and len(planes) < 8: return
1046 if scifiOneEff and len(planes) < 10: return
1047 rc = h['scaler'].Fill(0)
1048 if not prevEvent: rc = h['scaler'].Fill(1)
1049
1050 for l in range(systemAndPlanes[1]-1):
1051 rc = h['hitVeto_'+str(l)].Fill(hits[str(l)+'L'],hits[str(l)+'R'])
1052 rc = h['hitVeto_01'].Fill(hits[0],hits[1])
1053 if systemAndPlanes[1] > 2 :
1054 rc = h['hitVeto_02'].Fill(hits[0],hits[2])
1055 rc = h['hitVeto_12'].Fill(hits[1],hits[2])
1056
1057 for muTrack in self.M.Reco_MuonTracks:
1058 if not muTrack.GetUniqueID()==3: continue
1059 fstate = muTrack.getFittedState()
1060 posT,momT = fstate.getPos(),fstate.getMom()
1061 lam = (self.zEx-posT.z())/momT.z()
1062 yExTag = posT.y()+lam*momT.y()
1063 xExTag = posT.x()+lam*momT.x()
1064 sstate = scifiTrack.getFittedState()
1065 pos = sstate.getPos()
1066 delX = xExTag - pos.x()
1067 delY = yExTag - pos.y()
1068 rc = h['X/Y'].Fill(delX,delY)
1069 if abs(delX)>10 or abs(delY)>10: return
1070
1071 for aTrack in self.M.Reco_MuonTracks:
1072 if not aTrack.GetUniqueID()==1: continue
1073 fitStatus = aTrack.getFitStatus()
1074 if not fitStatus.isFitConverged(): continue
1075 fstate = aTrack.getFittedState()
1076 pos,mom = [fstate.getPos(),fstate.getMom()]
1077 beam = False
1078 if abs(mom.x()/mom.z())<0.1 and abs(mom.y()/mom.z())<0.1: beam = True
1079# extrapolate to veto
1080# check timing, remove backward tracks
1081# 1: get early DS time in horizontal plane 2:
1082 dsHitTimes = []
1083 for aHit in event.Digi_MuFilterHits:
1084 if not aHit.isValid(): continue
1085 detID = aHit.GetDetectorID()
1086 Minfo = self.M.MuFilter_PlaneBars(detID)
1087 s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
1088 if s==3 and l==2:
1089 self.M.MuFilter.GetPosition(detID,A,B)
1090 timeLeft = aHit.GetTime(0)
1091 timeRight = aHit.GetTime(1)
1092 if timeLeft>0 and timeRight>0:
1093 dL = abs(A[0]-B[0])
1094 avTime = self.M.MuFilter.GetCorrectedTime(detID, 0, timeLeft*self.M.TDC2ns,0) + \
1095 self.M.MuFilter.GetCorrectedTime(detID, 1, timeRight*self.M.TDC2ns,0)
1096 dsHitTimes.append( (avTime-abs(A[0]-B[0])/15)/2) # L/2 / 15cm/ns
1097 dsHitTimes.sort()
1098 scifiHitTimes = {1:[],2:[],3:[],4:[],5:[]}
1099 deltaT = -100
1100 if len(dsHitTimes)>0:
1101 for scifiHit in event.Digi_ScifiHits:
1102 detID = scifiHit.GetDetectorID()
1103 s = int(scifiHit.GetDetectorID()/1000000)
1104 if s>1.5: continue
1105 scifiHitTimes[s].append(self.M.Scifi.GetCorrectedTime(detID,scifiHit.GetTime()*self.M.TDC2ns,0))
1106 for s in scifiHitTimes:
1107 if len(scifiHitTimes[s])<1: continue
1108 scifiHitTimes[s].sort()
1109 deltaT = dsHitTimes[0] - scifiHitTimes[s][0] - (self.M.zPos['MuFilter'][34]-self.M.zPos['Scifi'][s*10])/u.speedOfLight
1110 rc = h['deltaT'].Fill(deltaT)
1111 if deltaT < -10: continue
1112 #look for previous event time
1113 T1 = event.EventHeader.GetEventTime()
1114 N1 = event.EventHeader.GetEventNumber()
1115 if prevAdded: rc = event.GetEvent(N1-2)
1116 else: rc = event.GetEvent(N1-1)
1117 T0 = event.EventHeader.GetEventTime()
1118 rc = event.GetEvent(N1+1)
1119 T2 = event.EventHeader.GetEventTime()
1120 rc = event.GetEvent(N1)
1121 if (T1-T0) < 100 and self.M.options.postScale < 2:
1122 if not prevEvent: print('what is going on?',N1,T1,T0,N1-1,tmpN,tmpT)
1123 prevEvent = True
1124 s = 1
1125 xEx = {}
1126 yEx = {}
1127 for l in range(systemAndPlanes[1]):
1128 zEx = self.M.zPos['MuFilter'][s*10+l]
1129 lam = (zEx-pos.z())/mom.z()
1130 xEx[l] = pos.x()+lam*mom.x()
1131 yEx[l] = pos.y()+lam*mom.y()
1132 for l in range(systemAndPlanes[1]):
1133 for noiseCut in self.noiseCuts:
1134 c=''
1135 if not prevEvent: c='NoPrev'
1136 nc = 'T'+c+str(noiseCut)
1137 ncL = 'T'+'TiNoFi'+str(noiseCut)
1138 if hits[l] > noiseCut or (l==2 and hits[l] > int(noiseCut/2+0.5)):
1139 rc = h[nc+'PosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1140 if tightNoiseFilter: rc = h[ncL+'PosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1141 if beam: rc = h[nc+'beamPosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1142 else:
1143 rc = h[nc+'XPosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1144 rc = h[nc+'XPosVetoXL_'+str(l)].Fill(xEx[l],yEx[l])
1145 if tightNoiseFilter: h[ncL+'XPosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1146 if beam: rc = h[nc+'beamXPosVeto_'+str(l)].Fill(xEx[l],yEx[l])
1147 if l==0:
1148 if -45<xEx[l] and xEx[l]<-10 and 27<yEx[l] and yEx[l]<54:
1149 rc = h['timeDiffPrev_'+str(noiseCut)].Fill(T1-T0)
1150 rc = h['timeDiffNext_'+str(noiseCut)].Fill(T2-T1)
1151 if ( systemAndPlanes[1]==2 and hits[0] > noiseCut and hits[1] > noiseCut) or \
1152 ( systemAndPlanes[1]==3 and hits[0] > noiseCut and hits[1] > noiseCut \
1153 and hits[2] > int(noiseCut/2+0.5) ):
1154 for l1 in range(systemAndPlanes[1]):
1155 rc = h[nc+'PosVeto_111'+str(l1)].Fill(xEx[l1],yEx[l1])
1156 if tightNoiseFilter:
1157 rc = h[ncL+'PosVeto_111'+str(l1)].Fill(xEx[l1],yEx[l1])
1158 if beam: rc = h[nc+'beamPosVeto_1110'].Fill(xEx[l],yEx[l])
1159 if ( systemAndPlanes[1]==2 and (hits[0] > noiseCut or hits[1] > noiseCut)) or \
1160 ( systemAndPlanes[1]==3 and (hits[0] > noiseCut or hits[1] > noiseCut \
1161 or hits[2] > int(noiseCut/2+0.5)) ):
1162 rc = h[nc+'PosVeto_000'].Fill(xEx[l],yEx[l])
1163 if tightNoiseFilter: h[ncL+'PosVeto_000'].Fill(xEx[l],yEx[l])
1164 if beam: rc = h[nc+'beamPosVeto_000'].Fill(xEx[l],yEx[l])
1165 else:
1166 if -45<xEx[l] and xEx[l]<-10 and 27<yEx[l] and yEx[l]<54:
1167 rc = h['XtimeDiffPrev_'+str(noiseCut)].Fill(T1-T0)
1168 rc = h['XtimeDiffNext_'+str(noiseCut)].Fill(T2-T1)
1169 if not prevEvent or (prevEvent and not tightNoiseFilter):
1170 if self.debug: print('no hits',noiseCut,prevEvent,beam,N1,tightNoiseFilter,otherFastTrigger,otherAdvTrigger)
1171 for l1 in range(systemAndPlanes[1]):
1172 rc = h[nc+'XPosVeto_111'+str(l1)].Fill(xEx[l1],yEx[l1])
1173 if tightNoiseFilter:
1174 rc = h[ncL+'XPosVeto_111'+str(l1)].Fill(xEx[l1],yEx[l1])
1175 if beam: rc = h[nc+'beamXPosVeto_1110'].Fill(xEx[l],yEx[l])
1176 # also for the 2 planea
1177 if hits[0] > noiseCut and hits[1] > noiseCut:
1178 rc = h[nc+'PosVeto_11'].Fill(xEx[l],yEx[l])
1179 if not (hits[0] > noiseCut or hits[1] > noiseCut):
1180 rc = h[nc+'XPosVeto_11'].Fill(xEx[l],yEx[l])
1181
1182 def checkOtherTriggers(self,event,debug=False):
1183 T0 = event.EventHeader.GetEventTime()
1184 N = event.EventHeader.GetEventNumber()
1185 otherFastTrigger = False
1186 otherAdvTrigger = False
1187 tightNoiseFilter = False
1188 Nprev = 1
1189 if N<Nprev: return otherFastTrigger, otherAdvTrigger, tightNoiseFilter, -1, 0
1190 rc = event.GetEvent(N-Nprev)
1191 dt = T0 - event.EventHeader.GetEventTime()
1192 while dt < self.deadTime and N>Nprev:
1193 otherFastTrigger = False
1194 for x in event.EventHeader.GetFastNoiseFilters():
1195 if debug: print('fast:', x.first, x.second )
1196 if x.second and not x.first == 'Veto_Total': otherFastTrigger = True
1197 otherAdvTrigger = False
1198 for x in event.EventHeader.GetAdvNoiseFilters():
1199 if debug: print('adv:', x.first, x.second )
1200 if x.second and not x.first == 'VETO_Planes': otherAdvTrigger = True
1201 if debug: print('pre event ',Nprev,dt,otherFastTrigger,otherAdvTrigger)
1202 if otherFastTrigger and otherAdvTrigger:
1203 rc = event.GetEvent(N)
1204 return otherFastTrigger, otherAdvTrigger, tightNoiseFilter, Nprev, dt
1205 Nprev+=1
1206 rc = event.GetEvent(N-Nprev)
1207 dt = T0 - event.EventHeader.GetEventTime()
1208 Nprev = 1
1209 rc = event.GetEvent(N-Nprev)
1210 dt = T0 - event.EventHeader.GetEventTime()
1211 while dt < self.deadTime and Nprev>N:
1212 hits = {2:0,1:0,0:0}
1213 for aHit in event.Digi_MuFilterHits:
1214 Minfo = self.M.MuFilter_PlaneBars(aHit.GetDetectorID())
1215 s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
1216 if s>1: continue
1217 allChannels = aHit.GetAllSignals(False,False)
1218 hits[l]+=len(allChannels)
1219 noiseFilter0 = sum(hits)>4.5
1220 noiseFilter1 = all(p > 0 for p in hits)
1221 if debug: print('veto hits:',hits)
1222 if noiseFilter0 and noiseFilter1:
1223 tightNoiseFilter = True
1224 rc = event.GetEvent(N)
1225 return otherFastTrigger, otherAdvTrigger, tightNoiseFilter, Nprev-1, dt
1226 Nprev+=1
1227 rc = event.GetEvent(N-Nprev)
1228 dt = T0 - event.EventHeader.GetEventTime()
1229 if Nprev>1:
1230 rc = event.GetEvent(N-Nprev+1)
1231 dt = T0 - event.EventHeader.GetEventTime()
1232 rc = event.GetEvent(N)
1233 return otherFastTrigger, otherAdvTrigger, tightNoiseFilter, Nprev-1, dt
1234
1235 def Plot(self,beamOnly=False):
1236 h = self.M.h
1237 nVetoPlanes =self.M.systemAndPlanes[1]
1238 if nVetoPlanes > 2 : hist_idx = ['0','1','2','000','1110']
1239 else: hist_idx = ['0','1','000','1110']
1240 if beamOnly: b='beam'
1241 else: b=''
1242 for c in ['','NoPrev']:
1243 allTracks = h['T'+c+'1PosVeto_0'].Clone('tmp')
1244 allTracks.Add(h['T'+c+'1XPosVeto_0'])
1245 for noiseCut in self.noiseCuts:
1246 nc = 'T'+c+str(noiseCut)+b
1247 h[nc+'XPosVeto_000']=allTracks.Clone(nc+'XPosVeto_000')
1248 h[nc+'XPosVeto_000'].Add(h[nc+'PosVeto_000'],-1)
1249 for l in hist_idx:
1250 h[nc+'Veto_ineff'+l] = h[nc+'PosVeto_'+l].Clone(nc+'Veto_ineff'+l)
1251 if l != '2': h[nc+'Veto_ineff'+l].SetTitle('Veto inefficiency '+l+' noise cut='+str(noiseCut))
1252 else: h[nc+'Veto_ineff'+l].SetTitle('Veto inefficiency '+l+' noise cut='+str(int(noiseCut/2+0.5)))
1253 h[nc+'Veto_ineff'+l].SetMinimum(0)
1254 h[nc+'Veto_ineff'+l].SetMaximum(1)
1255 for ix in range(allTracks.GetNbinsX()):
1256 for iy in range(allTracks.GetNbinsY()):
1257 for l in hist_idx:
1258 bc = allTracks.GetBinContent(ix,iy)
1259 if bc < 100:
1260 h[nc+'Veto_ineff'+l].SetBinContent(ix,iy,-1)
1261 h[nc+'Veto_ineff'+l].SetBinError(ix,iy,0)
1262 else:
1263 h[nc+'Veto_ineff'+l].SetBinContent(ix,iy,max(h[nc+'XPosVeto_'+l].GetBinContent(ix+1,iy+1)/bc, 2.7/bc))
1264 h[nc+'Veto_ineff'+l].SetBinError(ix,iy,h[nc+'XPosVeto_'+l].GetBinError(ix+1,iy+1)/bc)
1265 ut.bookCanvas(h,nc+'VetoEff','',1400,1800, nVetoPlanes, 2*nVetoPlanes)
1266 for p in range(2*nVetoPlanes**2):
1267 tc = h[nc+'VetoEff'].cd(p+1)
1268 if p < nVetoPlanes:
1269 h[nc+'PosVeto_'+str(p%nVetoPlanes)].Draw('colz')
1270 if p in range(nVetoPlanes, 2*nVetoPlanes):
1271 h[nc+'PosVeto_111'+str(p%nVetoPlanes)].Draw('colz')
1272 if p in range(2*nVetoPlanes, 3*nVetoPlanes):
1273 h[nc+'XPosVeto_'+str(p%nVetoPlanes)].Draw('colz')
1274 if p in range(3*nVetoPlanes, 4*nVetoPlanes):
1275 h[nc+'XPosVeto_111'+str(p%nVetoPlanes)].Draw('colz')
1276 tc = h[nc+'VetoEff'].cd(4*nVetoPlanes+1)
1277 h[nc+'PosVeto_000'].Draw('colz')
1278
1279 ut.bookCanvas(h,nc+'VetoInEff','',1800,1400,nVetoPlanes,2)
1280 for p in range(nVetoPlanes):
1281 tc = h[nc+'VetoInEff'].cd(p+1)
1282 tc.SetLogz(1)
1283 h[nc+'Veto_ineff'+str(p)].Draw('colz')
1284 tc = h[nc+'VetoInEff'].cd(nVetoPlanes+1)
1285 tc.SetLogz(1)
1286 h[nc+'Veto_ineff1110'].Draw('colz')
1287 tc = h[nc+'VetoInEff'].cd(nVetoPlanes+2)
1288 tc.SetLogz(1)
1289 h[nc+'Veto_ineff000'].Draw('colz')
1290# make some printout
1291 Ntot = h[nc+'PosVeto_0'].Clone('Ntot')
1292 Ntot.Add(h[nc+'XPosVeto_0'])
1293 ineff0 = h[nc+'XPosVeto_0'].GetEntries()/(Ntot.GetEntries()+1E-20)
1294 ineff1 = h[nc+'XPosVeto_1'].GetEntries()/(Ntot.GetEntries()+1E-20)
1295 ineff2, ineffAND_old, ineffOR_old = 0, 0, 0
1296 if nVetoPlanes > 2:
1297 ineff2 = h[nc+'XPosVeto_2'].GetEntries()/(Ntot.GetEntries()+1E-20)
1298 ineffOR_old = h[nc+'XPosVeto_11'].GetEntries()/(Ntot.GetEntries()+1E-20)
1299 ineffAND_old = 1.-h[nc+'PosVeto_11'].GetEntries()/(Ntot.GetEntries()+1E-20)
1300 ineffOR = h[nc+'XPosVeto_1110'].GetEntries()/(Ntot.GetEntries()+1E-20)
1301 ineffAND = 1.-h[nc+'PosVeto_1110'].GetEntries()/(Ntot.GetEntries()+1E-20)
1302 region = [21,91,34,89]
1303 xax = h[nc+'PosVeto_0'].GetXaxis()
1304 yax = h[nc+'PosVeto_0'].GetYaxis()
1305 Ntot_r = Ntot.Integral(region[0],region[1],region[2],region[3])+1E-20
1306 ineff0_r = h[nc+'XPosVeto_0'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1307 ineff1_r = h[nc+'XPosVeto_1'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1308 ineff2_r, ineffAND_old_r, ineffOR_old_r = 0, 0, 0
1309 if nVetoPlanes > 2:
1310 ineff2_r = h[nc+'XPosVeto_2'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1311 ineffOR_r_old = h[nc+'XPosVeto_11'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1312 ineffAND_r_old = 1.-h[nc+'PosVeto_11'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1313 ineffOR_r = h[nc+'XPosVeto_1110'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1314 ineffAND_r = 1.-h[nc+'PosVeto_1110'].Integral(region[0],region[1],region[2],region[3])/Ntot_r
1315 print('noise cut = ',noiseCut, 'previous event:',c)
1316 print('global inefficiency veto0: %5.2F%% veto1: %5.2F%% veto2: %5.2F%%'
1317 %(ineff0*100,ineff1*100,ineff2*100),
1318 'vetoAND: %5.2F%% vetoOR: %5.2F%% veto0AND1: %5.2F%% veto0OR1: %5.2F%%'
1319 %(ineffAND*100,ineffOR*100,ineffAND_old*100,ineffOR_old*100))
1320 print('region %5.2F < X < %5.2F and %5.2F < Y < %5.2F '%(xax.GetBinCenter(region[0]),
1321 xax.GetBinCenter(region[1]),yax.GetBinCenter(region[0]),yax.GetBinCenter(region[1])))
1322 print('veto0: %5.2F%% veto1: %5.2F%% veto2: %5.2F%%'
1323 %(ineff0_r*100,ineff1_r*100,ineff2_r*100),
1324 'vetoAND: %5.2F%% vetoOR: %5.2F%% veto0AND1: %5.2F%% veto0OR1: %5.2F%%'
1325 %(ineffAND_r*100,ineffOR_r*100,ineffAND_old*100,ineffOR_old*100))
1326#
1327 for p in range(1,nVetoPlanes):
1328 if p == 1:
1329 h['hitVeto_X'] = h['hitVeto_0'+str(p)].ProjectionX('hitVeto_X')
1330 h['hitVeto_X'].SetStats(0)
1331 h['hitVeto_X'].SetLineColor(ROOT.kOrange)
1332 h['hitVeto_X'].SetLineWidth(3)
1333 h['hitVeto_Y'+str(p)] = h['hitVeto_0'+str(p)].ProjectionY('hitVeto_Y'+str(p))
1334 h['hitVeto_Y'+str(p)].SetLineColor(ROOT.kBlue-2*p)
1335 h['hitVeto_Y'+str(p)].SetStats(0)
1336 ut.bookCanvas(h,'ThitVeto','',900,600,1,1)
1337 tc = h['ThitVeto'].cd()
1338 tc.SetLogy(1)
1339 h['hitVeto_X'].Draw('hist')
1340 for p in range(1,nVetoPlanes): h['hitVeto_Y'+str(p)].Draw('histsame')
1341
1342 #save the Veto inefficiency plots to file
1343 if self.M.options.postScale<2:
1344 self.M.presenterFile.mkdir('mufilter/VetoIneff')
1345 for item in h:
1346 if isinstance(h[item], ROOT.TCanvas) and \
1347 (item.find('Eff')>0 or item.find('ThitVeto')>0):
1348 self.M.myPrint(h[item],item,subdir='mufilter/expert/VetoIneff')
MuFilter::GetPosition
void GetPosition(Int_t id, TVector3 &vLeft, TVector3 &vRight)
Definition MuFilter.cxx:639
MuFilter::GetCorrectedTime
Float_t GetCorrectedTime(Int_t id, Int_t c, Double_t t, Double_t L)
Definition MuFilter.cxx:579
MuFilter::GetConfParI
Int_t GetConfParI(TString name)
Definition MuFilter.h:47
Mufi_monitoring.Mufi_hitMaps
Definition Mufi_monitoring.py:17
Mufi_monitoring.Mufi_hitMaps.OT
OT
Definition Mufi_monitoring.py:34
Mufi_monitoring.Mufi_hitMaps.trackTask
trackTask
Definition Mufi_monitoring.py:31
Mufi_monitoring.Mufi_hitMaps.Init
Init(self, options, monitor)
Definition Mufi_monitoring.py:26
Mufi_monitoring.Mufi_hitMaps.Plot
Plot(self)
Definition Mufi_monitoring.py:338
Mufi_monitoring.Mufi_hitMaps.ExecuteEvent
ExecuteEvent(self, event)
Definition Mufi_monitoring.py:143
Mufi_monitoring.Mufi_hitMaps.xing
xing
Definition Mufi_monitoring.py:42
Mufi_monitoring.Mufi_hitMaps.M
M
Definition Mufi_monitoring.py:27
Mufi_monitoring.Mufi_hitMaps.mufi_vsignal
mufi_vsignal
Definition Mufi_monitoring.py:35
Mufi_monitoring.Mufi_hitMaps.listOfHits
listOfHits
Definition Mufi_monitoring.py:142
Mufi_monitoring.Mufi_hitMaps.beamSpot
beamSpot(self, event)
Definition Mufi_monitoring.py:317
Mufi_monitoring.Mufi_largeVSsmall
Definition Mufi_monitoring.py:807
Mufi_monitoring.Mufi_largeVSsmall.ExecuteEvent
ExecuteEvent(self, event)
Definition Mufi_monitoring.py:832
Mufi_monitoring.Mufi_largeVSsmall.Plot
Plot(self)
Definition Mufi_monitoring.py:877
Mufi_monitoring.Mufi_largeVSsmall.Init
Init(self, options, monitor)
Definition Mufi_monitoring.py:811
Mufi_monitoring.Mufi_largeVSsmall.M
M
Definition Mufi_monitoring.py:812
Mufi_monitoring.Veto_Efficiency
Definition Mufi_monitoring.py:922
Mufi_monitoring.Veto_Efficiency.checkOtherTriggers
checkOtherTriggers(self, event, debug=False)
Definition Mufi_monitoring.py:1182
Mufi_monitoring.Veto_Efficiency.noiseCuts
noiseCuts
Definition Mufi_monitoring.py:937
Mufi_monitoring.Veto_Efficiency.zEx
zEx
Definition Mufi_monitoring.py:938
Mufi_monitoring.Veto_Efficiency.trackTask
trackTask
Definition Mufi_monitoring.py:932
Mufi_monitoring.Veto_Efficiency.eventBefore
eventBefore
Definition Mufi_monitoring.py:929
Mufi_monitoring.Veto_Efficiency.OT
OT
Definition Mufi_monitoring.py:935
Mufi_monitoring.Veto_Efficiency.Init
Init(self, options, monitor)
Definition Mufi_monitoring.py:924
Mufi_monitoring.Veto_Efficiency.debug
debug
Definition Mufi_monitoring.py:925
Mufi_monitoring.Veto_Efficiency.M
M
Definition Mufi_monitoring.py:927
Mufi_monitoring.Veto_Efficiency.deadTime
deadTime
Definition Mufi_monitoring.py:926
Mufi_monitoring.Veto_Efficiency.ExecuteEvent
ExecuteEvent(self, event)
Definition Mufi_monitoring.py:971
Mufi_monitoring.Veto_Efficiency.Plot
Plot(self, beamOnly=False)
Definition Mufi_monitoring.py:1235
Scifi::GetConfParF
Float_t GetConfParF(TString name)
Definition Scifi.h:49
Scifi::GetCorrectedTime
Double_t GetCorrectedTime(Int_t fDetectorID, Double_t rawTime, Double_t L)
Definition Scifi.cxx:517