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