SND@LHC Software
Loading...
Searching...
No Matches
scifiHitMaps.py
Go to the documentation of this file.
1import ROOT,os
2import rootUtils as ut
3import shipunit as u
4h={}
5from argparse import ArgumentParser
6parser = ArgumentParser()
7parser.add_argument("-r", "--runNumber", dest="runNumber", help="run number", type=int,required=True)
8parser.add_argument("-p", "--path", dest="path", help="run number",required=False,default="")
9parser.add_argument("-f", "--inputFile", dest="fname", help="file name for MC", type=str,default=None,required=False)
10parser.add_argument("-g", "--geoFile", dest="geoFile", help="geofile", required=True)
11options = parser.parse_args()
12
13fgeo = ROOT.TFile.Open(options.path+options.geoFile)
14from ShipGeoConfig import ConfigRegistry
15from rootpyPickler import Unpickler
16#load geo dictionary
17upkl = Unpickler(fgeo)
18snd_geo = upkl.load('ShipGeo')
19# -----Create geometry----------------------------------------------
20import shipLHC_conf as sndDet_conf
21run = ROOT.FairRunSim()
22modules = sndDet_conf.configure(run,snd_geo)
23sGeo = fgeo.FAIRGeom
24modules['Scifi'].SiPMmapping()
25
26def xPos(detID):
27 orientation = (detID//100000)%10
28 nStation = 2*(detID//1000000-1)+orientation
29 mat = (detID%100000)//10000
30 X = detID%1000+(detID%10000)//1000*128
31 return [nStation,mat,X] # even numbers are Y (horizontal plane), odd numbers X (vertical plane)
32
33if options.runNumber>0:
34 f=ROOT.TFile.Open(options.path+'sndsw_raw_'+str(options.runNumber).zfill(6)+'.root')
35 eventTree = f.rawConv
36else:
37 f=ROOT.TFile.Open(options.fname)
38 eventTree = f.cbmsim
39
40def slopes(Nev=-1):
41 A,B = ROOT.TVector3(),ROOT.TVector3()
42 ut.bookHist(h,'slopesX','slope diffs',1000,-1.0,1.0)
43 ut.bookHist(h,'slopesY','slope diffs',1000,-1.0,1.0)
44 ut.bookHist(h,'clX','cluster size',10,0.5,10.5)
45 ut.bookHist(h,'clY','cluster size',10,0.5,10.5)
46# assuming cosmics make straight line
47 if Nev < 0 : Nev = eventTree.GetEntries()
48 for event in eventTree:
49 if Nev<0: break
50 Nev=Nev-1
51 clusters = []
52 hitDict = {}
53 for k in range(event.Digi_ScifiHits.GetEntries()):
54 d = event.Digi_ScifiHits[k]
55 if not d.isValid(): continue
56 hitDict[d.GetDetectorID()] = k
57 hitList = list(hitDict.keys())
58 if len(hitList)>0:
59 hitList.sort()
60 tmp = [ hitList[0] ]
61 cprev = hitList[0]
62 ncl = 0
63 last = len(hitList)-1
64 hitlist = ROOT.std.vector("sndScifiHit*")()
65 for i in range(len(hitList)):
66 if i==0 and len(hitList)>1: continue
67 c=hitList[i]
68 if (c-cprev)==1:
69 tmp.append(c)
70 if (c-cprev)!=1 or c==hitList[last]:
71 first = tmp[0]
72 N = len(tmp)
73 hitlist.clear()
74 for aHit in tmp: hitlist.push_back( event.Digi_ScifiHits[hitDict[aHit]])
75 aCluster = ROOT.sndCluster(first,N,hitlist,scifiDet)
76 clusters.append(aCluster)
77 if c!=hitList[last]:
78 ncl+=1
79 tmp = [c]
80 cprev = c
81 xHits = {}
82 yHits = {}
83 for s in range(5):
84 xHits[s]=[]
85 yHits[s]=[]
86 for aCl in clusters:
87 aCl.GetPosition(A,B)
88 vertical = int(aCl.GetFirst()/100000)%10==1
89 s = int(aCl.GetFirst()/1000000)-1
90 if vertical:
91 xHits[s].append(ROOT.TVector3(A))
92 rc = h['clX'].Fill(aCl.GetN())
93 else:
94 yHits[s].append(ROOT.TVector3(A))
95 rc = h['clY'].Fill(aCl.GetN())
96 proj = {'X':xHits,'Y':yHits}
97 for p in proj:
98 sls = []
99 for s1 in range(0,5):
100 if len(proj[p][s1]) !=1: continue
101 cl1 = proj[p][s1][0]
102 for s2 in range(s1+1,5):
103 if len(proj[p][s2]) !=1: continue
104 cl2 = proj[p][s2][0]
105 dz = abs(cl1[2]-cl2[2])
106 if dz < 5: continue
107 dzRep = 1./dz
108 m = dzRep*(cl2-cl1)
109 sls.append( m )
110 for ix1 in range(0,len(sls)-1):
111 for ix2 in range(ix1+1,len(sls)):
112 if p=="X": rc = h['slopes'+p].Fill( sls[ix2][0]-sls[ix1][0])
113 if p=="Y": rc = h['slopes'+p].Fill( sls[ix2][1]-sls[ix1][1])
114 ut.bookCanvas(h,'slopes',' ',1024,768,1,2)
115 h['slopes'].cd(1)
116 h['slopesX'].GetXaxis().SetRangeUser(-0.2,0.2)
117 h['slopesX'].SetTitle('x projection; delta slope [rad]')
118 h['slopesX'].Draw()
119 h['slopesX'].Fit('gaus','S','',-0.02,0.02)
120 h['slopes'].Update()
121 stats = h['slopesX'].FindObject('stats')
122 stats.SetOptFit(111)
123 h['slopes'].cd(2)
124 h['slopesY'].GetXaxis().SetRangeUser(-0.2,0.2)
125 h['slopesY'].SetTitle('y projection; delta slope [rad]')
126 h['slopesY'].Draw()
127 h['slopesY'].Fit('gaus','S','',-0.02,0.02)
128 h['slopes'].Update()
129 stats = h['slopesY'].FindObject('stats')
130 stats.SetOptFit(111)
131 for event in eventTree:
132 if Nev<0: break
133 Nev=Nev-1
134 clusters = []
135 hitDict = {}
136 for k in range(event.Digi_ScifiHits.GetEntries()):
137 d = event.Digi_ScifiHits[k]
138 if not d.isValid(): continue
139 hitDict[d.GetDetectorID()] = k
140 hitList = list(hitDict.keys())
141 if len(hitList)>0:
142 hitList.sort()
143 tmp = [ hitList[0] ]
144 cprev = hitList[0]
145 ncl = 0
146 last = len(hitList)-1
147 hitlist = ROOT.std.vector("sndScifiHit*")()
148 for i in range(len(hitList)):
149 if i==0 and len(hitList)>1: continue
150 c=hitList[i]
151 if (c-cprev)==1:
152 tmp.append(c)
153 if (c-cprev)!=1 or c==hitList[last]:
154 first = tmp[0]
155 N = len(tmp)
156 hitlist.clear()
157 for aHit in tmp: hitlist.push_back( event.Digi_ScifiHits[hitDict[aHit]])
158 aCluster = ROOT.sndCluster(first,N,hitlist,scifiDet)
159 clusters.append(aCluster)
160 if c!=hitList[last]:
161 ncl+=1
162 tmp = [c]
163 cprev = c
164 xHits = {}
165 yHits = {}
166 for s in range(5):
167 xHits[s]=[]
168 yHits[s]=[]
169 for aCl in clusters:
170 aCl.GetPosition(A,B)
171 vertical = int(aCl.GetFirst()/100000)%10==1
172 s = int(aCl.GetFirst()/1000000)-1
173 if vertical:
174 xHits[s].append(ROOT.TVector3(A))
175 rc = h['clX'].Fill(aCl.GetN())
176 else:
177 yHits[s].append(ROOT.TVector3(A))
178 rc = h['clY'].Fill(aCl.GetN())
179 proj = {'X':xHits,'Y':yHits}
180 for p in proj:
181 sls = []
182 for s1 in range(0,5):
183 if len(proj[p][s1]) !=1: continue
184 cl1 = proj[p][s1][0]
185 for s2 in range(s1+1,5):
186 if len(proj[p][s2]) !=1: continue
187
188def hitMaps(Nev=-1):
189 for mat in range(30):
190 ut.bookHist(h,'mat_'+str(mat),'hit map / mat',512,-0.5,511.5)
191 ut.bookHist(h,'sig_'+str(mat),'signal / mat',150,0.0,150.)
192 N=-1
193 if Nev < 0 : Nev = eventTree.GetEntries()
194 for event in eventTree:
195 N+=1
196 if N>Nev: break
197 for aHit in event.Digi_ScifiHits:
198 if not aHit.isValid(): continue
199 X = xPos(aHit.GetDetectorID())
200 rc = h['mat_'+str(X[0]*3+X[1])].Fill(X[2])
201 rc = h['sig_'+str(X[0]*3+X[1])].Fill(aHit.GetSignal(0))
202 ut.bookCanvas(h,'hitmaps',' ',1024,768,6,5)
203 ut.bookCanvas(h,'signal',' ',1024,768,6,5)
204 for mat in range(30):
205 tc = h['hitmaps'].cd(mat+1)
206 A = h['mat_'+str(mat)].GetSumOfWeights()/512.
207 if h['mat_'+str(mat)].GetMaximum()>10*A: h['mat_'+str(mat)].SetMaximum(10*A)
208 h['mat_'+str(mat)].Draw()
209 tc = h['signal'].cd(mat+1)
210 h['sig_'+str(mat)].Draw()
211def eventTime():
212 Tprev = -1
213 freq = 160.316E6
214 ut.bookHist(h,'Etime','delta event time; dt [s]',100,0.0,1.)
215 ut.bookCanvas(h,'E',' ',1024,2*768,1,2)
216 eventTree.GetEvent(0)
217 t0 = eventTree.EventHeader.GetEventTime()/160.E6
218 eventTree.GetEvent(eventTree.GetEntries()-1)
219 tmax = eventTree.EventHeader.GetEventTime()/160.E6
220 ut.bookHist(h,'time','elapsed time; t [s]',1000,t0,tmax)
221 for event in eventTree:
222 T = event.EventHeader.GetEventTime()
223 dT = 0
224 if Tprev >0: dT = T-Tprev
225 Tprev = T
226 rc = h['Etime'].Fill(dT/freq)
227 rc = h['time'].Fill( (T/freq-t0))
228 tc = h['E'].cd(1)
229 rc = h['Etime'].Fit('expo')
230 tc.Update()
231 stats = h['Etime'].FindObject("stats")
232 stats.SetOptFit(111)
233 tc = h['E'].cd(2)
234 h['time'].Draw()
235
236def mergeSignals(hstore):
237 ut.bookHist(hstore,'signalAll','signal all mat',150,0.0,150.)
238 for mat in range(30):
239 hstore['signalAll'].Add(hstore['sig_'+str(mat)])
240 hstore['signalAll'].Scale(1./hstore['signalAll'].GetSumOfWeights())
241
242def signalZoom(smax):
243 for mat in range(30):
244 h['sig_'+str(mat)].GetXaxis().SetRangeUser(0.,smax)
245 tc = h['signal'].cd(mat+1)
246 tc.Update()
247
248def beamSpot():
249 A,B = ROOT.TVector3(),ROOT.TVector3()
250 ut.bookHist(h,'bs','beam spot',100,-100.,10.,100,0.,80.)
251 for event in eventTree:
252 xMean = 0
253 yMean = 0
254 w=0
255 for d in event.Digi_ScifiHits:
256 detID = d.GetDetectorID()
257 s = int(detID/1000000)
258 modules['Scifi'].GetSiPMPosition(detID,A,B)
259 vertical = int(detID/100000)%10==1
260 if vertical: xMean+=A[0]
261 else: yMean+=A[1]
262 w+=1
263 rc = h['bs'].Fill(xMean/w,yMean/w)
scifiHitMaps.mergeSignals
mergeSignals(hstore)
Definition scifiHitMaps.py:236
scifiHitMaps.hitMaps
hitMaps(Nev=-1)
Definition scifiHitMaps.py:188
scifiHitMaps.slopes
slopes(Nev=-1)
Definition scifiHitMaps.py:40