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