Belle II Software development
TRGValidation.py
1#!/usr/bin/env python3
2
9
10"""
11<header>
12 <input>TRGValidationGen.root</input>
13 <output>TRGValidation.root</output>
14 <contact>yinjh@nankai.edu.cn</contact>
15 <description>makes validation plots for TRG</description>
16</header>
17
18"""
19
20import basf2
21import modularAnalysis as ma # a shorthand for the analysis tools namespace
22import stdCharged as stdc
23import variables.utils as vu
24
25import ROOT
26from ROOT import Belle2, TH1F, TFile, TNamed, TEfficiency, TMath
27from math import pi as PI
28
29
30inputBits = ["t3_0", "ty_0", "t2_0", "ts_0", "ta_0", "typ", "ehigh", "elow", "elum", "ecl_3dbha", "cdc_open90",
31 "clst_0", "clst_1", "clst_2", "clst_3", "klm_hit", "klm_0", "klm_1", "klm_2", "eklm_0", "eklm_1", "eklm_2"]
32outputBits = ["fff", "ffy", "ffz", "fzo", "fso", "fyo", "ffb", "fsb", "ssb", "stt", "hie", "c4", "bha3d",
33 "lml1", "mu_b2b", "mu_eb2b", "eklm2", "beklm", "cdcklm1", "cdcklm2", "seklm1", "seklm2", "ecleklm1"]
34moniInbits = []
35moniOutbits = []
36
37inputBits_expert = [
38 't3_0',
39 't3_1',
40 't3_2',
41 't3_3',
42 'ty_0',
43 'ty_1',
44 'ty_2',
45 'ty_3',
46 't2_0',
47 't2_1',
48 't2_2',
49 't2_3',
50 'ts_0',
51 'ts_1',
52 'ts_2',
53 'ts_3',
54 'ta_0',
55 'ta_1',
56 'ta_2',
57 'ta_3',
58 'typ',
59 'typ4',
60 'typ5',
61 'cdc_open90',
62 'cdc_active',
63 'cdc_b2b3',
64 'cdc_b2b5',
65 'cdc_b2b7',
66 'cdc_b2b9',
67 'itsfb2b',
68 'ti',
69 'i2io',
70 'i2fo',
71 'f2f30',
72 's2f30',
73 's2s30',
74 's2s3',
75 's2s5',
76 's2so',
77 's2f3',
78 's2f5',
79 's2fo',
80 'fwd_s',
81 'bwd_s',
82 'track',
83 'trkflt',
84 'ehigh',
85 'elow',
86 'elum',
87 'ecl_bha',
88 'ecl_3dbha',
89 'bha_veto',
90 'bha_type_0',
91 'bha_type_1',
92 'bha_intrk',
93 'bha_theta_0',
94 'bha_theta_1',
95 'clst_0',
96 'clst_1',
97 'clst_2',
98 'clst_3',
99 'ecl_active',
100 'ecl_timing_fwd',
101 'ecl_timing_brl',
102 'ecl_timing_bwd',
103 'ecl_phys',
104 'ecl_oflo',
105 'ecl_lml_0',
106 'ecl_lml_1',
107 'ecl_lml_2',
108 'ecl_lml_3',
109 'ecl_lml_4',
110 'ecl_lml_5',
111 'ecl_lml_6',
112 'ecl_lml_7',
113 'ecl_lml_8',
114 'ecl_lml_9',
115 'ecl_lml_10',
116 'ecl_lml_12',
117 'ecl_lml_13',
118 'ecl_mumu',
119 'ecl_bhapur',
120 'ecl_bst',
121 'top_0',
122 'top_1',
123 'top_2',
124 'top_bb',
125 'top_active',
126 'klm_hit',
127 'klm_0',
128 'klm_1',
129 'klm_2',
130 'klmb2b',
131 'eklm_hit',
132 'eklm_0',
133 'eklm_1',
134 'eklm_2',
135 'eklmb2b',
136 'revo',
137 'her_kick',
138 'ler_kick',
139 'bha_delay',
140 'pseud_rand',
141 'plsin',
142 'poissonin',
143 'veto',
144 'injv',
145 'secl',
146 'iecl_0',
147 'iecl_1',
148 'samhem',
149 'opohem',
150 'd3',
151 'd5',
152 'd7',
153 'p3',
154 'p5',
155 'p7',
156 'typ6',
157 'cdcecl_0',
158 'cdcecl_1',
159 'cdcecl_2',
160 'cdcecl_3',
161 'c2gev_0',
162 'c2gev_1',
163 'c2gev_2',
164 'c2gev_3',
165 'cdctop_0',
166 'cdctop_1',
167 'cdctop_2',
168 'cdctop_3',
169 'cdcklm_0',
170 'cdcklm_1',
171 'seklm_0',
172 'seklm_1',
173 'ecleklm',
174 'ieklm',
175 'fwdsb',
176 'bwdsb',
177 'fwdnb',
178 'bwdnb',
179 'brlfb1',
180 'brlfb2',
181 'brlnb1',
182 'brlnb2',
183 'trkbha1',
184 'trkbha2',
185 'grlgg1',
186 'grlgg2',
187 'nimin0',
188 'nimin1',
189 'ecl_taub2b',
190 'ehigh1',
191 'ehigh2',
192 'ehigh3']
193outputBits_expert = [
194 'fff',
195 'ffs',
196 'fss',
197 'sss',
198 'ffz',
199 'fzz',
200 'zzz',
201 'ffy',
202 'fyy',
203 'yyy',
204 'ff',
205 'fs',
206 'ss',
207 'fz',
208 'zz',
209 'fy',
210 'yy',
211 'ffo',
212 'fso',
213 'sso',
214 'fzo',
215 'fyo',
216 'ffb',
217 'fsb',
218 'ssb',
219 'fzb',
220 'fyb',
221 'aaa',
222 'aaao',
223 'aao',
224 'aab',
225 'aa',
226 'hie',
227 'lowe',
228 'lume',
229 'hade',
230 'c2',
231 'c3',
232 'c4',
233 'c5',
234 'bha3d',
235 'bhabha',
236 'bhabha_trk',
237 'bhabha_brl',
238 'bhabha_ecp',
239 'bhapur',
240 'eclmumu',
241 'bhauni',
242 'ecloflo',
243 'eclbst',
244 'g_high',
245 'g_c1',
246 'gg',
247 'eed',
248 'fed',
249 'fp',
250 'sp',
251 'zp',
252 'yp',
253 'd_5',
254 'shem',
255 'ohem',
256 'toptiming',
257 'ecltiming',
258 'cdctiming',
259 'cdcbb',
260 'mu_pair',
261 'mu_b2b',
262 'klmhit',
263 'mu_epair',
264 'mu_eb2b',
265 'eklmhit',
266 'revolution',
267 'random',
268 'bg',
269 'pls',
270 'poisson',
271 'vetout',
272 'f',
273 's',
274 'z',
275 'y',
276 'a',
277 'n1gev1',
278 'n1gev2',
279 'n1gev3',
280 'n1gev4',
281 'n2gev1',
282 'n2gev2',
283 'n2gev3',
284 'n2gev4',
285 'c2gev1',
286 'c2gev2',
287 'c2gev3',
288 'c2gev4',
289 'cdcecl1',
290 'cdcecl2',
291 'cdcecl3',
292 'cdcecl4',
293 'cdcklm1',
294 'cdcklm2',
295 'cdcklm3',
296 'cdcklm4',
297 'cdctop1',
298 'cdctop2',
299 'cdctop3',
300 'cdctop4',
301 'c1hie',
302 'c1lume',
303 'n1hie',
304 'n1lume',
305 'c3hie',
306 'c3lume',
307 'n3hie',
308 'n3lume',
309 'lml0',
310 'lml1',
311 'lml2',
312 'lml3',
313 'lml4',
314 'lml5',
315 'lml6',
316 'lml7',
317 'lml8',
318 'lml9',
319 'lml10',
320 'lml12',
321 'lml13',
322 'zzzv',
323 'yyyv',
324 'fffv',
325 'zzv',
326 'yyv',
327 'ffov',
328 'fffov',
329 'hiev',
330 'lumev',
331 'c4v',
332 'bhabhav',
333 'mu_pairv',
334 'bha3dv',
335 'bffo',
336 'bhie',
337 'ffoc',
338 'ffoc2',
339 'fffo',
340 'ff30',
341 'fs30',
342 'ss30',
343 'itsf_b2b',
344 'gggrl',
345 'ggtsf',
346 'ggbrl',
347 'bhabrl',
348 'bhamtc1',
349 'bhamtc2',
350 'bhaf',
351 'nim0',
352 'nima01',
353 'nimo01']
354moniInbits_expert = []
355moniOutbits_expert = []
356
357
358class MakePlots(basf2.Module):
359 '''
360 Make validation histograms for trg ecl/cdc/klm
361 '''
362
363 def set_descr_shifter(self, histogram, description, check):
364 '''
365 Sets description, check and contact to validation histogram.
366 :param histogram validation histogram
367 :param description description text
368 :param check information on what to check in comparison with reference
369 '''
370
371 descr = TNamed("Description", description)
372 histogram.GetListOfFunctions().Add(descr)
373 Check = TNamed("Check", check)
374 histogram.GetListOfFunctions().Add(Check)
375 contact = TNamed("Contact", "yinjh@nankai.edu.cn")
376 histogram.GetListOfFunctions().Add(contact)
377 Meta = TNamed("MetaOptions", "shifter,pvalue-warn=0.02,pvalue-error=0.01")
378 histogram.GetListOfFunctions().Add(Meta)
379
380 def set_descr_expert(self, histogram, description, check):
381 '''
382 Sets description, check and contact to validation histogram.
383 :param histogram validation histogram
384 :param description description text
385 :param check information on what to check in comparison with reference
386 '''
387
388 descr = TNamed("Description", description)
389 histogram.GetListOfFunctions().Add(descr)
390 Check = TNamed("Check", check)
391 histogram.GetListOfFunctions().Add(Check)
392 contact = TNamed("Contact", "yinjh@nankai.edu.cn")
393 histogram.GetListOfFunctions().Add(contact)
394 Meta = TNamed("MetaOptions", "expert,pvalue-warn=0.02,pvalue-error=0.01")
395 histogram.GetListOfFunctions().Add(Meta)
396
397 def set_style(self, histogram, xtitle, ytitle):
398 '''
399 Sets x-y titles, and sets histogram style.
400 :param histogram validation histogram
401 :param xtitle X-axis title
402 :param xtitle Y-axis title
403 '''
404
405 histogram.GetXaxis().SetTitle(xtitle)
406 histogram.GetYaxis().SetTitle(ytitle)
407 histogram.GetYaxis().SetTitleOffset(1.4)
408 histogram.GetXaxis().SetTitleSize(0.045)
409 histogram.GetYaxis().SetLabelSize(0.020)
410 histogram.SetLineColor(ROOT.kBlack)
411
412 def get_relative(self, hist1, hist2, title, particle, trgbit):
413 '''
414 Get relative ratio between two hists.
415 :param hist1 numerator
416 :param hist2 denominator
417 :param title new histogram title
418 :param particle particle name
419 :param trgbit trigger bit
420 '''
421
422 bin1 = hist1.GetXaxis().GetNbins()
423 bin2 = hist2.GetXaxis().GetNbins()
424 Xmin1 = hist1.GetXaxis().GetXmin()
425 Xmax1 = hist1.GetXaxis().GetXmax()
426 Xmin2 = hist2.GetXaxis().GetXmin()
427 Xmax2 = hist2.GetXaxis().GetXmax()
428 if bin1 != bin2 or Xmin1 != Xmin2 or Xmax1 != Xmax2:
429 print("warning: two histograms cannot be compared!!")
430 return 0
431
432 teff = TEfficiency(hist1, hist2)
433 htitle = f'h_eff_{title}_{particle}_{trgbit}'
434 heff = TH1F(htitle, f"efficiency histogram of {particle} for {trgbit}", bin1, Xmin1, Xmax1)
435 for ibin in range(1, bin1+1):
436 binc0 = teff.GetEfficiency(ibin)
437 bine0 = (teff.GetEfficiencyErrorUp(ibin) + teff.GetEfficiencyErrorLow(ibin)) / 2.0
438 heff.SetBinContent(ibin, binc0)
439 heff.SetBinError(ibin, bine0)
440 heff.GetYaxis().SetRangeUser(0.0, 1.2)
441 return heff
442
443 def initialize(self):
444 '''
445 initialize class members
446 '''
447
448 print('start read')
449
450 self.tfile = TFile("./TRGValidation.root", "update")
451
452 self.Nevent = 0
453
454 m_dbinput = Belle2.PyDBObj("TRGGDLDBInputBits")
455 m_dbftdl = Belle2.PyDBObj("TRGGDLDBFTDLBits")
456 for i in range(320):
457 inbitname = m_dbinput.getinbitname(i)
458 outbitname = m_dbftdl.getoutbitname(i)
459 if inbitname in inputBits:
460 moniInbits.append({"name": inbitname, "index": i})
461 if outbitname in outputBits:
462 moniOutbits.append({"name": outbitname, "index": i})
463 if inbitname in inputBits_expert:
464 moniInbits_expert.append({"name": inbitname, "index": i})
465 if outbitname in outputBits_expert:
466 moniOutbits_expert.append({"name": outbitname, "index": i})
467
468 mc_px = "MCParticles.m_momentum_x"
469 mc_py = "MCParticles.m_momentum_y"
470 trk2d_omega = "TRGCDC2DFinderTracks.m_omega"
471 trk2d_phi = "TRGCDC2DFinderTracks.m_phi0"
472
473 ROOT.gROOT.SetBatch(True)
474 ROOT.gStyle.SetOptStat(ROOT.kFALSE)
475
476
477 n_inbit_test = len(inputBits)
478 self.hist_inbit = TH1F("hin", "trigger input bits", n_inbit_test, 0, n_inbit_test)
479 self.hist_inbit.GetXaxis().SetRangeUser(0, n_inbit_test)
480 self.hist_inbit.GetXaxis().SetLabelSize(0.04)
481 self.set_descr_shifter(self.hist_inbit, "trigger input bits for shifter",
482 "Efficiency should not drop significantly for any trigger bit")
483 self.set_style(self.hist_inbit, "", "Efficiency")
484
485
486 n_oubit_test = len(outputBits)
487 self.hist_outbit = TH1F("hout", "trigger output bits", n_oubit_test, 0, n_oubit_test)
488 self.hist_outbit.GetXaxis().SetRangeUser(0, n_oubit_test)
489 self.hist_outbit.GetXaxis().SetLabelSize(0.04)
490 self.set_descr_shifter(
491 self.hist_outbit,
492 "trigger ftdl bits for shifter",
493 "Efficiency should not drop significantly. For some output bits the efficiency is very low, close to zero.")
494 self.set_style(self.hist_outbit, "", "Efficiency")
495
496 for i in range(n_inbit_test):
497 self.hist_inbit.GetXaxis().SetBinLabel(
498 self.hist_inbit.GetXaxis().FindBin(i + 0.5), inputBits[i])
499 for i in range(n_oubit_test):
500 self.hist_outbit.GetXaxis().SetBinLabel(
501 self.hist_outbit.GetXaxis().FindBin(i + 0.5), outputBits[i])
502
503
504 self.h_E_ECL = TH1F("h_E_ECL", "ECL cluster energy [50 MeV]", 100, 0, 5)
505 self.set_descr_shifter(self.h_E_ECL, "TRG ECL cluster energy", "Exponentially falling distribution")
506 self.set_style(self.h_E_ECL, "ECL cluster energy (GeV)", "Events/(50 MeV)")
507
508
509 self.h_Esum_ECL = TH1F("h_Esum_ECL", "sum of ECL cluster energy [50 MeV]", 100, 0, 5)
510 self.set_descr_shifter(self.h_Esum_ECL, "sum of TRG ECL cluster energy", "Peak at 200 MeV with tail")
511 self.set_style(self.h_Esum_ECL, "sum of ECL cluster energy (GeV)", "Events/(50 MeV)")
512
513
514 self.h_theta_ECL = TH1F("h_theta_ECL", "TRG ECL cluster #theta [1.4 degrees]", 128, 0, 180)
515 self.set_descr_shifter(self.h_theta_ECL, "TRG ECL cluster polar angle", "uniform in barrel")
516 self.set_style(self.h_theta_ECL, "TRG ECL cluster polar angle (degree)", "Events/(1.4 degree)")
517
518
519 self.h_thetaID_ECL = TH1F("h_thetaID_ECL", "ECL cluster TC ID", 610, 0, 610)
520 self.set_descr_shifter(self.h_thetaID_ECL, "TRG ECL cluster theta ID", "uniform in barrel")
521 self.set_style(self.h_thetaID_ECL, "ECL cluster TC ID", "Events/(1)")
522
523
524 self.h_phi_ECL = TH1F("h_phi_ECL", "TRG ECL cluster phi [2.0 degrees]", 180, -180, 180)
525 self.set_descr_shifter(self.h_phi_ECL, "TRG ECL cluster phi distribution", "distribute uniformly")
526 self.set_style(self.h_phi_ECL, "ECL cluster #phi (degree) ", "Events/(2.0 degrees)")
527
528
529 self.h_sector_BKLM = TH1F("h_sector_BKLM", "BKLM TRG hit sector", 10, 0, 10)
530 self.set_descr_shifter(self.h_sector_BKLM, "# of BKLM TRG sector", "peak at 0")
531 self.set_style(self.h_sector_BKLM, "# of BKLM TRG sector", "Events/(1)")
532
533
534 self.h_sector_EKLM = TH1F("h_sector_EKLM", "EKLM TRG hit sector", 10, 0, 10)
535 self.set_descr_shifter(self.h_sector_EKLM, "# of EKLM TRG sector", "peak at 0")
536 self.set_style(self.h_sector_EKLM, "# of EKLM TRG sector", "Events/(1)")
537
538
539 n_inbit_expert = len(inputBits_expert)
540 self.hist_inbit_expert = TH1F("hin_expert", "trigger input bits", n_inbit_expert, 0, n_inbit_expert)
541 self.hist_inbit_expert.GetXaxis().SetRangeUser(0, n_inbit_expert)
542 self.set_descr_expert(self.hist_inbit_expert, "# of EKLM TRG sector", "peak at 0")
543 self.set_style(self.hist_inbit_expert, "", "Efficiency")
544
545
546 n_oubit_expert = len(outputBits_expert)
547 self.hist_outbit_expert = TH1F("hout_expert", "trigger output bits", n_oubit_expert, 0, n_oubit_expert)
548 self.hist_outbit_expert.GetXaxis().SetRangeUser(0, n_oubit_expert)
549 self.set_descr_expert(self.hist_outbit_expert, "# of EKLM TRG sector", "peak at 0")
550 self.set_style(self.hist_outbit_expert, "", "Efficiency")
551
552 for i in range(n_inbit_expert):
553 self.hist_inbit_expert.GetXaxis().SetBinLabel(
554 self.hist_inbit_expert.GetXaxis().FindBin(i + 0.5), inputBits_expert[i])
555 for i in range(n_oubit_expert):
556 self.hist_outbit_expert.GetXaxis().SetBinLabel(
557 self.hist_outbit_expert.GetXaxis().FindBin(i + 0.5), outputBits_expert[i])
558
559 mc = "abs(MCParticles.m_pdg)==11&&MCParticles.m_status==11"
560 tree = ROOT.TChain("tree")
561 tree.Add("../TRGValidationGen.root")
562
563
564 self.d_w = TH1F("d_w", "#Deltaw of CDC 2D finder, w = 0.00449/p_{t}", 50, -0.02, 0.02)
565 tree.Draw(f"({trk2d_omega} - 0.00449)/sqrt({mc_px}*{mc_px} + {mc_py}*{mc_py})>> d_w",
566 "MCParticles.m_pdg<0&&" + mc)
567
568 self.d_w_2 = TH1F("d_w_2", "d_w_2", 50, -0.02, 0.02)
569 tree.Draw(f"({trk2d_omega} + 0.00449)/sqrt({mc_px}*{mc_px} + {mc_py}*{mc_py}) >> d_w_2",
570 "MCParticles.m_pdg>0 && " + mc)
571 self.d_w.Add(self.d_w_2)
572 self.set_descr_shifter(self.d_w, "Comparison on w (=0.00449/pt) of a track between CDC 2D finder output and MC.",
573 "A clear peak at 0 with tail.")
574 self.set_style(self.d_w, "#Deltaw", "Events/(0.08)")
575
576
577 self.d_phi = TH1F("d_phi", "#Delta#phi of CDC 2D finder", 50, -0.5, 0.5)
578 tree.Draw(f"{trk2d_phi}-atan({mc_py}/{mc_px})>>d_phi",
579 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
580 f"fabs({trk2d_phi}-atan({mc_py}/{mc_px}))<{PI}&&" + mc)
581
582
583 self.d_phi_2 = TH1F("d_phi_2", "d_phi_2", 50, -0.5, 0.5)
584 tree.Draw(f"{trk2d_phi}-atan({mc_py}/{mc_px})-{PI}>>d_phi_2",
585 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
586 f"{trk2d_phi}-atan({mc_py}/{mc_px})>={PI} &&" + mc)
587
588
589 self.d_phi_3 = TH1F("d_phi_3", "d_phi_3", 50, -0.5, 0.5)
590 tree.Draw(f"{trk2d_phi}-atan({mc_py}/{mc_px})+{PI}>>d_phi_2",
591 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
592 f"{trk2d_phi}-atan({mc_py}/{mc_px})<=-{PI}&&" + mc)
593
594 self.d_phi.Add(self.d_phi_2)
595 self.d_phi.Add(self.d_phi_3)
596 self.set_descr_shifter(self.d_phi, "Comparison on phi_i of a track between CDC 2D finder output and MC.",
597 "A Gaussian peak at 0.")
598 self.set_style(self.d_phi, "#Delta#phi [rad]", "Events/(0.02 rad)")
599
600
601 self.d_z0_3d = TH1F("d_z0_3d", "#Deltaz0 of CDC 3D fitter", 60, -30, 30)
602 tree.Draw("TRGCDC3DFitterTracks.m_z0-MCParticles.m_productionVertex_z>>d_z0_3d", mc)
603 self.set_descr_shifter(self.d_z0_3d, "Comparison on z0 of a track between CDC 2D fitter output and MC.",
604 "A Gaussian peak at 0 with small tail.")
605 self.set_style(self.d_z0_3d, "#Deltaz0 [cm]", "Events/(1 cm)")
606
607
608 self.d_z0_nn = TH1F("d_z0_nn", "#Deltaz0 of CDC Neuro", 60, -30, 30)
609 tree.Draw("TRGCDCNeuroTracks.m_z0-MCParticles.m_productionVertex_z>>d_z0_nn", mc)
610 self.set_descr_shifter(self.d_z0_nn, "Comparison on z0 of a track between CDC Neuro output and MC.",
611 "A Gaussian peak at 0 with small tail.")
612 self.set_style(self.d_z0_nn, "#Deltaz0 [cm]", "Events/(1 cm)")
613
614
615 self.d_E_ECL = TH1F("d_E_ECL", "#DeltaE of ECL clustering", 100, -6, 6)
616 tree.Draw("TRGECLClusters.m_edep-MCParticles.m_energy>>d_E_ECL", mc)
617 self.set_descr_shifter(self.d_E_ECL, "Comparison on deposit energy between ECL cluster output and MC.",
618 "A peak around -0.5 ~ 0 with a tail toward -6.")
619 self.set_style(self.d_E_ECL, "#DeltaE [GeV]", "Events/(0.12 GeV)")
620
621 def beginRun(self):
622 '''
623 reset all histograms at the begin of a new run
624 '''
625
626 self.hist_inbit.Reset()
627 self.hist_outbit.Reset()
628 self.hist_inbit_expert.Reset()
629 self.hist_outbit_expert.Reset()
630 self.h_Esum_ECL.Reset()
631 self.h_E_ECL.Reset()
632 self.h_theta_ECL.Reset()
633 self.h_thetaID_ECL.Reset()
634 self.h_phi_ECL.Reset()
635
636 def event(self):
637 '''
638 event loop
639 '''
640
641 trg_summary = Belle2.PyStoreObj("TRGSummary")
642 clusters = Belle2.PyStoreArray("TRGECLClusters")
643 klmSummary = Belle2.PyStoreObj("KLMTrgSummary")
644
645 for bit in moniInbits:
646 binIndex = inputBits.index(bit["name"])
647 bitIndex = bit["index"]
648 if trg_summary.testInput(bitIndex):
649 self.hist_inbit.Fill(binIndex + 0.5)
650 for bit in moniOutbits:
651 binIndex = outputBits.index(bit["name"])
652 bitIndex = bit["index"]
653 if trg_summary.testFtdl(bitIndex):
654 self.hist_outbit.Fill(binIndex + 0.5)
655
656 for bit in moniInbits_expert:
657 binIndex = inputBits_expert.index(bit["name"])
658 bitIndex = bit["index"]
659 if trg_summary.testInput(bitIndex):
660 self.hist_inbit_expert.Fill(binIndex + 0.5)
661 for bit in moniOutbits_expert:
662 binIndex = outputBits_expert.index(bit["name"])
663 bitIndex = bit["index"]
664 if trg_summary.testFtdl(bitIndex):
665 self.hist_outbit_expert.Fill(binIndex + 0.5)
666
667 etot = 0
668 for cluster in clusters:
669 x = cluster.getPositionX()
670 y = cluster.getPositionY()
671 z = cluster.getPositionZ()
672 e = cluster.getEnergyDep()
673 self.h_E_ECL.Fill(e)
674 etot += e
675 vec = ROOT.Math.XYZVector(x, y, z)
676 self.h_theta_ECL.Fill(vec.Theta() * TMath.RadToDeg())
677 self.h_thetaID_ECL.Fill(cluster.getMaxTCId())
678 self.h_phi_ECL.Fill(vec.Phi() * TMath.RadToDeg())
679
680 if etot > 0:
681 self.h_Esum_ECL.Fill(etot)
682
683 self.h_sector_BKLM.Fill(klmSummary.getBKLM_n_trg_sectors())
684 self.h_sector_EKLM.Fill(klmSummary.getEKLM_n_trg_sectors())
685
686 self.Nevent = self.Nevent + 1
687
688 def endRun(self):
689 '''
690 end of run
691 '''
692
693 print("end")
694
695 def terminate(self):
696 '''
697 write histograms
698 '''
699
700 bits = ['ty_0', 'ehigh', 'clst_0']
701 part = 'e'
702 h_gen_p = self.tfile.Get(f"{part}_all/p")
703 h_gen_E = self.tfile.Get(f"{part}_all/clusterE")
704 h_gen_theta = self.tfile.Get(f"{part}_all/theta")
705 h_gen_phi = self.tfile.Get(f"{part}_all/phi")
706
707 for bit in bits:
708 h_p = self.tfile.Get(f"{part}_{bit}/p")
709 h_E = self.tfile.Get(f"{part}_{bit}/clusterE")
710 h_theta = self.tfile.Get(f"{part}_{bit}/theta")
711 h_phi = self.tfile.Get(f"{part}_{bit}/phi")
712
713 h_eff_p = self.get_relative(h_p, h_gen_p, 'p', part, bit)
714 self.set_style(h_eff_p, "Momentum", "Efficiency")
715 self.set_descr_expert(h_eff_p, "Momentum dependent efficiency for experts", "")
716
717 h_eff_E = self.get_relative(h_E, h_gen_E, 'E', part, bit)
718 self.set_style(h_eff_E, "Energy", "Efficiency")
719 if bit == 'ehigh':
720 self.set_descr_shifter(h_eff_E, "Energy dependent efficiency for shifter",
721 "Efficiency around 40% below 1 GeV, then jump up to 100%")
722 elif bit == 'clst_0':
723 self.set_descr_shifter(h_eff_E, "Energy dependent efficiency for shifter",
724 "Turn-on curve from 40% efficiency at 0.5 GeV to nearly 100% above 1.5 GeV")
725 else:
726 self.set_descr_expert(h_eff_E, "Energy dependent efficiency for experts", "")
727
728 h_eff_theta = self.get_relative(h_theta, h_gen_theta, 'theta', part, bit)
729 self.set_style(h_eff_theta, "Polar angle", "Efficiency")
730 self.set_descr_expert(h_eff_theta, "polar angle dependent efficiency for experts", "")
731
732 h_eff_phi = self.get_relative(h_phi, h_gen_phi, 'phi', part, bit)
733 self.set_style(h_eff_phi, "phi angle", "Efficiency")
734 self.set_descr_expert(h_eff_phi, "azimuth angle dependent efficiency for experts", "")
735
736 h_eff_p.Write()
737 h_eff_E.Write()
738 h_eff_theta.Write()
739 h_eff_phi.Write()
740
741 bits = ['ty_0', 'klm_0', "klm_hit"]
742 part = "mu"
743 h_gen_p = self.tfile.Get(f"{part}_all/p")
744 h_gen_E = self.tfile.Get(f"{part}_all/clusterE")
745 h_gen_theta = self.tfile.Get(f"{part}_all/theta")
746 h_gen_phi = self.tfile.Get(f"{part}_all/phi")
747
748 for bit in bits:
749 h_p = self.tfile.Get(f"{part}_{bit}/p")
750 h_theta = self.tfile.Get(f"{part}_{bit}/theta")
751 h_phi = self.tfile.Get(f"{part}_{bit}/phi")
752
753 h_eff_p = self.get_relative(h_p, h_gen_p, 'p', part, bit)
754 self.set_style(h_eff_p, "Momentum", "Efficiency")
755 if bit == 'ty_0':
756 self.set_descr_shifter(h_eff_p, "Momentum dependent efficiency for shifter",
757 "Efficiency should be above 90% for the whole momentum range")
758 elif bit == 'klm_0':
759 self.set_descr_shifter(h_eff_p, "Momentum dependent efficiency for shifter",
760 "Efficiency should rise to about 90% for momenta greater than 1 GeV")
761 else:
762 self.set_descr_expert(h_eff_p, "Momentum dependent efficiency for experts", "")
763
764 h_eff_theta = self.get_relative(h_theta, h_gen_theta, 'theta', part, bit)
765 self.set_style(h_eff_theta, "Polar angle", "Efficiency")
766 self.set_descr_expert(h_eff_theta, "polar angle dependent efficiency for experts", "")
767
768 h_eff_phi = self.get_relative(h_phi, h_gen_phi, 'phi', part, bit)
769 self.set_style(h_eff_phi, "phi angle", "Efficiency")
770 self.set_descr_expert(h_eff_phi, "azimuth angle dependent efficiency for experts", "")
771
772 h_eff_p.Write()
773 h_eff_theta.Write()
774 h_eff_phi.Write()
775
776 self.hist_inbit.Scale(1./self.Nevent)
777 self.hist_outbit.Scale(1.0/self.Nevent)
778 self.hist_inbit_expert.Scale(1./self.Nevent)
779 self.hist_outbit_expert.Scale(1.0/self.Nevent)
780
781 self.hist_inbit.Write()
782 self.hist_outbit.Write()
783 self.hist_inbit_expert.Write()
784 self.hist_outbit_expert.Write()
785 self.h_Esum_ECL.Write()
786 self.h_E_ECL.Write()
787 self.h_theta_ECL.Write()
788 self.h_thetaID_ECL.Write()
789 self.h_phi_ECL.Write()
790 self.h_sector_BKLM.Write()
791 self.h_sector_EKLM.Write()
792 self.d_w.Write()
793 self.d_phi.Write()
794 self.d_z0_3d.Write()
795 self.d_z0_nn.Write()
796 self.d_E_ECL.Write()
797 self.tfile.Close()
798
799
800
801
802mypath = basf2.Path() # create a new path
803
804# add input data and ParticleLoader modules to the path
805ma.inputMdst("../TRGValidationGen.root", path=mypath)
806
807
808stdc.stdMu(listtype='all', path=mypath)
809stdc.stdE(listtype='all', path=mypath)
810
811ma.matchMCTruth('mu+:all', path=mypath)
812ma.matchMCTruth('e+:all', path=mypath)
813
814ma.cutAndCopyList('mu+:true', 'mu+:all', cut='isSignal>0 and mcPrimary>0', path=mypath)
815ma.cutAndCopyList('e+:true', 'e+:all', cut='isSignal>0 and mcPrimary>0', path=mypath)
816
817
818pvars = ["p", "theta", "phi"]
819dvars = vu.create_aliases(pvars, 'daughter(0,{variable})', "emu")
820
821bits = ['ty_0', 'ehigh', 'clst_0', 'klm_0', "klm_hit"]
822varlists = []
823for bit in bits:
824 inb_cut = f"L1Input({bit})>0"
825 ma.cutAndCopyList(f'mu+:{bit}', 'mu+:true', cut=inb_cut, path=mypath)
826 ma.cutAndCopyList(f'e+:{bit}', 'e+:true', cut=inb_cut, path=mypath)
827
828 ma.variablesToHistogram(decayString=f'mu+:{bit}',
829 variables=[('p', 50, 0.5, 3.0),
830 ('clusterE', 50, 0.5, 3.0),
831 ('theta', 40, 0.5, 2.5),
832 ('phi', 40, -3.141593, 3.141593)],
833 filename='TRGValidation.root',
834 directory=f'mu_{bit}',
835 path=mypath
836 )
837
838 ma.variablesToHistogram(decayString=f'e+:{bit}',
839 variables=[('p', 50, 0.5, 3.0),
840 ('clusterE', 50, 0.5, 3.0),
841 ('theta', 40, 0.5, 2.5),
842 ('phi', 40, -3.141593, 3.141593)],
843 filename='TRGValidation.root',
844 directory=f'e_{bit}',
845 path=mypath
846 )
847
848
849ma.variablesToHistogram(decayString='mu+:true',
850 variables=[('p', 50, 0.5, 3.0),
851 ('clusterE', 50, 0.5, 3.0),
852 ('theta', 40, 0.5, 2.5),
853 ('phi', 40, -3.141593, 3.141593)],
854 filename='TRGValidation.root',
855 directory='mu_all',
856 path=mypath
857 )
858
859ma.variablesToHistogram(decayString='e+:true',
860 variables=[('p', 50, 0.5, 3.0),
861 ('clusterE', 50, 0.5, 3.0),
862 ('theta', 40, 0.5, 2.5),
863 ('phi', 40, -3.141593, 3.141593)],
864 filename='TRGValidation.root',
865 directory='e_all',
866 path=mypath
867 )
868
869
870# process the data
871basf2.process(mypath)
872
873# INput
874main = basf2.create_path()
875root_input = basf2.register_module("RootInput")
876root_input.param("inputFileName", "../TRGValidationGen.root")
877main.add_module(root_input)
878main.add_module(MakePlots())
879main.add_module('Progress')
880basf2.process(main)
Belle2::PyDBObj
Class to access a DBObjPtr from Python.
Definition: PyDBObj.h:50
Belle2::PyStoreArray
A (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:72
Belle2::PyStoreObj
a (simplified) python wrapper for StoreObjPtr.
Definition: PyStoreObj.h:67
TRGValidation.MakePlots.h_Esum_ECL
h_Esum_ECL
validation histogram
Definition: TRGValidation.py:509
TRGValidation.MakePlots.h_sector_BKLM
h_sector_BKLM
validation histogram
Definition: TRGValidation.py:529
TRGValidation.MakePlots.d_w
d_w
validation histogram
Definition: TRGValidation.py:564
TRGValidation.MakePlots.hist_outbit
hist_outbit
validation histogram
Definition: TRGValidation.py:487
TRGValidation.MakePlots.h_E_ECL
h_E_ECL
validation histogram
Definition: TRGValidation.py:504
TRGValidation.MakePlots.set_descr_expert
def set_descr_expert(self, histogram, description, check)
Definition: TRGValidation.py:380
TRGValidation.MakePlots.hist_outbit_expert
hist_outbit_expert
validation histogram for experts
Definition: TRGValidation.py:547
TRGValidation.MakePlots.d_E_ECL
d_E_ECL
validation histogram
Definition: TRGValidation.py:615
TRGValidation.MakePlots.h_thetaID_ECL
h_thetaID_ECL
validation histogram
Definition: TRGValidation.py:519
TRGValidation.MakePlots.set_style
def set_style(self, histogram, xtitle, ytitle)
Definition: TRGValidation.py:397
TRGValidation.MakePlots.set_descr_shifter
def set_descr_shifter(self, histogram, description, check)
Definition: TRGValidation.py:363
TRGValidation.MakePlots.hist_inbit_expert
hist_inbit_expert
validation histogram for experts
Definition: TRGValidation.py:540
TRGValidation.MakePlots.d_w_2
d_w_2
validation histogram
Definition: TRGValidation.py:568
TRGValidation.MakePlots.d_phi_3
d_phi_3
validation histogram
Definition: TRGValidation.py:589
TRGValidation.MakePlots.hist_inbit
hist_inbit
validation histogram
Definition: TRGValidation.py:478
TRGValidation.MakePlots.d_z0_nn
d_z0_nn
validation histogram
Definition: TRGValidation.py:608
TRGValidation.MakePlots.h_theta_ECL
h_theta_ECL
validation histogram
Definition: TRGValidation.py:514
TRGValidation.MakePlots.d_phi_2
d_phi_2
validation histogram
Definition: TRGValidation.py:583
TRGValidation.MakePlots.get_relative
def get_relative(self, hist1, hist2, title, particle, trgbit)
Definition: TRGValidation.py:412
TRGValidation.MakePlots.h_sector_EKLM
h_sector_EKLM
validation histogram
Definition: TRGValidation.py:534
TRGValidation.MakePlots.Nevent
Nevent
number of events
Definition: TRGValidation.py:452
TRGValidation.MakePlots.h_phi_ECL
h_phi_ECL
validation histogram
Definition: TRGValidation.py:524
TRGValidation.MakePlots.d_z0_3d
d_z0_3d
validation histogram
Definition: TRGValidation.py:601
TRGValidation.MakePlots.d_phi
d_phi
validation histogram
Definition: TRGValidation.py:577