Belle II Software development
testNewReco.py
1#!/usr/bin/env python3
2
3
10
11from simulation import add_simulation
12
13import basf2 as b2
14from svd import add_svd_reconstruction
15from ROOT import Belle2, TH1F, TH2F, TFile
16
17numEvents = 2000
18
19# bkgFiles = glob.glob('/sw/belle2/bkg/*.root') # Phase3 background
20bkgFiles = None # uncomment to remove background
21simulateJitter = False
22ROIfinding = False
23
24b2.set_random_seed(1234)
25# logging.log_level = LogLevel.DEBUG
26
27expList = [1003]
28
29
30class SVDClustersQuickCheck(b2.Module):
31 ''' quick check of cluster reconstruction'''
32
33 def initialize(self):
34 '''define histograms'''
35
36
37 self.test = []
38
39 self.testNew = []
40
41 self.size = TH1F("cl_size", "Cluster Size", 20, 0, 20)
42 self.sizeNew = TH1F("clNew_size", "New Cluster Size", 20, 0, 20)
43 self.time = TH1F("cl_time", "Cluster Time", 300, -100, 200)
44 self.timeNew = TH1F("clNew_time", "New Cluster Time", 300, -100, 200)
45 self.ff = TH1F("cl_ff", "Cluster FirstFrame", 4, -0.4, 3.5)
46 self.ffNew = TH1F("clNew_ff", "New Cluster FirstFrame", 4, -0.4, 3.5)
47 self.charge = TH1F("cl_charge", "Cluster Charge", 300, 0, 100)
48 self.chargeNew = TH1F("clNew_charge", "New Cluster Charge", 300, 0, 100)
49 self.SNR = TH1F("cl_SNR", "Cluster SNR", 100, 0, 100)
50 self.SNRNew = TH1F("clNew_SNR", "New Cluster SNR", 100, 0, 100)
51 self.position = TH1F("cl_position", "Cluster Position", 300, -6, 6)
52 self.positionNew = TH1F("clNew_position", "New Cluster Position", 300, -6, 6)
53 self.positionSigma = TH1F("cl_positionSigma", "Cluster Position Error", 300, 0, 100)
54 self.positionSigmaNew = TH1F("clNew_positionSigma", "New Cluster Position Error", 300, 0, 100)
55 self.positionS1 = TH1F("cl_positionS1", "Cluster Position Size 1", 300, -6, 6)
56 self.positionS1New = TH1F("clNew_positionS1", "New Cluster Position Size 1", 300, -6, 6)
57 self.positionS1Sigma = TH1F("cl_positionS1Sigma", "Cluster Position Error Size 1", 300, 0, 100)
58 self.positionS1SigmaNew = TH1F("clNew_positionS1Sigma", "New Cluster Position Error Size 1", 300, 0, 100)
59 self.positionS2 = TH1F("cl_positionS2", "Cluster Position Size 2", 300, -6, 6)
60 self.positionS2New = TH1F("clNew_positionS2", "New Cluster Position Size 2", 300, -6, 6)
61 self.positionS2Sigma = TH1F("cl_positionS2Sigma", "Cluster Position Size Error 2", 300, 0, 100)
62 self.positionS2SigmaNew = TH1F("clNew_positionS2Sigma", "New Cluster Position Error Size 2", 300, 0, 100)
63 self.positionS3 = TH1F("cl_positionS3", "Cluster Position Size >2", 300, -6, 6)
64 self.positionS3New = TH1F("clNew_positionS3", "New Cluster Position Size >2", 300, -6, 6)
65 self.positionS3Sigma = TH1F("cl_positionS3Sigma", "Cluster Position Error Size >2", 300, 0, 100)
66 self.positionS3SigmaNew = TH1F("clNew_positionS3Sigma", "New Cluster Position Error Size >2", 300, 0, 100)
67 self.positionPull = TH1F("cl_positionPull", "Cluster Position Pull", 200, -10, 10)
68 self.positionPullNew = TH1F("clNew_positionPull", "New Cluster Position Pull", 200, -10, 10)
69 self.positionPull1 = TH1F("cl_positionPull1", "Cluster Position Pull Size 1", 200, -10, 10)
70 self.positionPull1New = TH1F("clNew_positionPull1", "New Cluster Position Pull Size 1", 200, -10, 10)
71 self.positionPull2 = TH1F("cl_positionPull2", "Cluster Position Pull Size 2", 200, -10, 10)
72 self.positionPull2New = TH1F("clNew_positionPull2", "New Cluster Position Pull Size 2", 200, -10, 10)
73 self.positionPull3 = TH1F("cl_positionPull3", "Cluster Position Pull Size >2", 200, -10, 10)
74 self.positionPull3New = TH1F("clNew_positionPull3", "New Cluster Position Pull Size >2", 200, -10, 10)
75
76
77 geoCache = Belle2.VXD.GeoCache.getInstance()
78
79 ladderN = 0
80 sensorN = 0
81
82 for layer in geoCache.getLayers(Belle2.VXD.SensorInfoBase.SVD):
83 layerNumber = layer.getLayerNumber()
84 ladderN = 0
85 for ladder in geoCache.getLadders(layer):
86 ladderN = ladderN + 1
87 sensorN = 0
88 for sensor in geoCache.getSensors(ladder):
89 sensorN = sensorN + 1
90 self.test.append(TH2F("cl_layer" + str(layerNumber), "Layer " + str(layerNumber) +
91 " Ladder VS Sensor.Side", ladderN, 0.5, ladderN + 0.5, 2 * sensorN, +0.75, sensorN + 0.75))
92 self.testNew.append(TH2F("clNew_layer" + str(layerNumber), "Layer " + str(layerNumber) +
93 " Ladder VS Sensor.Side", ladderN, 0.5, ladderN + 0.5, 2 * sensorN, +0.75, sensorN + 0.75))
94
95 print(self.test)
96
97 def event(self):
98 ''' look at cluster and new clusters'''
99
100 clusterList = Belle2.PyStoreArray("SVDClusters")
101 clusterListNew = Belle2.PyStoreArray("SVDNewClusters")
102 print("number of clusters = " + str(clusterList.getEntries()) + " (old) VS " + str(clusterListNew.getEntries()) + " (new)")
103
104 for d in clusterList:
105 trueList = d.getRelationsTo('SVDTrueHits') # SVDClustersToSVDTrueHits
106 isU = 0
107 if (d.isUCluster()):
108 isU = 0.5
109 truePos = trueList[0].getU()
110 else:
111 truePos = trueList[0].getV()
112
113 self.size.Fill(d.getSize())
114 self.time.Fill(d.getClsTime())
115 self.ff.Fill(d.getFirstFrame())
116 self.charge.Fill(d.getCharge() / 1000)
117 self.SNR.Fill(d.getSNR())
118 self.position.Fill(d.getPosition())
119 self.positionSigma.Fill(d.getPositionSigma() * 1e4)
120 if d.getSize() == 1:
121 self.positionS1.Fill(d.getPosition())
122 self.positionS1Sigma.Fill(d.getPositionSigma() * 1e4)
123 self.positionPull1.Fill((d.getPosition() - truePos) / d.getPositionSigma())
124 if d.getSize() == 2:
125 self.positionS2.Fill(d.getPosition())
126 self.positionS2Sigma.Fill(d.getPositionSigma() * 1e4)
127 self.positionPull2.Fill((d.getPosition() - truePos) / d.getPositionSigma())
128 if d.getSize() > 2:
129 self.positionS3.Fill(d.getPosition())
130 self.positionS3Sigma.Fill(d.getPositionSigma() * 1e4)
131 self.positionPull3.Fill((d.getPosition() - truePos) / d.getPositionSigma())
132
133 self.test[d.getSensorID().getLayerNumber() - 3].Fill(d.getSensorID().getLadderNumber(),
134 d.getSensorID().getSensorNumber() + isU)
135 self.positionPull.Fill((d.getPosition() - truePos) / d.getPositionSigma())
136
137 for d in clusterListNew:
138
139 isU = 0
140 trueListNew = d.getRelationsTo('SVDTrueHits') # SVDClustersToSVDTrueHits
141 if (d.isUCluster()):
142 isU = 0.5
143 truePos = trueListNew[0].getU()
144 else:
145 truePos = trueListNew[0].getV()
146
147 self.sizeNew.Fill(d.getSize())
148 self.timeNew.Fill(d.getClsTime())
149 self.ffNew.Fill(d.getFirstFrame())
150 self.chargeNew.Fill(d.getCharge() / 1000)
151 self.SNRNew.Fill(d.getSNR())
152 self.positionNew.Fill(d.getPosition())
153 self.positionSigmaNew.Fill(d.getPositionSigma() * 1e4)
154 if d.getSize() == 1:
155 self.positionS1New.Fill(d.getPosition())
156 self.positionS1SigmaNew.Fill(d.getPositionSigma() * 1e4)
157 self.positionPull1New.Fill((d.getPosition() - truePos) / d.getPositionSigma())
158 if d.getSize() == 2:
159 self.positionS2New.Fill(d.getPosition())
160 self.positionS2SigmaNew.Fill(d.getPositionSigma() * 1e4)
161 self.positionPull2New.Fill((d.getPosition() - truePos) / d.getPositionSigma())
162 if d.getSize() > 2:
163 self.positionS3New.Fill(d.getPosition())
164 self.positionS3SigmaNew.Fill(d.getPositionSigma() * 1e4)
165 self.positionPull3New.Fill((d.getPosition() - truePos) / d.getPositionSigma())
166
167 self.testNew[d.getSensorID().getLayerNumber() - 3].Fill(d.getSensorID().getLadderNumber(),
168 d.getSensorID().getSensorNumber() + isU)
169 self.positionPullNew.Fill((d.getPosition() - truePos) / d.getPositionSigma())
170
171 def terminate(self):
172 '''write'''
173
174 f = TFile("quicktestSVDClusterOldDefault.root", "RECREATE")
175 for hist in self.test:
176 hist.GetXaxis().SetTitle("ladder #")
177 hist.GetYaxis().SetTitle("sensor # + 0.5 is isU")
178 hist.Write()
179 for hist in self.testNew:
180 hist.GetXaxis().SetTitle("ladder #")
181 hist.GetYaxis().SetTitle("sensor # + 0.5 is isU")
182 hist.Write()
183
184 self.size.GetXaxis().SetTitle("cluster size")
185 self.size.Write()
186 self.sizeNew.GetXaxis().SetTitle("cluster size")
187 self.sizeNew.Write()
188
189 self.time.GetXaxis().SetTitle("cluster time (ns)")
190 self.time.Write()
191 self.timeNew.GetXaxis().SetTitle("cluster time (ns)")
192 self.timeNew.Write()
193
194 self.ff.GetXaxis().SetTitle("cluster firstFrame")
195 self.ff.Write()
196 self.ffNew.GetXaxis().SetTitle("cluster firstFrame")
197 self.ffNew.Write()
198
199 self.charge.GetXaxis().SetTitle("cluster charge (ke-)")
200 self.charge.Write()
201 self.chargeNew.GetXaxis().SetTitle("cluster charge (ke-)")
202 self.chargeNew.Write()
203
204 self.SNR.GetXaxis().SetTitle("cluster SNR")
205 self.SNR.Write()
206 self.SNRNew.GetXaxis().SetTitle("cluster SNR")
207 self.SNRNew.Write()
208
209 self.position.GetXaxis().SetTitle("cluster position (cm)")
210 self.position.Write()
211 self.positionNew.GetXaxis().SetTitle("cluster position (cm)")
212 self.positionNew.Write()
213
214 self.positionSigma.GetXaxis().SetTitle("cluster position error (#mum)")
215 self.positionSigma.Write()
216 self.positionSigmaNew.GetXaxis().SetTitle("cluster position error (#mum)")
217 self.positionSigmaNew.Write()
218
219 self.positionS1.GetXaxis().SetTitle("cluster position (cm)")
220 self.positionS1.Write()
221 self.positionS1New.GetXaxis().SetTitle("cluster position (cm)")
222 self.positionS1New.Write()
223
224 self.positionS1Sigma.GetXaxis().SetTitle("cluster position error (#mum)")
225 self.positionS1Sigma.Write()
226 self.positionS1SigmaNew.GetXaxis().SetTitle("cluster position error (#mum)")
227 self.positionS1SigmaNew.Write()
228
229 self.positionS2.GetXaxis().SetTitle("cluster position (cm)")
230 self.positionS2.Write()
231 self.positionS2New.GetXaxis().SetTitle("cluster position (cm)")
232 self.positionS2New.Write()
233
234 self.positionS2Sigma.GetXaxis().SetTitle("cluster position error (#mum)")
235 self.positionS2Sigma.Write()
236 self.positionS2SigmaNew.GetXaxis().SetTitle("cluster position error (#mum)")
237 self.positionS2SigmaNew.Write()
238
239 self.positionS3.GetXaxis().SetTitle("cluster position (cm)")
240 self.positionS3.Write()
241 self.positionS3New.GetXaxis().SetTitle("cluster position (cm)")
242 self.positionS3New.Write()
243
244 self.positionS3Sigma.GetXaxis().SetTitle("cluster position error (#mum)")
245 self.positionS3Sigma.Write()
246 self.positionS3SigmaNew.GetXaxis().SetTitle("cluster position error (#mum)")
247 self.positionS3SigmaNew.Write()
248
249 self.positionPull.GetXaxis().SetTitle("cluster position pull")
250 self.positionPull.Write()
251 self.positionPullNew.GetXaxis().SetTitle("cluster position pull")
252 self.positionPullNew.Write()
253 self.positionPull1.GetXaxis().SetTitle("cluster position pull")
254 self.positionPull1.Write()
255 self.positionPull1New.GetXaxis().SetTitle("cluster position pull")
256 self.positionPull1New.Write()
257 self.positionPull2.GetXaxis().SetTitle("cluster position pull")
258 self.positionPull2.Write()
259 self.positionPull2New.GetXaxis().SetTitle("cluster position pull")
260 self.positionPull2New.Write()
261 self.positionPull3.GetXaxis().SetTitle("cluster position pull")
262 self.positionPull3.Write()
263 self.positionPull3New.GetXaxis().SetTitle("cluster position pull")
264 self.positionPull3New.Write()
265
266 f.Close()
267
268
269class SVDRecoDigitsQuickCheck(b2.Module):
270 '''quick check of SVDRecoDigits'''
271
272 def initialize(self):
273 '''define histograms'''
274
275
276 self.test = []
277
278 self.testNew = []
279
280 self.time = TH1F("rd_time", "RecoDigit Time", 300, -100, 200)
281 self.timeNew = TH1F("rdNew_time", "New RecoDigit Time", 300, -100, 200)
282 self.charge = TH1F("rd_charge", "RecoDigit Charge", 300, 0, 100000)
283 self.chargeNew = TH1F("rdNew_charge", "New RecoDigit Charge", 300, 0, 100000)
284
285
286 geoCache = Belle2.VXD.GeoCache.getInstance()
287
288 ladderN = 0
289 sensorN = 0
290
291 for layer in geoCache.getLayers(Belle2.VXD.SensorInfoBase.SVD):
292 layerNumber = layer.getLayerNumber()
293 ladderN = 0
294 for ladder in geoCache.getLadders(layer):
295 ladderN = ladderN + 1
296 sensorN = 0
297 for sensor in geoCache.getSensors(ladder):
298 sensorN = sensorN + 1
299 self.test.append(TH2F("rd_layer" + str(layerNumber), "Layer " + str(layerNumber) +
300 " Ladder VS Sensor.Side", ladderN, 0.5, ladderN + 0.5, 2 * sensorN, +0.75, sensorN + 0.75))
301 self.testNew.append(TH2F("rdNew_layer" + str(layerNumber), "Layer " + str(layerNumber) +
302 " Ladder VS Sensor.Side", ladderN, 0.5, ladderN + 0.5, 2 * sensorN, +0.75, sensorN + 0.75))
303
304 print(self.test)
305
306 def event(self):
307 '''look at old and new reco digits'''
308
309 recodigitList = Belle2.PyStoreArray("SVDRecoDigits")
310 recodigitListNew = Belle2.PyStoreArray("SVDNewRecoDigits")
311
312 print("number of recodigits = " + str(recodigitList.getEntries()) +
313 " (old) VS " + str(recodigitListNew.getEntries()) + " (new)")
314
315 for d in recodigitList:
316 self.time.Fill(d.getTime())
317 self.charge.Fill(d.getCharge())
318 isU = 0
319 if (d.isUStrip()):
320 isU = 0.5
321 self.test[d.getSensorID().getLayerNumber() - 3].Fill(d.getSensorID().getLadderNumber(),
322 d.getSensorID().getSensorNumber() + isU)
323
324 for d in recodigitListNew:
325 self.timeNew.Fill(d.getTime())
326 self.chargeNew.Fill(d.getCharge())
327 isU = 0
328 if (d.isUStrip()):
329 isU = 0.5
330 self.testNew[d.getSensorID().getLayerNumber() - 3].Fill(d.getSensorID().getLadderNumber(),
331 d.getSensorID().getSensorNumber() + isU)
332
333 def terminate(self):
334 '''write'''
335
336 f = TFile("quicktestSVDRecoDigitOldDefault.root", "RECREATE")
337 for hist in self.test:
338 hist.GetXaxis().SetTitle("ladder #")
339 hist.GetYaxis().SetTitle("sensor # + 0.5 is isU")
340 hist.Write()
341 for hist in self.testNew:
342 hist.GetXaxis().SetTitle("ladder #")
343 hist.GetYaxis().SetTitle("sensor # + 0.5 is isU")
344 hist.Write()
345
346 self.time.GetXaxis().SetTitle("recodigit time (ns)")
347 self.time.Write()
348 self.timeNew.GetXaxis().SetTitle("recodigit time (ns)")
349 self.timeNew.Write()
350
351 self.charge.GetXaxis().SetTitle("recodigit charge (ke-)")
352 self.charge.Write()
353 self.chargeNew.GetXaxis().SetTitle("recodigit charge (ke-)")
354 self.chargeNew.Write()
355
356 f.Close()
357
358# b2conditions.prepend_globaltag("svd_test_svdRecoConfiguration")
359
360
361main = b2.create_path()
362
363eventinfosetter = b2.register_module('EventInfoSetter')
364eventinfosetter.param('expList', expList)
365eventinfosetter.param('runList', [0])
366eventinfosetter.param('evtNumList', [numEvents])
367main.add_module(eventinfosetter)
368main.add_module('EventInfoPrinter')
369main.add_module('EvtGenInput')
370
371add_simulation(
372 main,
373 bkgfiles=bkgFiles,
374 usePXDDataReduction=ROIfinding,
375 simulateT0jitter=simulateJitter)
376
377add_svd_reconstruction(main)
378for mod in main.modules():
379 if (mod.name() == "SVDSimpleClusterizer"):
380 mod.param("timeAlgorithm", 0)
381 mod.param("HeadTailSize", 3)
382
383clusterizer = b2.register_module('SVDClusterizer')
384clusterizer.param('timeAlgorithm6Samples', "CoG6")
385clusterizer.param('timeAlgorithm3Samples', "CoG6")
386clusterizer.param('chargeAlgorithm6Samples', "MaxSample")
387clusterizer.param('chargeAlgorithm3Samples', "MaxSample")
388clusterizer.param('positionAlgorithm6Samples', "OldDefault")
389clusterizer.param('positionAlgorithm3Samples', "OldDefault")
390clusterizer.param('stripTimeAlgorithm6Samples', "dontdo")
391clusterizer.param('stripTimeAlgorithm3Samples', "dontdo")
392clusterizer.param('stripChargeAlgorithm6Samples', "MaxSample")
393clusterizer.param('stripChargeAlgorithm3Samples', "MaxSample")
394clusterizer.param('Clusters', "SVDNewClusters")
395clusterizer.param('useDB', False)
396main.add_module(clusterizer)
397
398recoDigitCreator = b2.register_module('SVDRecoDigitCreator')
399recoDigitCreator.param('timeAlgorithm6Samples', "CoG6")
400recoDigitCreator.param('timeAlgorithm3Samples', "CoG6")
401recoDigitCreator.param('chargeAlgorithm6Samples', "MaxSample")
402recoDigitCreator.param('chargeAlgorithm3Samples', "MaxSample")
403recoDigitCreator.param('RecoDigits', "SVDNewRecoDigits")
404recoDigitCreator.param('useDB', False)
405main.add_module(recoDigitCreator)
406
407main.add_module(SVDClustersQuickCheck())
408main.add_module(SVDRecoDigitsQuickCheck())
409
410# main.add_module('RootOutput')
411
412main.add_module('Progress')
413
414b2.print_path(main)
415
416b2.process(main, calculateStatistics=True)
417
418print(b2.statistics)
Belle2::PyStoreArray
A (simplified) python wrapper for StoreArray.
Definition PyStoreArray.h:73
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition GeoCache.cc:214
testNewReco.SVDClustersQuickCheck.test
list test
list of test histograms
Definition testNewReco.py:37
testNewReco.SVDClustersQuickCheck.testNew
list testNew
another list of test histograms
Definition testNewReco.py:39
testNewReco.SVDRecoDigitsQuickCheck.test
list test
list of test histograms
Definition testNewReco.py:276
testNewReco.SVDRecoDigitsQuickCheck.testNew
list testNew
another list of test histograms
Definition testNewReco.py:278