Belle II Software development
validatePhase2Cosmics.py
1#!/usr/bin/env python3
2
3
10
11# *****************************************************************************
12
13# title : 3_ValidationCosmics.py
14# description : Validation of cosmic tracks in phase II
15
16# *****************************************************************************
17
18import math
19import basf2 as b2
20
21import ROOT
22from ROOT import Belle2
23from ROOT import gROOT, AddressOf
24import reconstruction as reco
25
26gROOT.ProcessLine('struct TrackData {\
27 float chi2;\
28 float ndf;\
29 float chiSquaredOverNdf;\
30 float momentum;\
31};')
32
33gROOT.ProcessLine('struct DEDXData {\
34 float acdCount;\
35 float dE;\
36 float path;\
37 float momentumCDC;\
38 int numberOfHits;\
39 float distanceKLM;\
40};')
41
42from ROOT import TrackData, DEDXData # noqa
43
44
45class CosmicAnalysis(b2.Module):
46
47 """A module to analyse residuals in overlaps of ladders."""
48
49 def __init__(self):
50 """Initialize the module"""
51
52 super().__init__()
53
54
55 self.rootfile = ROOT.TFile('cosmicAnalysis.root', 'recreate')
56
57 self.tree_track = ROOT.TTree('track', '')
58
59 self.tree_DEDX = ROOT.TTree('dedx', '')
60
61 ROOT.gStyle.Reset()
62 ROOT.gStyle.SetCanvasColor(0)
63 ROOT.gStyle.SetStatBorderSize(1)
64 ROOT.gStyle.SetStatColor(0)
65 ROOT.gStyle.SetTitleColor(1)
66 ROOT.gStyle.SetCanvasColor(0)
67 ROOT.gStyle.SetPadColor(0)
68 ROOT.gStyle.SetPadBorderMode(0)
69 ROOT.gStyle.SetCanvasBorderMode(0)
70 ROOT.gStyle.SetFrameBorderMode(0)
71 ROOT.gStyle.SetOptStat(0)
72
73
74 self.trackData = TrackData()
75 # Declare tree branches
76 for key in TrackData.__dict__:
77 if '__' not in key:
78 formstring = '/F'
79 if isinstance(self.trackData.__getattribute__(key), int):
80 formstring = '/I'
81 self.tree_track.Branch(key, AddressOf(self.trackData, key), key + formstring)
82
83 self.dedxData = DEDXData()
84 # Declare tree branches
85 for key in DEDXData.__dict__:
86 if '__' not in key:
87 formstring = '/F'
88 if isinstance(self.dedxData.__getattribute__(key), int):
89 formstring = '/I'
90 self.tree_DEDX.Branch(key, AddressOf(self.dedxData, key), key + formstring)
91
92 self.TotalNumberOfHits = ROOT.TH1F('TotalNumberOfHits', '', 6, 0.5, 6.5)
93
94 self.HitsVsLayer = ROOT.TH2F('HitsVsLayer', '', 6, 0.5, 6.5, 6, 0.5, 6.5)
95
96 self.HitsVsSensor = ROOT.TH2F('HitsVsSensor', '', 6, 0.5, 6.5, 5, 0.5, 5.5)
97
98 self.LayerVsSensor = ROOT.TH2F('LayerVsSensor', '', 6, 0.5, 6.5, 5, 0.5, 5.5)
99
100 self.PXDClusterSize = ROOT.TH1F('PXDClusterSize', '', 20, 0.5, 20.5)
101
102 self.Chi2 = ROOT.TH1F('Chi2', '', 300, 0.0, 500)
103
104 self.NDF = ROOT.TH1F('NDF', '', 200, 0.0, 200)
105
106 self.Chi2OverNDF = ROOT.TH1F('Chi2OverNDF', '', 300, 0.0, 5)
107
108 self.Momentum = ROOT.TH1F('Momentum', '', 500, 0.0, 1000)
109
110 self.ADCCountOverNumberOfHitsInCDC = ROOT.TH1F('ADCCountOverNumberOfHitsInCDC', '', 200, 0.0, 300)
111
113 'ADCCountOverNumberOfHitsInCDCVsMomentum', '', 200, 0.0, 300, 0.0, 1000)
114
116 'MomentumVsADCCountOverNumberOfHitsInCDC', '', 100, 0.0, 100, 0.0, 1000)
117
118 def beginRun(self):
119 """Do nothing"""
120
121 def event(self):
122 """ Fill histograms """
123
124 # Study dEdx (CDC) as prediction of momentum
125 cdcDedxTracks = Belle2.PyStoreArray('CDCDedxTracks')
126 nCDCDedxTracks = cdcDedxTracks.getEntries()
127 for DEDXTrack_index in range(nCDCDedxTracks):
128 dedxTrack = cdcDedxTracks[DEDXTrack_index]
129 self.dedxData.momentumCDC = dedxTrack.getMomentum()
130 self.dedxData.acdCount = 0
131 self.dedxData.dE = 0
132 self.dedxData.path = 0
133 self.dedxData.numberOfHits = dedxTrack.size()
134 KLMClusters = Belle2.PyStoreArray('KLMClusters')
135 nKLMClusters = KLMClusters.getEntries()
136 if nKLMClusters == 2:
137 for nHit in range(dedxTrack.size()):
138 self.dedxData.acdCount = self.dedxData.acdCount + dedxTrack.getADCCount(nHit)
139 self.dedxData.dE = self.dedxData.dE + dedxTrack.getDE(nHit)
140 self.dedxData.path = self.dedxData.path + dedxTrack.getPath(nHit)
141 self.dedxData.distanceKLM = math.sqrt(
142 math.pow(
143 KLMClusters[0].getClusterPosition().Z() -
144 KLMClusters[1].getClusterPosition().Z(),
145 2))
146
147 self.ADCCountOverNumberOfHitsInCDC.Fill(self.dedxData.acdCount / dedxTrack.size())
148 self.MomentumVsADCCountOverNumberOfHitsInCDC.Fill(
149 self.dedxData.momentumCDC, self.dedxData.acdCount / dedxTrack.size())
150 self.ADCCountOverNumberOfHitsInCDCVsMomentum.Fill(
151 self.dedxData.acdCount / dedxTrack.size(), self.dedxData.momentumCDC)
152
153 self.rootfile.cd()
154 self.tree_DEDX.Fill()
155
156 # We are using RecoTracks to finding overlaps.
157 RecoTracks = Belle2.PyStoreArray('RecoTracks')
158 nRecoTracks = RecoTracks.getEntries()
160 for track_index in range(nRecoTracks):
161 track = RecoTracks[track_index]
162
163 if track.wasFitSuccessful():
164
165 TrackFitResults = Belle2.PyStoreArray('TrackFitResults')
166 nTrackFitResults = TrackFitResults.getEntries()
167 if nTrackFitResults == 1:
168 self.trackData.momentum = TrackFitResults[0].getMomentum().Mag()
169 self.Momentum.Fill(TrackFitResults[0].getMomentum().Mag())
170
171 if track.hasPXDHits() or track.hasSVDHits():
172
173 if track.getNumberOfSVDHits() % 2 == 0:
174 totalNumberOfHits = track.getNumberOfPXDHits() + track.getNumberOfSVDHits() / 2
175 else:
176 totalNumberOfHits = track.getNumberOfPXDHits() + (track.getNumberOfSVDHits() - 1) / 2
177
178 self.TotalNumberOfHits.Fill(totalNumberOfHits)
179 # print('Total number of hits:', totalNumberOfHits)
180
181 self.Chi2.Fill(track.getTrackFitStatus().getChi2())
182 self.NDF.Fill(track.getTrackFitStatus().getNdf())
183 self.Chi2OverNDF.Fill(track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf())
184
185 self.trackData.chi2 = track.getTrackFitStatus().getChi2()
186 self.trackData.ndf = track.getTrackFitStatus().getNdf()
187 self.trackData.chiSquaredOverNdf = track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf()
188 # print('Chi2/NDF:', self.trackData.chiSquaredOverNdf, 'Chi2:', self.trackData.chi2)
189 self.rootfile.cd()
190 self.tree_track.Fill()
191
192 # Check overlaps in PXD ladders & Check track in RecoTracks
193 if track.hasPXDHits():
194
195 # Print number of PXD hits
196 # EventMetaData = Belle2.PyStoreObj('EventMetaData')
197 # event = EventMetaData.getEvent()
198 # nPxdHits = track.getNumberOfPXDHits()
199 # print('Event', event, 'has PXD', nPxdHits, 'hit(s):')
200
201 # First loop over PXD Hits
202 for n in range(0, len(track.getPXDHitList())):
203 pxdHit = track.getPXDHitList()[n]
204 sensorID = Belle2.VxdID(pxdHit.getSensorID())
205 info = geoCache.get(sensorID)
206 layer = sensorID.getLayerNumber()
207 sensor = sensorID.getSensorNumber()
208 # print('Hit information: #hit:', n, 'sensor information:', layer, ladder, sensor)
209 self.HitsVsLayer.Fill(totalNumberOfHits, layer)
210 self.HitsVsSensor.Fill(totalNumberOfHits, sensor)
211 self.LayerVsSensor.Fill(layer, sensor)
212 self.PXDClusterSize.Fill(track.getPXDHitList()[n].getSize())
213
214 # Check overlaps in SVD ladders
215 if track.hasSVDHits():
216
217 # Print number of SVD hits
218 # EventMetaData = Belle2.PyStoreObj('EventMetaData')
219 # event = EventMetaData.getEvent()
220 # nSvdHits = track.getNumberOfSVDHits()
221 # print('Event', event, 'has SVD', nSvdHits, 'hit(s):')
222
223 # First loop over SVD Hits
224 for n in range(0, len(track.getSVDHitList())):
225 svdHit = track.getSVDHitList()[n]
226 sensorID = Belle2.VxdID(svdHit.getSensorID())
227 info = geoCache.get(sensorID)
228 layer = sensorID.getLayerNumber()
229 sensor = sensorID.getSensorNumber()
230
231 if svdHit.isUCluster():
232 # print('Hit information: #hit:', n, 'sensor information:', layer, ladder, sensor, 'isU')
233 for lst in range(0, len(track.getSVDHitList())):
234 svdHitA = track.getSVDHitList()[lst]
235 sensorIDA = Belle2.VxdID(svdHitA.getSensorID())
236 if info == geoCache.get(sensorIDA):
237 if svdHitA.isUCluster() == 0:
238 # print('Hit information: #hit:', n, 'sensor information:', layer, ladder, sensor)
239 self.HitsVsLayer.Fill(totalNumberOfHits, layer)
240 self.HitsVsSensor.Fill(totalNumberOfHits, sensor)
241 self.LayerVsSensor.Fill(layer, sensor)
242 # else:
243 # print('Hit information: #hit:', n, 'sensor information:', layer, ladder, sensor, 'isNotU')
244 else:
245 print('Fit was not successful')
246
247 def terminate(self):
248 """Close & write output files"""
249 self.rootfile.cd()
250
251 labelHits = ('', 'One', 'Two', 'Three', 'Four', 'Five', 'Six')
252 labelLayers = ('', 'First', 'Second', 'Third', 'Fourth', 'Fifth', 'Sixth')
253 labelSensors = ('', 'First', 'Second', 'Third', 'Fourth', 'Fifth')
254
255 for i in range(1, 7):
256 self.TotalNumberOfHits.GetXaxis().SetBinLabel(i, labelHits[i])
257 self.TotalNumberOfHits.GetXaxis().SetLabelSize(0.05)
258 self.TotalNumberOfHits.GetYaxis().SetTitle("Number of tracks in VXD")
259 self.TotalNumberOfHits.GetXaxis().SetTitle("Number of hits in VXD")
260 self.TotalNumberOfHits.GetYaxis().CenterTitle()
261 self.TotalNumberOfHits.GetXaxis().CenterTitle()
262 self.TotalNumberOfHits.GetYaxis().SetTitleOffset(1.3)
263 self.TotalNumberOfHits.GetXaxis().SetTitleOffset(1.3)
264 self.TotalNumberOfHits.SetFillStyle(3365)
265 self.TotalNumberOfHits.SetFillColor(9)
266 self.TotalNumberOfHits.SetLineColor(9)
267 self.TotalNumberOfHits.Draw()
268 # self.TotalNumberOfHits.Write()
269
270 for i in range(1, 7):
271 self.HitsVsLayer.GetXaxis().SetBinLabel(i, labelHits[i])
272 self.HitsVsLayer.GetYaxis().SetBinLabel(i, labelLayers[i])
273 self.HitsVsLayer.GetXaxis().SetLabelSize(0.05)
274 self.HitsVsLayer.GetYaxis().SetLabelSize(0.05)
275 self.HitsVsLayer.GetXaxis().SetTitle("Number hits in VXD")
276 self.HitsVsLayer.GetYaxis().SetTitle("Layer")
277 self.HitsVsLayer.GetYaxis().CenterTitle()
278 self.HitsVsLayer.GetXaxis().CenterTitle()
279 self.HitsVsLayer.GetYaxis().SetTitleOffset(1.6)
280 self.HitsVsLayer.GetXaxis().SetTitleOffset(1.3)
281 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
282 self.HitsVsLayer.Draw()
283 # self.HitsVsLayer.Write()
284
285 for i in range(1, 7):
286 self.HitsVsSensor.GetXaxis().SetBinLabel(i, labelHits[i])
287 for i in range(1, 6):
288 self.HitsVsSensor.GetYaxis().SetBinLabel(i, labelSensors[i])
289 self.HitsVsSensor.GetXaxis().SetLabelSize(0.05)
290 self.HitsVsSensor.GetYaxis().SetLabelSize(0.05)
291 self.HitsVsSensor.GetXaxis().SetTitle("Number hits in VXD")
292 self.HitsVsSensor.GetYaxis().SetTitle("Sensor")
293 self.HitsVsSensor.GetYaxis().CenterTitle()
294 self.HitsVsSensor.GetXaxis().CenterTitle()
295 self.HitsVsSensor.GetYaxis().SetTitleOffset(1.6)
296 self.HitsVsSensor.GetXaxis().SetTitleOffset(1.3)
297 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
298 self.HitsVsSensor.Draw()
299 # self.HitsVsSensor.Write()
300
301 for i in range(1, 7):
302 self.LayerVsSensor.GetXaxis().SetBinLabel(i, labelLayers[i])
303 for i in range(1, 6):
304 self.LayerVsSensor.GetYaxis().SetBinLabel(i, labelSensors[i])
305 self.LayerVsSensor.GetXaxis().SetLabelSize(0.05)
306 self.LayerVsSensor.GetYaxis().SetLabelSize(0.05)
307 self.LayerVsSensor.GetXaxis().SetTitle("Layer")
308 self.LayerVsSensor.GetYaxis().SetTitle("Sensor")
309 self.LayerVsSensor.GetYaxis().CenterTitle()
310 self.LayerVsSensor.GetXaxis().CenterTitle()
311 self.LayerVsSensor.GetYaxis().SetTitleOffset(1.6)
312 self.LayerVsSensor.GetXaxis().SetTitleOffset(1.3)
313 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
314 self.LayerVsSensor.Draw()
315 # self.LayerVsSensor.Write()
316
317 self.Chi2.GetXaxis().SetTitle("#Chi^{2}")
318 self.Chi2.GetYaxis().SetTitle("Number of tracks")
319 self.Chi2.GetYaxis().CenterTitle()
320 self.Chi2.GetXaxis().CenterTitle()
321 self.Chi2.GetYaxis().SetTitleOffset(1.3)
322 self.Chi2.GetXaxis().SetTitleOffset(1.3)
323 self.Chi2.SetFillStyle(3365)
324 self.Chi2.SetFillColor(9)
325 self.Chi2.SetLineColor(9)
326 self.Chi2.Draw()
327 # self.Chi2.Write()
328
329 self.NDF.GetXaxis().SetTitle("Degrees of freedom")
330 self.NDF.GetYaxis().SetTitle("Number of tracks")
331 self.NDF.GetYaxis().CenterTitle()
332 self.NDF.GetXaxis().CenterTitle()
333 self.NDF.GetYaxis().SetTitleOffset(1.3)
334 self.NDF.GetXaxis().SetTitleOffset(1.3)
335 self.NDF.SetFillStyle(3365)
336 self.NDF.SetFillColor(9)
337 self.NDF.SetLineColor(9)
338 self.NDF.Draw()
339 # self.NDF.Write()
340
341 self.Chi2OverNDF.SetTitle("#Chi^{2}/Degrees of freedom")
342 self.Chi2OverNDF.GetYaxis().SetTitle("Number of tracks")
343 self.Chi2OverNDF.GetYaxis().CenterTitle()
344 self.Chi2OverNDF.GetXaxis().CenterTitle()
345 self.Chi2OverNDF.GetYaxis().SetTitleOffset(1.3)
346 self.Chi2OverNDF.GetXaxis().SetTitleOffset(1.3)
347 self.Chi2OverNDF.SetFillStyle(3365)
348 self.Chi2OverNDF.SetFillColor(9)
349 self.Chi2OverNDF.SetLineColor(9)
350 self.Chi2OverNDF.Draw()
351 # self.Chi2OverNDF.Write()
352
353 self.Momentum.GetXaxis().SetTitle("Momentum [GeVc^{-1}]")
354 self.Momentum.GetYaxis().SetTitle("Number of tracks")
355 self.Momentum.GetYaxis().CenterTitle()
356 self.Momentum.GetXaxis().CenterTitle()
357 self.Momentum.GetYaxis().SetTitleOffset(1.3)
358 self.Momentum.GetXaxis().SetTitleOffset(1.3)
359 self.Momentum.SetFillStyle(3365)
360 self.Momentum.SetFillColor(9)
361 self.Momentum.SetLineColor(9)
362 self.Momentum.Draw()
363 # self.Momentum.Write()
364
365 self.ADCCountOverNumberOfHitsInCDC.GetXaxis().SetTitle("ADC Count [ADC] / number of hit")
366 self.ADCCountOverNumberOfHitsInCDC.GetYaxis().SetTitle("Number of tracks")
367 self.ADCCountOverNumberOfHitsInCDC.GetYaxis().CenterTitle()
368 self.ADCCountOverNumberOfHitsInCDC.GetXaxis().CenterTitle()
369 self.ADCCountOverNumberOfHitsInCDC.GetYaxis().SetTitleOffset(1.3)
370 self.ADCCountOverNumberOfHitsInCDC.GetXaxis().SetTitleOffset(1.3)
371 self.ADCCountOverNumberOfHitsInCDC.SetFillStyle(3365)
372 self.ADCCountOverNumberOfHitsInCDC.SetFillColor(9)
373 self.ADCCountOverNumberOfHitsInCDC.SetLineColor(9)
375 # self.ADCCountOverNumberOfHitsInCDC.Write()
376
377 function = ROOT.TF1("function", "[0]+[1]*log(x)")
378 function.SetParameters(0.0, 0.0)
379 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetYaxis().SetTitle("ADC Count [ADC] / number of hit")
380 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetXaxis().SetTitle("Momentum [GeVc^{-1}]")
381 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetYaxis().CenterTitle()
382 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetXaxis().CenterTitle()
383 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetYaxis().SetTitleOffset(1.3)
384 self.MomentumVsADCCountOverNumberOfHitsInCDC.GetXaxis().SetTitleOffset(1.3)
387 # self.MomentumVsADCCountOverNumberOfHitsInCDC.Write()
388
389 inverse_function = ROOT.TF1("inverse_function", "exp((-[0]+x)/[1])", 0, 115)
390 inverse_function.SetParameter(0, function.GetParameter(0))
391 inverse_function.SetParameter(1, function.GetParameter(1))
392 constant = ROOT.TF1("constant", "[0]", 115, 300)
393 constant.SetParameter(0, 20)
394 self.ADCCountOverNumberOfHitsInCDCVsMomentum.Fit(inverse_function, "R+")
395 self.ADCCountOverNumberOfHitsInCDCVsMomentum.Fit(constant, "R+")
396 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetXaxis().SetTitle("ADC Count [ADC] / number of hit")
397 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetYaxis().SetTitle("Momentum [GeVc^{-1}]")
398 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetYaxis().CenterTitle()
399 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetXaxis().CenterTitle()
400 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetYaxis().SetTitleOffset(1.3)
401 self.ADCCountOverNumberOfHitsInCDCVsMomentum.GetXaxis().SetTitleOffset(1.3)
403 # self.ADCCountOverNumberOfHitsInCDCVsMomentum.Write()
404
405 print('Parameters of fitted function y = exp((-p0+x)/p1) + p2, where x ~ ADCCounts/NumberOfCDCHits, y ~ momentum')
406 print('p0 =', inverse_function.GetParameter(0), 'p1 =', inverse_function.GetParameter(1), 'p2 =', constant.GetParameter(0))
407
408 self.PXDClusterSize.GetXaxis().SetTitle("Size of PXD Clusters")
409 self.PXDClusterSize.GetYaxis().SetTitle("Number of tracks")
410 self.PXDClusterSize.GetYaxis().CenterTitle()
411 self.PXDClusterSize.GetXaxis().CenterTitle()
412 self.PXDClusterSize.GetYaxis().SetTitleOffset(1.3)
413 self.PXDClusterSize.GetXaxis().SetTitleOffset(1.3)
414 self.PXDClusterSize.SetFillStyle(3365)
415 self.PXDClusterSize.SetFillColor(9)
416 self.PXDClusterSize.SetLineColor(9)
417 self.PXDClusterSize.Draw()
418
419 self.rootfile.Write()
420 self.rootfile.Close()
421
422
423main = b2.create_path()
424
425main.add_module('RootInput')
426main.add_module('Gearbox')
427main.add_module('Geometry')
428
429reco.add_cosmics_reconstruction(main, pruneTracks=False, merge_tracks=True)
430main.add_module('DAFRecoFitter', pdgCodesToUseForFitting=[13])
431main.add_module("CDCDedxPID")
432
433CosmicAnalysis = CosmicAnalysis()
434main.add_module(CosmicAnalysis)
435
436progress = b2.register_module('ProgressBar')
437main.add_module(progress)
438
439b2.process(main)
440
441# Print call statistics
442print(b2.statistics)
A (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:72
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:214
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
Momentum
Histogram with track momentum.
HitsVsSensor
Histogram with number of hits vs.
TotalNumberOfHits
Histogram with total number of hits on track.
HitsVsLayer
Histogram with number of hits vs.
Chi2OverNDF
Histogram with track Chi2/NDF.
ADCCountOverNumberOfHitsInCDC
Histogram with ADC count vs.
LayerVsSensor
Histogram with hitted layer vs.
PXDClusterSize
Histogram with PXD cluster size.
ADCCountOverNumberOfHitsInCDCVsMomentum
Profile with ADC count vs.
MomentumVsADCCountOverNumberOfHitsInCDC
Profile with ADC count vs.