22from ROOT
import Belle2
23from ROOT
import gROOT, AddressOf
24import reconstruction
as reco
26gROOT.ProcessLine(
'struct TrackData {\
29 float chiSquaredOverNdf;\
33gROOT.ProcessLine(
'struct DEDXData {\
42from ROOT
import TrackData, DEDXData
47 """A module to analyse residuals in overlaps of ladders."""
50 """Initialize the module"""
55 self.
rootfile = ROOT.TFile(
'cosmicAnalysis.root',
'recreate')
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)
76 for key
in TrackData.__dict__:
79 if isinstance(self.
trackData.__getattribute__(key), int):
85 for key
in DEDXData.__dict__:
88 if isinstance(self.
dedxData.__getattribute__(key), int):
94 self.
HitsVsLayer = ROOT.TH2F(
'HitsVsLayer',
'', 6, 0.5, 6.5, 6, 0.5, 6.5)
96 self.
HitsVsSensor = ROOT.TH2F(
'HitsVsSensor',
'', 6, 0.5, 6.5, 5, 0.5, 5.5)
98 self.
LayerVsSensor = ROOT.TH2F(
'LayerVsSensor',
'', 6, 0.5, 6.5, 5, 0.5, 5.5)
102 self.
Chi2 = ROOT.TH1F(
'Chi2',
'', 300, 0.0, 500)
104 self.
NDF = ROOT.TH1F(
'NDF',
'', 200, 0.0, 200)
108 self.
Momentum = ROOT.TH1F(
'Momentum',
'', 500, 0.0, 1000)
113 'ADCCountOverNumberOfHitsInCDCVsMomentum',
'', 200, 0.0, 300, 0.0, 1000)
116 'MomentumVsADCCountOverNumberOfHitsInCDC',
'', 100, 0.0, 100, 0.0, 1000)
122 """ Fill histograms """
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
132 self.dedxData.path = 0
133 self.dedxData.numberOfHits = dedxTrack.size()
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(
143 KLMClusters[0].getClusterPosition().Z() -
144 KLMClusters[1].getClusterPosition().Z(),
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)
154 self.tree_DEDX.Fill()
158 nRecoTracks = RecoTracks.getEntries()
160 for track_index
in range(nRecoTracks):
161 track = RecoTracks[track_index]
163 if track.wasFitSuccessful():
166 nTrackFitResults = TrackFitResults.getEntries()
167 if nTrackFitResults == 1:
168 self.trackData.momentum = TrackFitResults[0].getMomentum().Mag()
169 self.Momentum.Fill(TrackFitResults[0].getMomentum().Mag())
171 if track.hasPXDHits()
or track.hasSVDHits():
173 if track.getNumberOfSVDHits() % 2 == 0:
174 totalNumberOfHits = track.getNumberOfPXDHits() + track.getNumberOfSVDHits() / 2
176 totalNumberOfHits = track.getNumberOfPXDHits() + (track.getNumberOfSVDHits() - 1) / 2
178 self.TotalNumberOfHits.Fill(totalNumberOfHits)
181 self.Chi2.Fill(track.getTrackFitStatus().getChi2())
182 self.NDF.Fill(track.getTrackFitStatus().getNdf())
183 self.Chi2OverNDF.Fill(track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf())
185 self.trackData.chi2 = track.getTrackFitStatus().getChi2()
186 self.trackData.ndf = track.getTrackFitStatus().getNdf()
187 self.trackData.chiSquaredOverNdf = track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf()
190 self.tree_track.Fill()
193 if track.hasPXDHits():
202 for n
in range(0, len(track.getPXDHitList())):
203 pxdHit = track.getPXDHitList()[n]
205 info = geoCache.get(sensorID)
206 layer = sensorID.getLayerNumber()
207 sensor = sensorID.getSensorNumber()
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())
215 if track.hasSVDHits():
224 for n
in range(0, len(track.getSVDHitList())):
225 svdHit = track.getSVDHitList()[n]
227 info = geoCache.get(sensorID)
228 layer = sensorID.getLayerNumber()
229 sensor = sensorID.getSensorNumber()
231 if svdHit.isUCluster():
233 for lst
in range(0, len(track.getSVDHitList())):
234 svdHitA = track.getSVDHitList()[lst]
236 if info == geoCache.get(sensorIDA):
237 if svdHitA.isUCluster() == 0:
239 self.HitsVsLayer.Fill(totalNumberOfHits, layer)
240 self.HitsVsSensor.Fill(totalNumberOfHits, sensor)
241 self.LayerVsSensor.Fill(layer, sensor)
245 print(
'Fit was not successful')
248 """Close & write output files"""
251 labelHits = (
'',
'One',
'Two',
'Three',
'Four',
'Five',
'Six')
252 labelLayers = (
'',
'First',
'Second',
'Third',
'Fourth',
'Fifth',
'Sixth')
253 labelSensors = (
'',
'First',
'Second',
'Third',
'Fourth',
'Fifth')
255 for i
in range(1, 7):
270 for i
in range(1, 7):
271 self.
HitsVsLayer.GetXaxis().SetBinLabel(i, labelHits[i])
272 self.
HitsVsLayer.GetYaxis().SetBinLabel(i, labelLayers[i])
275 self.
HitsVsLayer.GetXaxis().SetTitle(
"Number hits in VXD")
281 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
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])
291 self.
HitsVsSensor.GetXaxis().SetTitle(
"Number hits in VXD")
297 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
301 for i
in range(1, 7):
303 for i
in range(1, 6):
304 self.
LayerVsSensor.GetYaxis().SetBinLabel(i, labelSensors[i])
313 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
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)
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)
341 self.
Chi2OverNDF.SetTitle(
"#Chi^{2}/Degrees of freedom")
342 self.
Chi2OverNDF.GetYaxis().SetTitle(
"Number of tracks")
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)
377 function = ROOT.TF1(
"function",
"[0]+[1]*log(x)")
378 function.SetParameters(0.0, 0.0)
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)
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))
423main = b2.create_path()
425main.add_module(
'RootInput')
426main.add_module(
'Gearbox')
427main.add_module(
'Geometry')
429reco.add_cosmics_reconstruction(main, pruneTracks=
False, merge_tracks=
True)
430main.add_module(
'DAFRecoFitter', pdgCodesToUseForFitting=[13])
431main.add_module(
"CDCDedxPID")
434main.add_module(CosmicAnalysis)
436progress = b2.register_module(
'ProgressBar')
437main.add_module(progress)
A (simplified) python wrapper for StoreArray.
static GeoCache & getInstance()
Return a reference to the singleton instance.
Class to uniquely identify a any structure of the PXD and SVD.
Momentum
Histogram with track momentum.
HitsVsSensor
Histogram with number of hits vs.
tree_track
Tree with track data.
TotalNumberOfHits
Histogram with total number of hits on track.
HitsVsLayer
Histogram with number of hits vs.
trackData
struct with track data
Chi2OverNDF
Histogram with track Chi2/NDF.
NDF
Histogram with track NDF.
Chi2
Histogram with track Chi2.
ADCCountOverNumberOfHitsInCDC
Histogram with ADC count vs.
LayerVsSensor
Histogram with hitted layer vs.
tree_DEDX
Tree with dE/dx data.
PXDClusterSize
Histogram with PXD cluster size.
ADCCountOverNumberOfHitsInCDCVsMomentum
Profile with ADC count vs.
dedxData
struct with dE/dx data
MomentumVsADCCountOverNumberOfHitsInCDC
Profile with ADC count vs.