Belle II Software development
extmuid_createPlots.py
1#!/usr/bin/env python3
2
3
10
11
17
18"""
19<header>
20 <input>muon-ExtMuidValidation.root,
21 pion-ExtMuidValidation.root</input>
22 <contact>piilonen@vt.edu</contact>
23 <description>Create validation plots for Ext and Muid</description>
24</header>
25"""
26
27import ROOT
28
29ROOT.PyConfig.IgnoreCommandLineOptions = True # noqa
30ROOT.gROOT.SetBatch(True) # noqa
31
32from ROOT import TChain, TFile, TH1F, TH2F, TNamed, gStyle
33import sys
34import math
35import numpy as np
36from optparse import OptionParser
37
38# contact person information
39# is added to the plot descriptions
40CONTACT_PERSON = {'Name': 'Leo Piilonen',
41 'Email': 'piilonen@vt.edu'}
42
43ACTIVE = True
44
45
46def main():
47 """Create validation plots for Ext and Muid"""
48
49 option_parser = OptionParser()
50 option_parser.add_option('-i', '--input-file', dest='input_file',
51 default='../*on-ExtMuidValidation.root',
52 help='Root file with Ext and Muid validation data.'
53 )
54 option_parser.add_option('-o', '--output-file', dest='output_file',
55 default='ExtMuidValidation.root',
56 help='Root file with Ext and Muid validation histograms.')
57
58 (options, args) = option_parser.parse_args()
59
60 # load chain of input files
61 file_chain = TChain('tree')
62 file_chain.Add(options.input_file)
63
64 # for file in file_chain:
65
66 number_entries = 0
67 try:
68 number_entries = file_chain.GetEntries()
69 except AttributeError:
70 print(f'Could not load input file(s) {options.input_file}.')
71
72 if number_entries == 0:
73 print(f'Data tree is empty or does not exist in file(s) {options.input_file}. Exit.')
74 sys.exit(0)
75
76 # open the output root file
77 output_root_file = TFile(options.output_file, 'recreate')
78
79 # create and draw histograms
80 gStyle.SetOptStat(0)
81 draw_exthits(file_chain)
82 draw_likelihoods(file_chain)
83
84 # close the output file
85 output_root_file.Write()
86 output_root_file.Close()
87
88
89def draw_exthits(file_chain):
90 """
91 Draw the ExtHit-related distributions.
92 """
93
94 contact_mail = CONTACT_PERSON["Email"]
95
96 # NOTE: *.Draw() must precede *.GetListOfFunctions().Add() or the latter will be discarded!
97 detectorID = TH1F('DetectorID', 'Detector ID for ExtHits', 8, -0.5, 7.5)
98 file_chain.Draw('ExtHits.m_DetectorID&0x0F>>DetectorID', '')
99 detectorID.GetXaxis().SetTitle('0=undefined, 1=PXD, 2=SVD, 3=CDC, 4=TOP, 5=ARICH, 6=ECL, 7=KLM')
100 detectorID.GetListOfFunctions().Add(TNamed('Description', "0=undefined, 1=PXD, 2=SVD, 3=CDC, 4=TOP, 5=ARICH, 6=ECL, 7=KLM"))
101 detectorID.GetListOfFunctions().Add(TNamed('Check', "ECL > (KLM ~ TOP) >> ARICH; others ~ 0"))
102 detectorID.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
103 detectorID.GetListOfFunctions().Add(TNamed('MetaOptions', 'shifter,pvalue-warn=0.50,pvalue-error=0.10'))
104 detectorID.Write()
105
106 tof = TH1F('TOF', 'Time of Flight for non-KLM ExtHits', 100, 0.0, 25.0)
107 file_chain.Draw('ExtHits.m_TOF>>TOF', '(ExtHits.m_DetectorID&0x0F)!=7')
108 tof.GetXaxis().SetTitle('t (ns)')
109 tof.GetListOfFunctions().Add(TNamed('Description', "Time of flight along track from IP to each ExtHit"))
110 tof.GetListOfFunctions().Add(TNamed('Check', "Dominant peak at ~4.5 ns, secondary peak at ~6 ns"))
111 tof.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
112 tof.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
113 tof.Write()
114
115 tofKLM = TH1F('TOFKLM', 'Time of Flight for KLM ExtHits', 100, 0.0, 25.0)
116 file_chain.Draw('ExtHits.m_TOF>>TOFKLM', '(ExtHits.m_DetectorID&0x0F)==7')
117 tofKLM.GetXaxis().SetTitle('t (ns)')
118 tofKLM.GetListOfFunctions().Add(TNamed('Description', "Time of flight along track from IP to each ExtHit"))
119 tofKLM.GetListOfFunctions().Add(TNamed('Check', "Broad rising distribution at 7-11 ns, then falling to ~17 ns"))
120 tofKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
121 tofKLM.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
122 tofKLM.Write()
123
124 r = TH1F('r', 'r for non-KLM ExtHits', 40, 0.0, 200.0)
125 file_chain.Draw('ExtHits.getPosition().Rho()>>r', '(ExtHits.m_DetectorID&0x0F)!=7')
126 r.GetXaxis().SetTitle('r (cm)')
127 r.GetListOfFunctions().Add(TNamed('Description', "Radial position in transverse plane of each ExtHit"))
128 r.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at 120 cm; broad peak at 140 cm"))
129 r.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
130 r.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
131 r.Write()
132
133 z = TH1F('z', 'z for non-KLM ExtHits', 100, -200.0, 300.0)
134 file_chain.Draw('ExtHits.getPosition().Z()>>z', '(ExtHits.m_DetectorID&0x0F)!=7')
135 z.GetXaxis().SetTitle('z (cm)')
136 z.GetListOfFunctions().Add(TNamed('Description', "Axial position of each ExtHit"))
137 z.GetListOfFunctions().Add(TNamed('Check', "Broad peak centered at 0 cm; sharp dip at 0; secondary peak at 170"))
138 z.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
139 z.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
140 z.Write()
141
142 rKLM = TH1F('rKLM', 'r for KLM ExtHits', 50, 100.0, 350.0)
143 file_chain.Draw('ExtHits.getPosition().Rho()>>rKLM', '(ExtHits.m_DetectorID&0x0F)==7')
144 rKLM.GetXaxis().SetTitle('r (cm)')
145 rKLM.GetListOfFunctions().Add(TNamed('Description', "Radial position in transverse plane of each ExtHit"))
146 rKLM.GetListOfFunctions().Add(TNamed('Check', "Low shoulder below 200 cm (EKLM); comb-like pattern above 200 cm (BKLM)"))
147 rKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
148 rKLM.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
149 rKLM.Write()
150
151 zKLM = TH1F('zKLM', 'z for KLM ExtHits', 140, -300.0, 400.0)
152 file_chain.Draw('ExtHits.getPosition().Z()>>zKLM', '(ExtHits.m_DetectorID&0x0F)==7')
153 zKLM.GetXaxis().SetTitle('z (cm)')
154 zKLM.GetListOfFunctions().Add(TNamed('Description', "Axial position of each ExtHit"))
155 zKLM.GetListOfFunctions().Add(TNamed('Check', "Comb-like pattern at z<-180 and z>280 cm (EKLM); broad peak at 0 (BKLM)"))
156 zKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
157 zKLM.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
158 zKLM.Write()
159
160 xy = TH2F('xy', 'y vs x for non-KLM ExtHits', 100, -200.0, 200.0, 100, -200.0, 200.0)
161 file_chain.Draw('ExtHits.getPosition().Y():ExtHits.getPosition().X()>>xy', '(ExtHits.m_DetectorID&0x0F)!=7')
162 xy.GetXaxis().SetTitle('x (cm)')
163 xy.GetYaxis().SetTitle('y (cm)')
164 xy.GetListOfFunctions().Add(TNamed('Description', "Position projected into transverse plane of each ExtHit"))
165 xy.GetListOfFunctions().Add(TNamed('Check', "Cylindrical, with most hits in ECL and TOP"))
166 xy.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
167 xy.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
168 xy.Write()
169
170 xz = TH2F('xz', 'x vs z for non-KLM ExtHits', 125, -200.0, 300.0, 100, -200.0, 200.0)
171 file_chain.Draw('ExtHits.getPosition().X():ExtHits.getPosition().Z()>>xz', '(ExtHits.m_DetectorID&0x0F)!=7')
172 xz.GetXaxis().SetTitle('z (cm)')
173 xz.GetYaxis().SetTitle('x (cm)')
174 xz.GetListOfFunctions().Add(TNamed('Description', "Position projected into x-z plane of each ExtHit"))
175 xz.GetListOfFunctions().Add(TNamed('Check', " "))
176 xz.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
177 xz.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
178 xz.Write()
179
180 yz = TH2F('yz', 'y vs z for non-KLM ExtHits', 125, -200.0, 300.0, 100, -200.0, 200.0)
181 file_chain.Draw('ExtHits.getPosition().Y():ExtHits.getPosition().Z()>>yz', '(ExtHits.m_DetectorID&0x0F)!=7')
182 yz.GetXaxis().SetTitle('z (cm)')
183 yz.GetYaxis().SetTitle('y (cm)')
184 yz.GetListOfFunctions().Add(TNamed('Description', "Position projected into y-z plane of each ExtHit"))
185 yz.GetListOfFunctions().Add(TNamed('Check', " "))
186 yz.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
187 yz.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
188 yz.Write()
189
190 xyKLM = TH2F('xyKLM', 'y vs x for KLM ExtHits', 140, -350.0, 350.0, 140, -350.0, 350.0)
191 file_chain.Draw('ExtHits.getPosition().Y():ExtHits.getPosition().X()>>xyKLM', '(ExtHits.m_DetectorID&0x0F)==7')
192 xyKLM.GetXaxis().SetTitle('x (cm)')
193 xyKLM.GetYaxis().SetTitle('y (cm)')
194 xyKLM.GetListOfFunctions().Add(TNamed('Description', "Position projected into transverse plane of each ExtHit"))
195 xyKLM.GetListOfFunctions().Add(TNamed('Check', "Octagonal (BKLM) and cylindrical (EKLM)"))
196 xyKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
197 xyKLM.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
198 xyKLM.Write()
199
200 xzKLM = TH2F('xzKLM', 'x vs z for KLM ExtHits', 140, -300.0, 400.0, 140, -350.0, 350.0)
201 file_chain.Draw('ExtHits.getPosition().X():ExtHits.getPosition().Z()>>xzKLM', '(ExtHits.m_DetectorID&0x0F)==7')
202 xzKLM.GetXaxis().SetTitle('z (cm)')
203 xzKLM.GetYaxis().SetTitle('x (cm)')
204 xzKLM.GetListOfFunctions().Add(TNamed('Description', "Position projected into x-z plane of each ExtHit"))
205 xzKLM.GetListOfFunctions().Add(TNamed('Check', " "))
206 xzKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
207 xzKLM.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
208 xzKLM.Write()
209
210 yzKLM = TH2F('yzKLM', 'y vs z for KLM ExtHits', 140, -300.0, 400.0, 140, -350.0, 350.0)
211 file_chain.Draw('ExtHits.getPosition().Y():ExtHits.getPosition().Z()>>yzKLM', '(ExtHits.m_DetectorID&0x0F)==7')
212 yzKLM.GetXaxis().SetTitle('z (cm)')
213 yzKLM.GetYaxis().SetTitle('y (cm)')
214 yzKLM.GetListOfFunctions().Add(TNamed('Description', "Position projected into y-z plane of each ExtHit"))
215 yzKLM.GetListOfFunctions().Add(TNamed('Check', " "))
216 yzKLM.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
217 yzKLM.GetListOfFunctions().Add(TNamed('MetaOptions', "box, expert"))
218 yzKLM.Write()
219
220
221def draw_likelihoods(file_chain):
222 """
223 Draw the Muid likelihood-based distributions.
224 """
225
226 contact_mail = CONTACT_PERSON["Email"]
227
228 outcome = TH1F('Outcome', 'Outcome', 5, -0.5, 4.5)
229
230 llBins = 120
231 llMin = -60.0
232 llMax = 0.0
233 llMu_mu = TH1F('LLmu-muons', 'logL(mu) for true muons', llBins, llMin, llMax)
234 llMu_pi = TH1F('LLmu-pions', 'logL(mu) for true pions', llBins, llMin, llMax)
235 llPi_mu = TH1F('LLpi-muons', 'logL(pi) for true muons', llBins, llMin, llMax)
236 llPi_pi = TH1F('LLpi-pions', 'logL(pi) for true pions', llBins, llMin, llMax)
237 llDiff_mu = TH1F('LLdiff-muons', 'logL(mu)-logL(pi) for true muons', llBins * 2, llMin, -llMin)
238 llDiff_pi = TH1F('LLdiff-pions', 'logL(mu)-logL(pi) for true pions', llBins * 2, llMin, -llMin)
239
240 layerBins = 15
241 layerMin = -0.5
242 layerMax = 14.5
243 barrelLayer_mu = TH1F('BarrelLayer-muons', 'Last barrel layer with hit for true muons', layerBins, layerMin, layerMax)
244 barrelLayer_pi = TH1F('BarrelLayer-pions', 'Last barrel layer with hit for true pions', layerBins, layerMin, layerMax)
245 endcapLayer_mu = TH1F('EndcapLayer-muons', 'Last endcap layer with hit for true muons', layerBins, layerMin, layerMax)
246 endcapLayer_pi = TH1F('EndcapLayer-pions', 'Last endcap layer with hit for true pions', layerBins, layerMin, layerMax)
247
248 diffLayerBins = 20
249 diffLayerMin = -0.5
250 diffLayerMax = 19.5
251 diffLayer_mu = TH1F('DiffLayer-muons', 'Range difference for true muons', diffLayerBins, diffLayerMin, diffLayerMax)
252 diffLayer_pi = TH1F('DiffLayer-pions', 'Range difference for true pions', diffLayerBins, diffLayerMin, diffLayerMax)
253
254 ndofBins = 32
255 ndofMin = -0.5
256 ndofMax = 31.5
257 ndof_mu = TH1F('Ndof-muons', 'Degrees of freedom for true muons', ndofBins, ndofMin, ndofMax)
258 ndof_pi = TH1F('Ndof-pions', 'Degrees of freedom for true pions', ndofBins, ndofMin, ndofMax)
259
260 rchisqBins = 50
261 rchisqMin = 0.0
262 rchisqMax = 10.0
263 rchisq_mu = TH1F('rchisq-muons', 'Reduced chi-squared for true muons', rchisqBins, rchisqMin, rchisqMax)
264 rchisq_pi = TH1F('rchisq-pions', 'Reduced chi-squared for true pions', rchisqBins, rchisqMin, rchisqMax)
265
266 momentumBins = 25 # for efficiency
267 momentumBins2 = 5 # for fake rate (low statistics)
268 momentumMin = 0.0
269 momentumMax = 5.0
270 efficiency_momentum = TH1F('Eff-momentum', 'Muon efficiency vs momentum', momentumBins, momentumMin, momentumMax)
271 efficiency_momentum_denom = TH1F('Eff-momentum-denom', 'Muon efficiency vs momentum', momentumBins, momentumMin, momentumMax)
272 fakerate_momentum = TH1F('FakeRate-momentum', 'Pion fake rate vs momentum', momentumBins2, momentumMin, momentumMax)
273 fakerate_momentum_denom = TH1F('FakeRate-momentum-denom', 'Pion fake rate vs momentum', momentumBins2, momentumMin, momentumMax)
274
275 thetaBins = 35 # for efficiency
276 thetaBins2 = 7 # for fake rate (low statistics)
277 thetaMin = 10.0
278 thetaMax = 150.0
279 efficiency_theta = TH1F('Eff-theta', 'Muon efficiency vs theta', thetaBins, thetaMin, thetaMax)
280 efficiency_theta_denom = TH1F('Eff-theta-denom', 'Muon efficiency vs theta', thetaBins, thetaMin, thetaMax)
281 fakerate_theta = TH1F('FakeRate-theta', 'Pion fake rate vs theta', thetaBins2, thetaMin, thetaMax)
282 fakerate_theta_denom = TH1F('FakeRate-theta-denom', 'Pion fake rate vs theta', thetaBins2, thetaMin, thetaMax)
283
284 phiBins = 36 # for efficiency
285 phiBins2 = 8 # for fake rate (low statistics)
286 phiMin = 0.0
287 phiMax = 360.0
288 efficiency_phi = TH1F('Eff-phi', 'Muon efficiency vs phi', phiBins, phiMin, phiMax)
289 efficiency_phi_denom = TH1F('Eff-phi-denom', 'Muon efficiency vs phi', phiBins, phiMin, phiMax)
290 fakerate_phi = TH1F('FakeRate-phi', 'Pion fake rate vs phi', phiBins2, phiMin, phiMax)
291 fakerate_phi_denom = TH1F('FakeRate-phi-denom', 'Pion fake rate vs phi', phiBins2, phiMin, phiMax)
292
293 for chainIndex in range(0, file_chain.GetEntries()):
294 file_chain.GetEntry(chainIndex)
295 if file_chain.GetBranch('MCParticles').GetNdata() == 0:
296 continue
297 # only first (ParticleGun) entry
298 px = file_chain.GetLeaf('MCParticles.m_momentum_x').GetValue(0)
299 py = file_chain.GetLeaf('MCParticles.m_momentum_y').GetValue(0)
300 pz = file_chain.GetLeaf('MCParticles.m_momentum_z').GetValue(0)
301 momentum = ROOT.Math.XYZVector(px, py, pz)
302 isMuon = 'muon' in file_chain.GetFile().GetName()
303 for i in range(0, file_chain.GetBranch('KLMMuidLikelihoods').GetNdata()):
304 o = file_chain.GetLeaf('KLMMuidLikelihoods.m_Outcome').GetValue(i)
305 outcome.Fill(o)
306 if o > 0:
307 # ChargedStable indices: 0=electron, 1=muon, 2=pion, 3=kaon, 4=proton, 5=deuteron
308 llMu = file_chain.GetLeaf('KLMMuidLikelihoods.m_LogL').GetValue(6*i+1)
309 llPi = file_chain.GetLeaf('KLMMuidLikelihoods.m_LogL').GetValue(6*i+2)
310 blayer = file_chain.GetLeaf('KLMMuidLikelihoods.m_BarrelHitLayer').GetValue(i)
311 elayer = file_chain.GetLeaf('KLMMuidLikelihoods.m_EndcapHitLayer').GetValue(i)
312 diffLayer = file_chain.GetLeaf('KLMMuidLikelihoods.m_ExtLayer').GetValue(i) - \
313 file_chain.GetLeaf('KLMMuidLikelihoods.m_HitLayer').GetValue(i)
314 ndof = file_chain.GetLeaf('KLMMuidLikelihoods.m_DegreesOfFreedom').GetValue(i)
315 chisq = file_chain.GetLeaf('KLMMuidLikelihoods.m_ChiSquared').GetValue(i)
316 rchisq = -1.0
317 if ndof > 0:
318 rchisq = chisq / ndof
319 p = momentum.R()
320 theta = momentum.Theta() * 180.0 / np.pi
321 phi = momentum.Phi() * 180.0 / np.pi
322 if phi < 0.0:
323 phi = phi + 360.0
324 if isMuon:
325 llMu_mu.Fill(llMu)
326 llPi_mu.Fill(llPi)
327 llDiff_mu.Fill(llMu - llPi)
328 barrelLayer_mu.Fill(blayer)
329 endcapLayer_mu.Fill(elayer)
330 diffLayer_mu.Fill(diffLayer)
331 ndof_mu.Fill(ndof)
332 rchisq_mu.Fill(rchisq)
333 efficiency_momentum_denom.Fill(p)
334 efficiency_theta_denom.Fill(theta)
335 efficiency_phi_denom.Fill(phi)
336 if llMu > llPi:
337 efficiency_momentum.Fill(p)
338 efficiency_theta.Fill(theta)
339 efficiency_phi.Fill(phi)
340 else:
341 llMu_pi.Fill(llMu)
342 llPi_pi.Fill(llPi)
343 llDiff_pi.Fill(llMu - llPi)
344 barrelLayer_pi.Fill(blayer)
345 endcapLayer_pi.Fill(elayer)
346 diffLayer_pi.Fill(diffLayer)
347 ndof_pi.Fill(ndof)
348 rchisq_pi.Fill(rchisq)
349 fakerate_momentum_denom.Fill(p)
350 fakerate_theta_denom.Fill(theta)
351 fakerate_phi_denom.Fill(phi)
352 if llMu > llPi:
353 fakerate_momentum.Fill(p)
354 fakerate_theta.Fill(theta)
355 fakerate_phi.Fill(phi)
356
357 for j in range(efficiency_momentum_denom.GetNbinsX()):
358 num = efficiency_momentum.GetBinContent(j + 1)
359 denom = efficiency_momentum_denom.GetBinContent(j + 1)
360 if denom > 0:
361 efficiency_momentum.SetBinContent(j + 1, num / denom)
362 efficiency_momentum.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
363 for j in range(efficiency_theta_denom.GetNbinsX()):
364 num = efficiency_theta.GetBinContent(j + 1)
365 denom = efficiency_theta_denom.GetBinContent(j + 1)
366 if denom > 0:
367 efficiency_theta.SetBinContent(j + 1, num / denom)
368 efficiency_theta.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
369 for j in range(efficiency_phi_denom.GetNbinsX()):
370 num = efficiency_phi.GetBinContent(j + 1)
371 denom = efficiency_phi_denom.GetBinContent(j + 1)
372 if denom > 0:
373 efficiency_phi.SetBinContent(j + 1, num / denom)
374 efficiency_phi.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
375 for j in range(fakerate_momentum_denom.GetNbinsX()):
376 num = fakerate_momentum.GetBinContent(j + 1)
377 denom = fakerate_momentum_denom.GetBinContent(j + 1)
378 if denom > 0:
379 fakerate_momentum.SetBinContent(j + 1, num / denom)
380 fakerate_momentum.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
381 for j in range(fakerate_theta_denom.GetNbinsX()):
382 num = fakerate_theta.GetBinContent(j + 1)
383 denom = fakerate_theta_denom.GetBinContent(j + 1)
384 if denom > 0:
385 fakerate_theta.SetBinContent(j + 1, num / denom)
386 fakerate_theta.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
387 for j in range(fakerate_phi_denom.GetNbinsX()):
388 num = fakerate_phi.GetBinContent(j + 1)
389 denom = fakerate_phi_denom.GetBinContent(j + 1)
390 if denom > 0:
391 fakerate_phi.SetBinContent(j + 1, num / denom)
392 fakerate_phi.SetBinError(j + 1, math.sqrt(num * (denom - num) / (denom * denom * denom)))
393
394 # NOTE: *.Fill() must precede *.GetListOfFunctions().Add() or the latter will be discarded!
395 outcome.GetListOfFunctions().Add(TNamed('Description', "0=not in KLM, 1/2=barrel/endcap stop, 3/4=barrel/endcap exit"))
396 outcome.GetListOfFunctions().Add(TNamed('Check', "Peak at 3; rest ~ flat"))
397 outcome.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
398 outcome.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
399 outcome.SetMinimum(0.0)
400
401 llMu_mu.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood of muon hypothesis for true muons"))
402 llMu_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at 0 with tail on negative side"))
403 llMu_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
404 llMu_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
405 llMu_pi.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood of muon hypothesis for true pions"))
406 llMu_pi.GetListOfFunctions().Add(TNamed('Check', "Sharp double peak at ~52 with tail up to 0"))
407 llMu_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
408 llMu_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
409 llPi_mu.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood of pion hypothesis for true muons"))
410 llPi_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at ~-42 with tail up to 0"))
411 llPi_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
412 llPi_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
413 llPi_pi.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood of pion hypothesis for true pions"))
414 llPi_pi.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at 0 with tail on negative side"))
415 llPi_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
416 llPi_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
417 llDiff_mu.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood difference of muon-pion hypotheses for true muons"))
418 llDiff_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at ~+42 with tail on lower side"))
419 llDiff_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
420 llDiff_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'shifter,pvalue-warn=0.50,pvalue-error=0.10'))
421 llDiff_pi.GetListOfFunctions().Add(TNamed('Description', "Log-likelihood difference of muon-pion hypotheses for true pions"))
422 llDiff_pi.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at ~-52 with tail on higher side"))
423 llDiff_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
424 llDiff_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'shifter,pvalue-warn=0.50,pvalue-error=0.10'))
425
426 barrelLayer_mu.GetListOfFunctions().Add(TNamed('Description', "Outermost barrel layer with matching hit for true muons"))
427 barrelLayer_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at layer 14"))
428 barrelLayer_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
429 barrelLayer_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
430 barrelLayer_pi.GetListOfFunctions().Add(TNamed('Description', "Outermost barrel layer with matching hit for true pions"))
431 barrelLayer_pi.GetListOfFunctions().Add(TNamed('Check', "Peak at layers 0-1 with tail above that"))
432 barrelLayer_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
433 barrelLayer_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
434 endcapLayer_mu.GetListOfFunctions().Add(TNamed('Description', "Outermost endcap layer with matching hit for true muons"))
435 endcapLayer_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peaks at layers 11 (backward) and 13 (forward)"))
436 endcapLayer_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
437 endcapLayer_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
438 endcapLayer_pi.GetListOfFunctions().Add(TNamed('Description', "Outermost endcap layer with matching hit for true pions"))
439 endcapLayer_pi.GetListOfFunctions().Add(TNamed('Check', "Peak at layers 0-1 with tail above that"))
440 endcapLayer_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
441 endcapLayer_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
442
443 diffLayer_mu.GetListOfFunctions().Add(TNamed('Description', "Outermost extrapolated-vs-hit layer difference for true muons"))
444 diffLayer_mu.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at 0"))
445 diffLayer_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
446 diffLayer_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'shifter,pvalue-warn=0.50,pvalue-error=0.10'))
447 diffLayer_mu.GetXaxis().SetTitle('Last Ext Layer - Last Hit Layer')
448 diffLayer_pi.GetListOfFunctions().Add(TNamed('Description', "Outermost extrapolated-vs-hit layer difference for true pions"))
449 diffLayer_pi.GetListOfFunctions().Add(TNamed('Check', "Sharp peak at 15 with tail below that"))
450 diffLayer_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
451 diffLayer_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'shifter,pvalue-warn=0.50,pvalue-error=0.10'))
452 diffLayer_pi.GetXaxis().SetTitle('Last Ext Layer - Last Hit Layer')
453
454 ndof_mu.GetListOfFunctions().Add(TNamed('Description', "(Number of matching-hit layers) x 2 for true muons"))
455 ndof_mu.GetListOfFunctions().Add(TNamed('Check', "Even-only entries; peak at 30 with tail below"))
456 ndof_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
457 ndof_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
458 ndof_mu.GetXaxis().SetTitle('#dof')
459 ndof_pi.GetListOfFunctions().Add(TNamed('Description', "(Number of matching-hit layers) x 2 for true pions"))
460 ndof_pi.GetListOfFunctions().Add(TNamed('Check', "Even-only entries; sharp peak at 0 with tail above"))
461 ndof_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
462 ndof_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
463 ndof_pi.GetXaxis().SetTitle('#dof')
464
465 rchisq_mu.GetListOfFunctions().Add(TNamed('Description', "Transverse-deviation chi-squared/#dof for true muons"))
466 rchisq_mu.GetListOfFunctions().Add(TNamed('Check', "Narrow peak at 1"))
467 rchisq_mu.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
468 rchisq_mu.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
469 rchisq_mu.GetXaxis().SetTitle('#chi^{2} / #ndof')
470 rchisq_pi.GetListOfFunctions().Add(TNamed('Description', "Transverse-deviation chi-squared/#dof for true pions"))
471 rchisq_pi.GetListOfFunctions().Add(TNamed('Check', "Broad peak at ~1"))
472 rchisq_pi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
473 rchisq_pi.GetListOfFunctions().Add(TNamed('MetaOptions', 'expert,pvalue-warn=0.50,pvalue-error=0.10'))
474 rchisq_pi.GetXaxis().SetTitle('#chi^{2} / #ndof')
475
476 efficiency_momentum.GetListOfFunctions().Add(TNamed('Description', "Muon ID efficiency vs p for logL(mu) > logL(pi)"))
477 efficiency_momentum.GetListOfFunctions().Add(TNamed('Check', "90-95% for p > 0.6 GeV/c"))
478 efficiency_momentum.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
479 efficiency_momentum.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
480 efficiency_momentum.GetXaxis().SetTitle('p (GeV/#it{c})')
481 efficiency_momentum.SetMinimum(0.0)
482 fakerate_momentum.GetListOfFunctions().Add(TNamed('Description', "Pion fake rate vs p for logL(mu) > logL(pi)"))
483 fakerate_momentum.GetListOfFunctions().Add(TNamed('Check', "Peak of 11% below 1 GeV/c, 2-6% above 1 GeV/c"))
484 fakerate_momentum.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
485 fakerate_momentum.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
486 fakerate_momentum.GetXaxis().SetTitle('p (GeV/#it{c})')
487 fakerate_momentum.SetMinimum(0.0)
488
489 efficiency_theta.GetListOfFunctions().Add(TNamed('Description', "Muon ID efficiency vs theta for logL(mu) > logL(pi)"))
490 efficiency_theta.GetListOfFunctions().Add(TNamed('Check', "Somewhat flat (80-100%) in KLM range (17 < theta < 145 degrees)"))
491 efficiency_theta.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
492 efficiency_theta.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
493 efficiency_theta.GetXaxis().SetTitle('#theta (deg)')
494 efficiency_theta.SetMinimum(0.0)
495 fakerate_theta.GetListOfFunctions().Add(TNamed('Description', "Pion fake rate vs theta for logL(mu) > logL(pi)"))
496 fakerate_theta.GetListOfFunctions().Add(TNamed('Check', "Roughly flat (3-7%) in KLM range (17 < theta < 145 degrees)"))
497 fakerate_theta.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
498 fakerate_theta.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
499 fakerate_theta.GetXaxis().SetTitle('#theta (deg)')
500 fakerate_theta.SetMinimum(0.0)
501
502 efficiency_phi.GetListOfFunctions().Add(TNamed('Description', "Muon identification efficiency vs phi for logL(mu) > logL(pi)"))
503 efficiency_phi.GetListOfFunctions().Add(TNamed('Check', "Somewhat flat distribution (85-100%)"))
504 efficiency_phi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
505 efficiency_phi.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
506 efficiency_phi.GetXaxis().SetTitle('#phi (deg)')
507 efficiency_phi.SetMinimum(0.0)
508 fakerate_phi.GetListOfFunctions().Add(TNamed('Description', "Pion fake rate vs phi for logL(mu) > logL(pi)"))
509 fakerate_phi.GetListOfFunctions().Add(TNamed('Check', "Roughly flat distribution (2-7%)"))
510 fakerate_phi.GetListOfFunctions().Add(TNamed('Contact', contact_mail))
511 fakerate_phi.GetListOfFunctions().Add(TNamed('MetaOptions', 'pvalue-warn=0.50,pvalue-error=0.10'))
512 fakerate_phi.GetXaxis().SetTitle('#phi (deg)')
513 fakerate_phi.SetMinimum(0.0)
514
515 outcome.Write()
516 llMu_mu.Write()
517 llMu_pi.Write()
518 llPi_mu.Write()
519 llPi_pi.Write()
520 llDiff_mu.Write()
521 llDiff_pi.Write()
522 barrelLayer_mu.Write()
523 barrelLayer_pi.Write()
524 endcapLayer_mu.Write()
525 endcapLayer_pi.Write()
526 diffLayer_mu.Write()
527 diffLayer_pi.Write()
528 ndof_mu.Write()
529 ndof_pi.Write()
530 rchisq_mu.Write()
531 rchisq_pi.Write()
532 efficiency_momentum.Write()
533 efficiency_theta.Write()
534 efficiency_phi.Write()
535 fakerate_momentum.Write()
536 fakerate_theta.Write()
537 fakerate_phi.Write()
538
539
540
543if __name__ == '__main__':
544 if ACTIVE:
545 main()
546 else:
547 print("This validation deactivated and thus basf2 is not executed.\n"
548 "If you want to run this validation, please set the 'ACTIVE' flag above to 'True'.\n"
549 "Exiting.")
main
Definition main.py:1