19 <input>TRGValidationGen.root</input>
20 <output>TRGValidation.root</output>
21 <contact>yinjh2012@korea.ac.kr</contact>
22 <description>makes validation plots for TRG</description>
28 import modularAnalysis
as ma
29 import stdCharged
as stdc
33 from ROOT
import Belle2, TH1F, TFile, TNamed, TEfficiency
34 from math
import pi
as 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"]
200 outputBits_expert = [
361 moniInbits_expert = []
362 moniOutbits_expert = []
367 Make validation histograms for trg ecl/cdc/klm
372 Sets description, check and contact to validation histogram.
373 :param h validation histogram
374 :param Descr description text
377 descr = TNamed(
"Description", description)
378 histogram.GetListOfFunctions().Add(descr)
379 Check = TNamed(
"Check", check)
380 histogram.GetListOfFunctions().Add(Check)
381 contact = TNamed(
"Contact",
"yinjh2012@korea.ac.kr")
382 histogram.GetListOfFunctions().Add(contact)
383 Meta = TNamed(
"MetaOptions",
"shifter")
384 histogram.GetListOfFunctions().Add(Meta)
388 Sets x-y titles, and sets histogram style.
389 :param xtitle X-axis title
390 :param xtitle Y-axis title
393 histogram.GetXaxis().SetTitle(xtitle)
394 histogram.GetYaxis().SetTitle(ytitle)
395 histogram.GetYaxis().SetTitleOffset(1.4)
396 histogram.GetXaxis().SetTitleSize(0.045)
397 histogram.GetYaxis().SetLabelSize(0.020)
398 histogram.SetLineColor(ROOT.kBlack)
402 Get relative ratio between two hists.
403 :param hist1 numerator
404 :param hist2 denominator
405 :param title new histogram title
408 bin1 = hist1.GetXaxis().GetNbins()
409 bin2 = hist2.GetXaxis().GetNbins()
410 Xmin1 = hist1.GetXaxis().GetXmin()
411 Xmax1 = hist1.GetXaxis().GetXmax()
412 Xmin2 = hist2.GetXaxis().GetXmin()
413 Xmax2 = hist2.GetXaxis().GetXmax()
414 if bin1 != bin2
or Xmin1 != Xmin2
or Xmax1 != Xmax2:
415 print(
"warning: two histograms cannot be compared!!")
418 teff = TEfficiency(hist1, hist2)
419 htitle = f
'h_eff_{title}_{particle}_{trgbit}'
420 heff = TH1F(htitle, f
"efficiency histogram of {particle} for {trgbit}", bin1, Xmin1, Xmax1)
421 for ibin
in range(1, bin1+1):
422 binc0 = teff.GetEfficiency(ibin)
423 bine0 = (teff.GetEfficiencyErrorUp(ibin) + teff.GetEfficiencyErrorLow(ibin)) / 2.0
424 heff.SetBinContent(ibin, binc0)
425 heff.SetBinError(ibin, bine0)
426 heff.GetYaxis().SetRangeUser(0.0, 1.2)
429 def initialize(self):
432 self.
tfiletfile = TFile(
"./TRGValidation.root",
"update")
438 inbitname = m_dbinput.getinbitname(i)
439 outbitname = m_dbftdl.getoutbitname(i)
440 if inbitname
in inputBits:
441 moniInbits.append({
"name": inbitname,
"index": i})
442 if outbitname
in outputBits:
443 moniOutbits.append({
"name": outbitname,
"index": i})
444 if inbitname
in inputBits_expert:
445 moniInbits_expert.append({
"name": inbitname,
"index": i})
446 if outbitname
in outputBits_expert:
447 moniOutbits_expert.append({
"name": outbitname,
"index": i})
449 mc_px =
"MCParticles.m_momentum_x"
450 mc_py =
"MCParticles.m_momentum_y"
451 trk2d_omega =
"TRGCDC2DFinderTracks.m_omega"
452 trk2d_phi =
"TRGCDC2DFinderTracks.m_phi0"
454 ROOT.gStyle.SetOptStat(ROOT.kFALSE)
457 n_inbit_test = len(inputBits)
458 self.
hist_inbithist_inbit = TH1F(
"hin",
"trigger input bits", n_inbit_test, 0, n_inbit_test)
459 self.
hist_inbithist_inbit.GetXaxis().SetRangeUser(0, n_inbit_test)
460 self.
hist_inbithist_inbit.GetXaxis().SetLabelSize(0.04)
465 n_oubit_test = len(outputBits)
466 self.
hist_outbithist_outbit = TH1F(
"hout",
"trigger output bits", n_oubit_test, 0, n_oubit_test)
467 self.
hist_outbithist_outbit.GetXaxis().SetRangeUser(0, n_oubit_test)
468 self.
hist_outbithist_outbit.GetXaxis().SetLabelSize(0.04)
472 for i
in range(n_inbit_test):
473 self.
hist_inbithist_inbit.GetXaxis().SetBinLabel(
474 self.
hist_inbithist_inbit.GetXaxis().FindBin(i + 0.5), inputBits[i])
475 for i
in range(n_oubit_test):
476 self.
hist_outbithist_outbit.GetXaxis().SetBinLabel(
477 self.
hist_outbithist_outbit.GetXaxis().FindBin(i + 0.5), outputBits[i])
480 self.
h_E_ECLh_E_ECL = TH1F(
"h_E_ECL",
"ECL cluster energy [50 MeV]", 200, 0, 10)
482 self.
set_styleset_style(self.
h_E_ECLh_E_ECL,
"ECL cluster energy (GeV)",
"Events/(50 MeV)")
485 self.
h_Esum_ECLh_Esum_ECL = TH1F(
"h_Esum_ECL",
"sum of ECL cluster energy [50 MeV]", 200, 0, 10)
487 self.
set_styleset_style(self.
h_Esum_ECLh_Esum_ECL,
"sum of ECL cluster energy (GeV)",
"Events/(50 MeV)")
490 self.
h_theta_ECLh_theta_ECL = TH1F(
"h_theta_ECL",
"TRG ECL cluster #theta [1.4 degrees]", 128, 0, 180)
491 self.
set_descrset_descr(self.
h_theta_ECLh_theta_ECL,
"TRG ECL cluster polar angle",
"unform in barrel")
492 self.
set_styleset_style(self.
h_theta_ECLh_theta_ECL,
"TRG ECL cluster polar angle (degree)",
"Events/(1.4 degree)")
495 self.
h_thetaID_ECLh_thetaID_ECL = TH1F(
"h_thetaID_ECL",
"ECL cluster TC ID", 610, 0, 610)
500 self.
h_phi_ECLh_phi_ECL = TH1F(
"h_phi_ECL",
"TRG ECL cluster phi [2.0 degrees]", 180, -180, 180)
501 self.
set_descrset_descr(self.
h_phi_ECLh_phi_ECL,
"TRG ECL cluster phi distribution",
"distribute unformly")
502 self.
set_styleset_style(self.
h_phi_ECLh_phi_ECL,
"ECL cluster #phi (degree) ",
"Events/(2.0 degrees)")
505 self.
h_sector_BKLMh_sector_BKLM = TH1F(
"h_sector_BKLM",
"BKLM TRG hit sector", 10, 0, 10)
510 self.
h_sector_EKLMh_sector_EKLM = TH1F(
"h_sector_EKLM",
"EKLM TRG hit sector", 10, 0, 10)
515 n_inbit_expert = len(inputBits_expert)
516 self.
hist_inbit_experthist_inbit_expert = TH1F(
"hin_expert",
"trigger input bits", n_inbit_expert, 0, n_inbit_expert)
517 self.
hist_inbit_experthist_inbit_expert.GetXaxis().SetRangeUser(0, n_inbit_expert)
521 n_oubit_expert = len(outputBits_expert)
522 self.
hist_outbit_experthist_outbit_expert = TH1F(
"hout_expert",
"trigger output bits", n_oubit_expert, 0, n_oubit_expert)
526 for i
in range(n_inbit_expert):
528 self.
hist_inbit_experthist_inbit_expert.GetXaxis().FindBin(i + 0.5), inputBits_expert[i])
529 for i
in range(n_oubit_expert):
531 self.
hist_outbit_experthist_outbit_expert.GetXaxis().FindBin(i + 0.5), outputBits_expert[i])
533 mc =
"abs(MCParticles.m_pdg)==11&&MCParticles.m_status==11"
534 tree = ROOT.TChain(
"tree")
535 tree.Add(
"../TRGValidationGen.root")
538 self.
d_wd_w = TH1F(
"d_w",
"#Deltaw of CDC 2D finder, w = 0.00449/p_{t}", 50, -0.02, 0.02)
539 self.
d_w_2d_w_2 = TH1F(
"d_w_2",
"d_w_2", 50, -0.02, 0.02)
540 tree.Draw(
"({0} - 0.00449)/sqrt({1}*{1} + {2}*{2})>> d_w".format(trk2d_omega, mc_px, mc_py),
541 "MCParticles.m_pdg<0&&" + mc)
542 tree.Draw(
"({0} + 0.00449)/sqrt({1}*{1} + {2}*{2}) >> d_w_2".format(trk2d_omega, mc_px, mc_py),
543 "MCParticles.m_pdg>0 && " + mc)
545 self.
set_descrset_descr(self.
d_wd_w,
"Comparison on w (=0.00449/pt) of a track between CDC 2D finder output and MC.",
546 "A clear peak at 0 with tail.")
547 self.
set_styleset_style(self.
d_wd_w,
"#Deltaw",
"Events/(0.08)")
550 self.
d_phid_phi = TH1F(
"d_phi",
"#Delta#phi of CDC 2D finder", 50, -0.5, 0.5)
551 tree.Draw(
"{0}-atan({1}/{2})>>d_phi".format(trk2d_phi, mc_py, mc_px),
552 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
553 "fabs({0}-atan({1}/{2}))<{3}&&".format(trk2d_phi, mc_py, mc_px, PI) + mc)
555 self.
d_phi_2d_phi_2 = TH1F(
"d_phi_2",
"d_phi_2", 50, -0.5, 0.5)
556 tree.Draw(
"{0}-atan({1}/{2})-{3}>>d_phi_2".format(trk2d_phi, mc_py, mc_px, PI),
557 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
558 "{0}-atan({1}/{2})>={3} &&".format(trk2d_phi, mc_py, mc_px, PI) + mc)
560 self.
d_phi_3d_phi_3 = TH1F(
"d_phi_3",
"d_phi_3", 50, -0.5, 0.5)
561 tree.Draw(
"{0}-atan({1}/{2})+{3}>>d_phi_2".format(trk2d_phi, mc_py, mc_px, PI),
562 "MCParticles.m_status==11&&abs(MCParticles.m_pdg)==11",
563 "{0}-atan({1}/{2})<=-{3}&&".format(trk2d_phi, mc_py, mc_px, PI) + mc)
567 self.
set_descrset_descr(self.
d_phid_phi,
"Comparison on phi_i of a track between CDC 2D finder output and MC.",
568 "A Gaussian peak at 0.")
569 self.
set_styleset_style(self.
d_phid_phi,
"#Delta#phi [rad]",
"Events/(0.02 rad)")
572 self.
d_z0_3dd_z0_3d = TH1F(
"d_z0_3d",
"#Deltaz0 of CDC 3D fitter", 60, -30, 30)
573 tree.Draw(
"TRGCDC3DFitterTracks.m_z0-MCParticles.m_productionVertex_z>>d_z0_3d", mc)
574 self.
set_descrset_descr(self.
d_z0_3dd_z0_3d,
"Comparison on z0 of a track between CDC 2D fitter output and MC.",
575 "A Gaussian peak at 0 with small tail.")
576 self.
set_styleset_style(self.
d_z0_3dd_z0_3d,
"#Deltaz0 [cm]",
"Events/(1 cm)")
579 self.
d_z0_nnd_z0_nn = TH1F(
"d_z0_nn",
"#Deltaz0 of CDC Neuro", 60, -30, 30)
580 tree.Draw(
"TRGCDCNeuroTracks.m_z0-MCParticles.m_productionVertex_z>>d_z0_nn", mc)
581 self.
set_descrset_descr(self.
d_z0_nnd_z0_nn,
"Comparison on z0 of a track between CDC Neuro output and MC.",
582 "A Gaussian peak at 0 with small tail.")
583 self.
set_styleset_style(self.
d_z0_nnd_z0_nn,
"#Deltaz0 [cm]",
"Events/(1 cm)")
586 self.
d_E_ECLd_E_ECL = TH1F(
"d_E_ECL",
"#DeltaE of ECL clustering", 100, -6, 6)
587 tree.Draw(
"TRGECLClusters.m_edep-MCParticles.m_energy>>d_E_ECL", mc)
588 self.
set_descrset_descr(self.
d_E_ECLd_E_ECL,
"Comparison on deposit energy between ECL cluster output and MC.",
589 "A peak around -0.5 ~ 0 with a tail toward -6.")
590 self.
set_styleset_style(self.
d_E_ECLd_E_ECL,
"#DeltaE [GeV]",
"Events/(0.12 GeV)")
609 for bit
in moniInbits:
610 binIndex = inputBits.index(bit[
"name"])
611 bitIndex = bit[
"index"]
612 if trg_summary.testInput(bitIndex):
613 self.
hist_inbithist_inbit.Fill(binIndex + 0.5)
614 for bit
in moniOutbits:
615 binIndex = outputBits.index(bit[
"name"])
616 bitIndex = bit[
"index"]
617 if trg_summary.testFtdl(bitIndex):
620 for bit
in moniInbits_expert:
621 binIndex = inputBits_expert.index(bit[
"name"])
622 bitIndex = bit[
"index"]
623 if trg_summary.testInput(bitIndex):
625 for bit
in moniOutbits_expert:
626 binIndex = outputBits_expert.index(bit[
"name"])
627 bitIndex = bit[
"index"]
628 if trg_summary.testFtdl(bitIndex):
632 for cluster
in clusters:
633 x = cluster.getPositionX()
634 y = cluster.getPositionY()
635 z = cluster.getPositionZ()
636 e = cluster.getEnergyDep()
639 vec = ROOT.TVector3(x, y, z)
640 self.
h_theta_ECLh_theta_ECL.Fill(vec.Theta() * Fac)
642 self.
h_phi_ECLh_phi_ECL.Fill(vec.Phi() * Fac)
647 self.
h_sector_BKLMh_sector_BKLM.Fill(klmSummary.getBKLM_n_trg_sectors())
648 self.
h_sector_EKLMh_sector_EKLM.Fill(klmSummary.getEKLM_n_trg_sectors())
657 bits = [
'ty_0',
'ehigh',
'clst_0',
'klm_0']
658 for part
in [
"e",
"mu"]:
659 h_gen_p = self.
tfiletfile.Get(f
"{part}_all/p")
660 h_gen_E = self.
tfiletfile.Get(f
"{part}_all/clusterE")
661 h_gen_theta = self.
tfiletfile.Get(f
"{part}_all/theta")
662 h_gen_phi = self.
tfiletfile.Get(f
"{part}_all/phi")
665 h_p = self.
tfiletfile.Get(f
"{part}_{bit}/p")
666 h_E = self.
tfiletfile.Get(f
"{part}_{bit}/clusterE")
667 h_theta = self.
tfiletfile.Get(f
"{part}_{bit}/theta")
668 h_phi = self.
tfiletfile.Get(f
"{part}_{bit}/phi")
670 h_eff_p = self.
get_relativeget_relative(h_p, h_gen_p,
'p', part, bit)
671 self.
set_styleset_style(h_eff_p,
"Momentum",
"Efficiency")
672 if (bit ==
'ty_0' or bit ==
'klm_0')
and part ==
'mu':
673 self.
set_descrset_descr(h_eff_p,
"Momentum dependent efficiency for shifters",
"")
675 h_eff_E = self.
get_relativeget_relative(h_E, h_gen_E,
'E', part, bit)
676 self.
set_styleset_style(h_eff_E,
"Energy",
"Efficiency")
677 if (bit ==
'ehigh' or bit ==
'clst_0')
and part ==
'e':
678 self.
set_descrset_descr(h_eff_E,
"Energy dependent efficiency for shifters",
"")
680 h_eff_theta = self.
get_relativeget_relative(h_theta, h_gen_theta,
'theta', part, bit)
681 self.
set_styleset_style(h_eff_theta,
"Polar angle",
"Efficiency")
683 h_eff_phi = self.
get_relativeget_relative(h_phi, h_gen_phi,
'phi', part, bit)
684 self.
set_styleset_style(h_eff_phi,
"phi angle",
"Efficiency")
708 self.
d_phid_phi.Write()
712 self.
tfiletfile.Close()
717 mypath = basf2.Path()
720 ma.inputMdst(
"../TRGValidationGen.root", path=mypath)
723 stdc.stdMu(listtype=
'all', path=mypath)
724 stdc.stdE(listtype=
'all', path=mypath)
726 ma.matchMCTruth(
'mu+:all', path=mypath)
727 ma.matchMCTruth(
'e+:all', path=mypath)
729 ma.cutAndCopyList(
'mu+:true',
'mu+:all', cut=
'isSignal>0 and mcPrimary>0', path=mypath)
730 ma.cutAndCopyList(
'e+:true',
'e+:all', cut=
'isSignal>0 and mcPrimary>0', path=mypath)
733 pvars = [
"p",
"theta",
"phi"]
734 dvars = vu.create_aliases(pvars,
'daughter(0,{variable})',
"emu")
736 bits = [
'ty_0',
'ehigh',
'clst_0',
'klm_0']
739 inb_cut = f
"L1Input({bit})>0"
740 ma.cutAndCopyList(f
'mu+:{bit}',
'mu+:true', cut=inb_cut, path=mypath)
741 ma.cutAndCopyList(f
'e+:{bit}',
'e+:true', cut=inb_cut, path=mypath)
743 ma.variablesToHistogram(decayString=f
'mu+:{bit}',
744 variables=[(
'p', 40, 1.0, 3.0),
745 (
'clusterE', 40, 1.0, 3.0),
746 (
'theta', 40, 0.5, 2.5),
747 (
'phi', 40, -3.141593, 3.141593)],
748 filename=
'TRGValidation.root',
749 directory=f
'mu_{bit}',
753 ma.variablesToHistogram(decayString=f
'e+:{bit}',
754 variables=[(
'p', 40, 1.0, 3.0),
755 (
'clusterE', 40, 1.0, 3.0),
756 (
'theta', 40, 0.5, 2.5),
757 (
'phi', 40, -3.141593, 3.141593)],
758 filename=
'TRGValidation.root',
759 directory=f
'e_{bit}',
764 ma.variablesToHistogram(decayString=
'mu+:true',
765 variables=[(
'p', 40, 1.0, 3.0),
766 (
'clusterE', 40, 1.0, 3.0),
767 (
'theta', 40, 0.5, 2.5),
768 (
'phi', 40, -3.141593, 3.141593)],
769 filename=
'TRGValidation.root',
774 ma.variablesToHistogram(decayString=
'e+:true',
775 variables=[(
'p', 40, 1.0, 3.0),
776 (
'clusterE', 40, 1.0, 3.0),
777 (
'theta', 40, 0.5, 2.5),
778 (
'phi', 40, -3.141593, 3.141593)],
779 filename=
'TRGValidation.root',
786 basf2.process(mypath)
789 main = basf2.create_path()
790 root_input = basf2.register_module(
"RootInput")
791 root_input.param(
"inputFileName",
"../TRGValidationGen.root")
792 main.add_module(root_input)
794 main.add_module(
'Progress')
Class to access a DBObjPtr from Python.
A (simplified) python wrapper for StoreArray.
a (simplified) python wrapper for StoreObjPtr.
h_Esum_ECL
validation histogram
h_sector_BKLM
validation histogram
hist_outbit
validation histogram
h_E_ECL
validation histogram
hist_outbit_expert
validation histogram for experts
d_E_ECL
validation histogram
h_thetaID_ECL
validation histogram
def set_style(self, histogram, xtitle, ytitle)
hist_inbit_expert
validation histogram for experts
hist_inbit
validation histogram
d_z0_nn
validation histogram
h_theta_ECL
validation histogram
def set_descr(self, histogram, description, check)
def get_relative(self, hist1, hist2, title, particle, trgbit)
h_sector_EKLM
validation histogram
h_phi_ECL
validation histogram
d_z0_3d
validation histogram
d_phi
validation histogram