13 from modularAnalysis
import *
18 from ROOT
import Belle2
19 from ROOT
import gROOT, AddressOf
21 gROOT.ProcessLine(
'struct TrackData {\
24 float chiSquaredOverNdf;\
28 gROOT.ProcessLine(
'struct DEDXData {\
37 from ROOT
import TrackData, DEDXData
42 """A module to analyse residuals in overlaps of ladders."""
45 """Initialize the module"""
47 super(CosmicAnalysis, self).
__init__()
50 self.
rootfile = ROOT.TFile(
'cosmicAnalysis.root',
'recreate')
57 ROOT.gStyle.SetCanvasColor(0)
58 ROOT.gStyle.SetStatBorderSize(1)
59 ROOT.gStyle.SetStatColor(0)
60 ROOT.gStyle.SetTitleColor(1)
61 ROOT.gStyle.SetCanvasColor(0)
62 ROOT.gStyle.SetPadColor(0)
63 ROOT.gStyle.SetPadBorderMode(0)
64 ROOT.gStyle.SetCanvasBorderMode(0)
65 ROOT.gStyle.SetFrameBorderMode(0)
66 ROOT.gStyle.SetOptStat(0)
71 for key
in TrackData.__dict__:
74 if isinstance(self.
trackData.__getattribute__(key), int):
80 for key
in DEDXData.__dict__:
83 if isinstance(self.
dedxData.__getattribute__(key), int):
89 self.
HitsVsLayer = ROOT.TH2F(
'HitsVsLayer',
'', 6, 0.5, 6.5, 6, 0.5, 6.5)
91 self.
HitsVsSensor = ROOT.TH2F(
'HitsVsSensor',
'', 6, 0.5, 6.5, 5, 0.5, 5.5)
93 self.
LayerVsSensor = ROOT.TH2F(
'LayerVsSensor',
'', 6, 0.5, 6.5, 5, 0.5, 5.5)
97 self.
Chi2 = ROOT.TH1F(
'Chi2',
'', 300, 0.0, 500)
99 self.
NDF = ROOT.TH1F(
'NDF',
'', 200, 0.0, 200)
103 self.
Momentum = ROOT.TH1F(
'Momentum',
'', 500, 0.0, 1000)
108 'ADCCountOverNumberOfHitsInCDCVsMomentum',
'', 200, 0.0, 300, 0.0, 1000)
111 'MomentumVsADCCountOverNumberOfHitsInCDC',
'', 100, 0.0, 100, 0.0, 1000)
117 """ Fill histograms """
121 nCDCDedxTracks = cdcDedxTracks.getEntries()
122 for DEDXTrack_index
in range(nCDCDedxTracks):
123 dedxTrack = cdcDedxTracks[DEDXTrack_index]
124 self.dedxData.momentumCDC = dedxTrack.getMomentum()
125 self.dedxData.acdCount = 0
127 self.dedxData.path = 0
128 self.dedxData.numberOfHits = dedxTrack.size()
130 nKLMClusters = KLMClusters.getEntries()
131 if nKLMClusters == 2:
132 for nHit
in range(dedxTrack.size()):
133 self.dedxData.acdCount = self.dedxData.acdCount + dedxTrack.getADCCount(nHit)
134 self.dedxData.dE = self.dedxData.dE + dedxTrack.getDE(nHit)
135 self.dedxData.path = self.dedxData.path + dedxTrack.getPath(nHit)
136 self.dedxData.distanceKLM = math.sqrt(
138 KLMClusters[0].getClusterPosition().Z() -
139 KLMClusters[1].getClusterPosition().Z(),
142 self.ADCCountOverNumberOfHitsInCDC.Fill(self.dedxData.acdCount / dedxTrack.size())
143 self.MomentumVsADCCountOverNumberOfHitsInCDC.Fill(
144 self.dedxData.momentumCDC, self.dedxData.acdCount / dedxTrack.size())
145 self.ADCCountOverNumberOfHitsInCDCVsMomentum.Fill(
146 self.dedxData.acdCount / dedxTrack.size(), self.dedxData.momentumCDC)
149 self.tree_DEDX.Fill()
153 nRecoTracks = RecoTracks.getEntries()
155 for track_index
in range(nRecoTracks):
156 track = RecoTracks[track_index]
158 if track.wasFitSuccessful():
161 nTrackFitResults = TrackFitResults.getEntries()
162 if nTrackFitResults == 1:
163 self.trackData.momentum = TrackFitResults[0].getMomentum().Mag()
164 self.Momentum.Fill(TrackFitResults[0].getMomentum().Mag())
166 if track.hasPXDHits()
or track.hasSVDHits():
168 if track.getNumberOfSVDHits() % 2 == 0:
169 totalNumberOfHits = track.getNumberOfPXDHits() + track.getNumberOfSVDHits() / 2
171 totalNumberOfHits = track.getNumberOfPXDHits() + (track.getNumberOfSVDHits() - 1) / 2
173 self.TotalNumberOfHits.Fill(totalNumberOfHits)
176 self.Chi2.Fill(track.getTrackFitStatus().getChi2())
177 self.NDF.Fill(track.getTrackFitStatus().getNdf())
178 self.Chi2OverNDF.Fill(track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf())
180 self.trackData.chi2 = track.getTrackFitStatus().getChi2()
181 self.trackData.ndf = track.getTrackFitStatus().getNdf()
182 self.trackData.chiSquaredOverNdf = track.getTrackFitStatus().getChi2() / track.getTrackFitStatus().getNdf()
185 self.tree_track.Fill()
188 if track.hasPXDHits():
192 event = EventMetaData.getEvent()
193 nPxdHits = track.getNumberOfPXDHits()
197 for n
in range(0, len(track.getPXDHitList())):
198 pxdHit = track.getPXDHitList()[n]
200 info = geoCache.get(sensorID)
201 layer = sensorID.getLayerNumber()
202 ladder = sensorID.getLadderNumber()
203 sensor = sensorID.getSensorNumber()
205 self.HitsVsLayer.Fill(totalNumberOfHits, layer)
206 self.HitsVsSensor.Fill(totalNumberOfHits, sensor)
207 self.LayerVsSensor.Fill(layer, sensor)
208 self.PXDClusterSize.Fill(track.getPXDHitList()[n].getSize())
211 if track.hasSVDHits():
215 event = EventMetaData.getEvent()
216 nSvdHits = track.getNumberOfSVDHits()
220 measured = ROOT.TVector3(0, 0, 0)
221 expected = ROOT.TVector3(0, 0, 0)
222 for n
in range(0, len(track.getSVDHitList())):
223 svdHit = track.getSVDHitList()[n]
225 info = geoCache.get(sensorID)
226 layer = sensorID.getLayerNumber()
227 ladder = sensorID.getLadderNumber()
228 sensor = sensorID.getSensorNumber()
230 if svdHit.isUCluster():
232 for l
in range(0, len(track.getSVDHitList())):
233 svdHitA = track.getSVDHitList()[l]
235 if info == geoCache.get(sensorIDA):
236 if svdHitA.isUCluster() == 0:
238 self.HitsVsLayer.Fill(totalNumberOfHits, layer)
239 self.HitsVsSensor.Fill(totalNumberOfHits, sensor)
240 self.LayerVsSensor.Fill(layer, sensor)
244 print(
'Fit was not successful')
247 """Close & write output files"""
250 labelHits = (
'',
'One',
'Two',
'Three',
'Four',
'Five',
'Six')
251 labelLayers = (
'',
'First',
'Second',
'Third',
'Fourth',
'Fifth',
'Sixth')
252 labelSensors = (
'',
'First',
'Second',
'Third',
'Fourth',
'Fifth')
254 for i
in range(1, 7):
269 for i
in range(1, 7):
270 self.
HitsVsLayer.GetXaxis().SetBinLabel(i, labelHits[i])
271 self.
HitsVsLayer.GetYaxis().SetBinLabel(i, labelLayers[i])
274 self.
HitsVsLayer.GetXaxis().SetTitle(
"Number hits in VXD")
280 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
284 for i
in range(1, 7):
285 self.
HitsVsSensor.GetXaxis().SetBinLabel(i, labelHits[i])
286 for i
in range(1, 6):
287 self.
HitsVsSensor.GetYaxis().SetBinLabel(i, labelSensors[i])
290 self.
HitsVsSensor.GetXaxis().SetTitle(
"Number hits in VXD")
296 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
300 for i
in range(1, 7):
302 for i
in range(1, 6):
303 self.
LayerVsSensor.GetYaxis().SetBinLabel(i, labelSensors[i])
312 ROOT.gPad.SetMargin(0.11, 0.11, 0.1, 0.1)
316 self.
Chi2.GetXaxis().SetTitle(
"#Chi^{2}")
317 self.
Chi2.GetYaxis().SetTitle(
"Number of tracks")
318 self.
Chi2.GetYaxis().CenterTitle()
319 self.
Chi2.GetXaxis().CenterTitle()
320 self.
Chi2.GetYaxis().SetTitleOffset(1.3)
321 self.
Chi2.GetXaxis().SetTitleOffset(1.3)
322 self.
Chi2.SetFillStyle(3365)
323 self.
Chi2.SetFillColor(9)
324 self.
Chi2.SetLineColor(9)
328 self.
NDF.GetXaxis().SetTitle(
"Degrees of freedom")
329 self.
NDF.GetYaxis().SetTitle(
"Number of tracks")
330 self.
NDF.GetYaxis().CenterTitle()
331 self.
NDF.GetXaxis().CenterTitle()
332 self.
NDF.GetYaxis().SetTitleOffset(1.3)
333 self.
NDF.GetXaxis().SetTitleOffset(1.3)
334 self.
NDF.SetFillStyle(3365)
335 self.
NDF.SetFillColor(9)
336 self.
NDF.SetLineColor(9)
340 self.
Chi2OverNDF.SetTitle(
"#Chi^{2}/Degrees of freedom")
341 self.
Chi2OverNDF.GetYaxis().SetTitle(
"Number of tracks")
352 self.
Momentum.GetXaxis().SetTitle(
"Momentum [GeVc^{-1}]")
353 self.
Momentum.GetYaxis().SetTitle(
"Number of tracks")
354 self.
Momentum.GetYaxis().CenterTitle()
355 self.
Momentum.GetXaxis().CenterTitle()
356 self.
Momentum.GetYaxis().SetTitleOffset(1.3)
357 self.
Momentum.GetXaxis().SetTitleOffset(1.3)
376 function = ROOT.TF1(
"function",
"[0]+[1]*log(x)")
377 function.SetParameters(0.0, 0.0)
388 inverse_function = ROOT.TF1(
"inverse_function",
"exp((-[0]+x)/[1])", 0, 115)
389 inverse_function.SetParameter(0, function.GetParameter(0))
390 inverse_function.SetParameter(1, function.GetParameter(1))
391 constant = ROOT.TF1(
"constant",
"[0]", 115, 300)
392 constant.SetParameter(0, 20)
404 print(
'Parameters of fitted function y = exp((-p0+x)/p1) + p2, where x ~ ADCCounts/NumberOfCDCHits, y ~ momentum')
405 print(
'p0 =', inverse_function.GetParameter(0),
'p1 =', inverse_function.GetParameter(1),
'p2 =', constant.GetParameter(0))
423 main.add_module(
'RootInput')
424 main.add_module(
'Gearbox')
425 main.add_module(
'Geometry')
427 import reconstruction
as reco
428 reco.add_cosmics_reconstruction(main, pruneTracks=
False, merge_tracks=
True)
429 main.add_module(
'DAFRecoFitter', pdgCodesToUseForFitting=[13])
430 main.add_module(
"CDCDedxPID")
433 main.add_module(CosmicAnalysis)
435 progress = register_module(
'ProgressBar')
436 main.add_module(progress)