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