File indexing completed on 2025-08-27 22:40:50
0001
0002
0003 import os
0004 import argparse
0005 import numpy as np
0006 import hist
0007 import matplotlib.pyplot as plt
0008 from matplotlib.transforms import blended_transform_factory
0009 import array
0010 import ROOT
0011 import mplhep as hep
0012 hep.style.use("CMS")
0013
0014 import warnings
0015 warnings.filterwarnings("ignore", message="The value of the smallest subnormal")
0016
0017 class dotdict(dict):
0018 """dot.notation access to dictionary attributes"""
0019 __getattr__ = dict.get
0020 __setattr__ = dict.__setitem__
0021 __delattr__ = dict.__delitem__
0022
0023 def CheckRootFile(hname, rebin=None):
0024 hist_orig = file.Get(hname)
0025 if not hist_orig:
0026 raise RuntimeError(f"WARNING: Histogram {hname} not found.")
0027
0028 hist = hist_orig.Clone(hname + "_clone")
0029 hist.SetDirectory(0)
0030
0031 if rebin is not None:
0032 if isinstance(rebin, (int, float)):
0033 hist = hist.Rebin(int(rebin), hname + "_rebin")
0034 elif hasattr(rebin, '__iter__'):
0035 bin_edges_c = array.array('d', rebin)
0036 hist = hist.Rebin(len(bin_edges_c) - 1, hname + "_rebin", bin_edges_c)
0037 else:
0038 raise ValueError(f"Unknown type for rebin: {type(rebin)}")
0039
0040 return hist
0041
0042 def findBestGaussianCoreFit(histo, quiet=True):
0043 mean, rms = histo.GetMean(), histo.GetRMS()
0044
0045 HalfMaxBinLow = histo.FindFirstBinAbove(histo.GetMaximum()/2)
0046 HalfMaxBinHigh = histo.FindLastBinAbove(histo.GetMaximum()/2)
0047 WidthAtHalfMaximum = 0.5*(histo.GetBinCenter(HalfMaxBinHigh) - histo.GetBinCenter(HalfMaxBinLow))
0048 Xmax = histo.GetXaxis().GetBinCenter(histo.GetMaximumBin())
0049
0050 gausTF1 = ROOT.TF1()
0051
0052 Pvalue = 0.
0053 RangeLow = histo.GetBinLowEdge(2)
0054 RangeUp = histo.GetBinLowEdge(histo.GetNbinsX())
0055
0056 PvalueBest = 0.
0057 RangeLowBest = 0.
0058 RangeUpBest = 0.
0059
0060 if Xmax < mean:
0061 meanForRange = Xmax
0062 else:
0063 meanForRange = mean
0064
0065 if WidthAtHalfMaximum < rms and WidthAtHalfMaximum>0:
0066 spreadForRange = WidthAtHalfMaximum
0067 else:
0068 spreadForRange = rms
0069
0070 rms_step_minus = 2.2
0071 while rms_step_minus>1.1:
0072 RangeLow = meanForRange - rms_step_minus*spreadForRange
0073 rms_step_plus = rms_step_minus
0074
0075 while rms_step_plus>0.7:
0076 RangeUp = meanForRange + rms_step_plus*spreadForRange
0077 if quiet:
0078 histo.Fit("gaus", "0Q", "0", RangeLow, RangeUp)
0079 else:
0080 histo.Fit("gaus", "0", "0", RangeLow, RangeUp)
0081 gausTF1 = histo.GetListOfFunctions().FindObject("gaus")
0082 ChiSquare = gausTF1.GetChisquare()
0083 ndf = gausTF1.GetNDF()
0084 Pvalue = ROOT.TMath.Prob(ChiSquare, ndf)
0085
0086 if Pvalue > PvalueBest:
0087 PvalueBest = Pvalue
0088 RangeLowBest = RangeLow
0089 RangeUpBest = RangeUp
0090 ndfBest = ndf
0091 ChiSquareBest = ChiSquare
0092 StepMinusBest = rms_step_minus
0093 StepPlusBest = rms_step_plus
0094 meanForRange = gausTF1.GetParameter(1)
0095
0096 if not quiet:
0097 print("\n\nFitting range used: [Mean - " + str(rms_step_minus) + " sigma, Mean + " + str(rms_step_plus) + " sigma ] ")
0098 print("ChiSquare = " + str(ChiSquare) + " NDF = " + str(ndf) + " Prob = " + str(Pvalue) + " Best Prob so far = " + str(PvalueBest))
0099
0100 rms_step_plus = rms_step_plus - 0.1
0101
0102 rms_step_minus = rms_step_minus - 0.1
0103
0104 if quiet:
0105 histo.Fit("gaus", "0Q", "0", RangeLowBest, RangeUpBest)
0106 else:
0107 histo.Fit("gaus","0","0",RangeLowBest, RangeUpBest)
0108 print("\n\n\nMean = " + str(mean) + " Xmax = " + str(Xmax) + " RMS = " + str(rms) + " WidthAtHalfMax = " + str(WidthAtHalfMaximum))
0109 print("Fit found!")
0110 print("Final fitting range used: [Mean(Xmax) - " + str(StepMinusBest) + " rms(WHF), Mean(Xmax) + " + str(StepPlusBest) + " rms(WHF) ] ")
0111 print("ChiSquare = " + str(ChiSquareBest) + " NDF = " + str(ndfBest) + " Prob = " + str(PvalueBest) + "\n\n")
0112
0113 return histo.GetListOfFunctions().FindObject("gaus")
0114
0115 def tails_errors(n1, n2):
0116 """Compute and return the lower and upper errors."""
0117 if n1 == 0:
0118 return 0, 0
0119 elif n1 < 0 or n2 < 0:
0120 raise RuntimeError(f"Negative number of entries! n1={n1}, n2={n2}")
0121 elif n1 < n2:
0122 raise RuntimeError(f"n1 is smaller than n2! n1={n1}, n2={n2}")
0123 else:
0124 return ( n2/n1 - (n2/n1 + 0.5/n1 - np.sqrt(n2/pow(n1,2)*(1-n2/n1) + 0.25/pow(n1,2))) / (1+1.0/n1),
0125 (n2/n1 + 0.5/n1 + np.sqrt(n2/pow(n1,2)*(1-n2/n1) + 0.25/pow(n1,2))) / (1+1.0/n1) - n2/n1 )
0126
0127 def rate_errorbar_declutter(eff, err, yaxmin, frac=0.01):
0128 """
0129 Filter uncertainties if they lie below the minimum (vertical) axis value.
0130 Used to plot the filtered points differently, for instance by displaying only the upper uncertainty.
0131 Uncertainties are trimmed to 1. or 0. if these values are crossed
0132 (the retrieved Clopper-Pearson uncertainties have been symmetrized in the DQMGenericClient).
0133 """
0134 filt = eff-err/2 <= yaxmin
0135 eff_filt = np.where(filt, np.nan, eff)
0136 err_filt = np.where(filt, np.nan, err)
0137
0138 up_error = eff_filt + err_filt/2
0139 transform = blended_transform_factory(plotter.ax.transData, plotter.ax.transAxes)
0140 up_error = np.where(np.isnan(up_error), frac, up_error)
0141 up_error = np.where(up_error != frac, np.nan, up_error)
0142
0143 filt_limit_one = eff_filt+err_filt/2 > 1.
0144 filt_limit_zero = eff_filt-err_filt/2 < 0.
0145 err_hi = np.where(filt_limit_one, 0., err_filt)
0146 err_lo = np.where(filt_limit_zero, 0., err_filt)
0147 return eff_filt, (err_lo/2,err_hi/2), up_error, transform
0148
0149 def define_bins(h):
0150 """
0151 Computes the number of bins, edges, centers and widths of a histogram.
0152 """
0153 N = h.GetNbinsX()
0154 edges = np.array([h.GetBinLowEdge(i+1) for i in range(N)])
0155 edges = np.append(edges, h.GetBinLowEdge(N+1))
0156 return N, edges, 0.5*(edges[:-1]+edges[1:]), np.diff(edges)
0157
0158 def define_bins_2D(h):
0159 Nx = h.GetNbinsX()
0160 Ny = h.GetNbinsY()
0161
0162 x_edges = np.array([h.GetXaxis().GetBinLowEdge(i+1) for i in range(Nx)])
0163 x_edges = np.append(x_edges, h.GetXaxis().GetBinUpEdge(Nx))
0164
0165 y_edges = np.array([h.GetYaxis().GetBinLowEdge(j+1) for j in range(Ny)])
0166 y_edges = np.append(y_edges, h.GetYaxis().GetBinUpEdge(Ny))
0167
0168 return Nx, Ny, x_edges, y_edges
0169
0170 def histo_values_errors(h):
0171 N = h.GetNbinsX()
0172 values = np.array([h.GetBinContent(i+1) for i in range(N)])
0173 errors = np.array([h.GetBinError(i+1) for i in range(N)])
0174 return values, errors
0175
0176 def histo_values_2D(h, error=False):
0177 Nx = h.GetNbinsX()
0178 Ny = h.GetNbinsY()
0179 values = np.array([
0180 [h.GetBinContent(i+1, j+1) for i in range(Nx)]
0181 for j in range(Ny)
0182 ])
0183 return values
0184
0185 class Plotter:
0186 def __init__(self, label, fontsize=18, grid_color='grey'):
0187 self._fig, self._ax = plt.subplots(figsize=(10, 10))
0188 self.fontsize = fontsize
0189
0190 hep.cms.text(' Phase-2 Simulation Preliminary', ax=self._ax, fontsize=fontsize)
0191 hep.cms.lumitext(label + " | 14 TeV", ax=self._ax, fontsize=fontsize)
0192 if grid_color:
0193 self._ax.grid(which='major', color=grid_color)
0194
0195 self.extensions = ('png', 'pdf')
0196
0197 @property
0198 def fig(self):
0199 return self._fig
0200
0201 @property
0202 def ax(self):
0203 return self._ax
0204
0205 def labels(self, x, y, legend_title=None, legend_loc='upper right'):
0206 self._ax.set_xlabel(x)
0207 self._ax.set_ylabel(y)
0208 if legend_title is not None:
0209 self._ax.legend(title=legend_title,
0210 title_fontsize=self.fontsize, fontsize=self.fontsize,
0211 loc=legend_loc)
0212
0213 def limits(self, x=None, y=None, logY=False, logX=False):
0214 if x:
0215 self._ax.set_xlim(x)
0216 if y:
0217 self._ax.set_ylim(y)
0218 if logX:
0219 self._ax.set_xscale('log')
0220 if logY:
0221 self._ax.set_yscale('log')
0222
0223 def save(self, name):
0224 for ext in self.extensions:
0225 print(" ### INFO: Saving " + name + '.' + ext)
0226 plt.savefig(name + '.' + ext)
0227 plt.close()
0228
0229 class HLabels:
0230 _eff_types = ('Efficiency', 'Fake Rate', 'Gen Duplicates Rate', 'Reco Duplicates Rate')
0231 _resol_types = ('RecoOverGen', 'CorrOverGen', 'CorrOverReco')
0232
0233 def __init__(self, atype):
0234 assert atype in self._eff_types or atype in self._resol_types, f"Invalid type: {atype}"
0235 self.mytype = atype
0236
0237 @staticmethod
0238 def eff_types():
0239 return HLabels._eff_types
0240
0241 @staticmethod
0242 def resol_types():
0243 return HLabels._resol_types
0244
0245 @staticmethod
0246 def fraction_label(labtype):
0247 return {
0248 'neHad': 'Neutral Hadron Energy Fraction',
0249 'chEm': 'Charged EM Energy Fraction',
0250 'chHad': 'Charged Hadron Energy Fraction',
0251 'neEm': 'Neutral EM Energy Fraction',
0252 'nCost': '# Jet Constituents',
0253 }[labtype]
0254
0255 @staticmethod
0256 def multiplicity_label(labtype):
0257 return {
0258 'neHad': 'Neutral Hadron Multiplicity',
0259 'chEm': 'Charged EM Multiplicity',
0260 'chHad': 'Charged Hadron Multiplicity',
0261 'neEm': 'Neutral EM Multiplicity',
0262 'chMult': 'Charged Multiplicity',
0263 'neMult': 'Neutral Multiplicity',
0264 }[labtype]
0265
0266 @staticmethod
0267 def pt_label(labtype, average=False):
0268 labels = {
0269 "pt" : "$p_{T}\,$ [GeV]",
0270 "gen" : "$p_{T}^{gen}\,$ [GeV]",
0271 "reco" : "$p_{T}^{reco}\,$ [GeV]",
0272 "pt/gen" : "$p_{T}/p_{T}^{gen}\,$",
0273 "reco/gen" : "$p_{T}^{reco}/p_{T}^{gen}\,$",
0274 "corr/gen" : "$p_{T}^{corr}/p_{T}^{gen}\,$",
0275 "corr/reco" : "$p_{T}^{corr}/p_{T}^{reco}\,$",
0276 }
0277 if average:
0278 return f'<' + labels[labtype] + '>'
0279 else:
0280 return labels[labtype]
0281
0282 @staticmethod
0283 def resol_label(labtype):
0284 return f"$\sigma$({HLabels.pt_label(labtype)}) / {HLabels.pt_label(labtype, average=True)}"
0285
0286 def ytitle(self, resol=False):
0287 return {'Efficiency' : 'Jet Efficiency',
0288 'Fake Rate' : 'Jet Fake Rate',
0289 'Gen Duplicates Rate' : 'Jet Gen Duplicate Rate',
0290 'Reco Duplicates Rate' : 'Jet Reco Duplicate Rate',
0291 'RecoOverGen' : HLabels.resol_label('reco/gen') if resol else HLabels.pt_label('reco/gen', average=True),
0292 'CorrOverGen' : HLabels.resol_label('corr/gen') if resol else HLabels.pt_label('corr/gen', average=True),
0293 'CorrOverReco' : HLabels.resol_label('corr/reco') if resol else HLabels.pt_label('corr/reco', average=True),
0294 }[self.mytype]
0295
0296 @property
0297 def savename(self):
0298 return self.ytitle().replace(' ', '_')
0299
0300 @property
0301 def xvars(self):
0302 return ('Eta', 'Phi', 'Pt')
0303
0304 def nhisto(self, var, basedir):
0305 """Numerator"""
0306 return {'Efficiency' : basedir + '/MatchedGen' + var,
0307 'Fake Rate' : basedir + '/MatchedJet' + var,
0308 'Gen Duplicates Rate' : basedir + '/DuplicatesGen' + var,
0309 'Reco Duplicates Rate' : basedir + '/DuplicatesJet' + var
0310 }[self.mytype]
0311
0312 def dhisto(self, var, basedir):
0313 """Denominator"""
0314 return {'Efficiency' : basedir + '/Gen' + var,
0315 'Fake Rate' : basedir + '/Jet' + var,
0316 'Gen Duplicates Rate' : basedir + '/Gen' + var,
0317 'Reco Duplicates Rate' : basedir + '/Jet' + var
0318 }[self.mytype]
0319
0320 def rhisto(self, var, basedir):
0321 """Ratio"""
0322 return {'Efficiency': basedir + '/Eff_vs_' + var,
0323 'Fake Rate': basedir + '/Fake_vs_' + var,
0324 'Gen Duplicates Rate': basedir + '/DupGen_vs_' + var,
0325 'Reco Duplicates Rate': basedir + '/Dup_vs_' + var
0326 }[self.mytype]
0327
0328 def mhisto(self, var, basedir):
0329 if self.mytype in self._resol_types:
0330 var_type = 'Gen' if 'Gen' in self.mytype else ''
0331 name = basedir + f'/Res_{self.mytype}_{var_type}{var}_Mean'
0332 else:
0333 raise ValueError(f"Unavailable HLabels.mhisto for {self.mytype}")
0334 return name
0335
0336 def shisto(self, var, basedir):
0337 if self.mytype in self._resol_types:
0338 var_type = 'Gen' if 'Gen' in self.mytype else ''
0339 name = basedir + f'/Res_{self.mytype}_{var_type}{var}_Sigma'
0340 else:
0341 raise ValueError(f"Unavailable HLabels.shisto for {self.mytype}")
0342 return name
0343
0344 def leglabel(self, isNum):
0345 HowMany = ''
0346 HowMany = '' if 'Duplicates' not in self.mytype else ' multiple'
0347 if self.mytype in ('Efficiency', 'Gen Duplicates Rate'):
0348 if isNum:
0349 label = f'Gen Jets $p_T > 20$ GeV matched to{HowMany} HLT jets'
0350 else:
0351 label = 'Gen Jets $p_T > 20$ GeV'
0352 elif self.mytype in ('Fake Rate', 'Reco Duplicates Rate'):
0353 if isNum:
0354 label = f'HLT Jets $p_T > 30$ GeV matched to{HowMany} gen jets'
0355 else:
0356 label = 'HLT Jets $p_T > 30$ GeV'
0357 return label
0358
0359 class EtaInfo:
0360 """
0361 Manage information related with the eta region.
0362 """
0363 info = {
0364 'B': ('black', 'o', r'$|\eta|<1.5$'),
0365 'E': ('red', 's', r'$1.5\leq|\eta|<3$'),
0366 'F': ('blue', 'd', r'$3\leq|\eta|<6$')
0367 }
0368
0369 @staticmethod
0370 def color(x):
0371 return EtaInfo.info.get(x)[0]
0372
0373 @staticmethod
0374 def marker(x):
0375 return EtaInfo.info.get(x)[1]
0376
0377 @staticmethod
0378 def label(x):
0379 return EtaInfo.info.get(x)[2]
0380
0381 @staticmethod
0382 def regions():
0383 return ('B', 'E', 'F')
0384
0385
0386 if __name__ == '__main__':
0387 parser = argparse.ArgumentParser(description='Make HLT Jet validation plots.')
0388 parser.add_argument('-f', '--file', type=str, required=True, help='Paths to the DQM ROOT file.')
0389 parser.add_argument('-j', '--jet', type=str, default='hltAK4PFPuppiJets', help='Name of the jet collection')
0390 parser.add_argument('-o', '--odir', type=str, default="HLTJetValidationPlots", required=False, help='Path to the output directory.')
0391 parser.add_argument('-l', '--sample_label', type=str, default="QCD (200 PU)", required=False, help='Path to the output directory.')
0392 args = parser.parse_args()
0393
0394 if not os.path.exists(args.odir):
0395 os.makedirs(args.odir)
0396 histo2d_dir = os.path.join(args.odir, 'histo2D')
0397
0398 file = ROOT.TFile.Open(args.file)
0399 dqm_dir = f"DQMData/Run 1/HLT/Run summary/JetMET/JetValidation/{args.jet}"
0400 if not file.Get(dqm_dir):
0401 raise RuntimeError(f"Directory '{dqm_dir}' not found in {args.file}")
0402
0403 fontsize = 16
0404
0405 tprofile_rebinning = {'B': (30, 40, 50, 80, 100, 120, 140, 160, 200, 250, 300, 350, 400, 500, 600),
0406 'E': (30, 40, 50, 80, 100, 120, 140, 160, 200, 250, 300, 350, 400, 500, 600),
0407 'F': (30, 40, 50, 80, 120, 240, 600)}
0408
0409 if args.jet == 'hltAK4PFPuppiJets': JetType = "AK4 PF Puppi Jets"
0410 elif args.jet == 'hltAK4PFClusterJets': JetType = "AK4 PF Cluster Jets"
0411 elif args.jet == 'hltAK4PFJets': JetType = "AK4 PF Jets"
0412 elif args.jet == 'hltAK4PFCHSJets': JetType = "AK4 PF CHS Jets"
0413 else: JetType = args.jet
0414
0415 colors = hep.style.CMS['axes.prop_cycle'].by_key()['color']
0416 markers = ('o', 's', 'd')
0417 errorbar_kwargs = dict(capsize=3, elinewidth=0.8, capthick=2, linewidth=2, linestyle='')
0418
0419
0420
0421
0422
0423 Var1DList = {
0424 'HLT Jets' : 'JetPt',
0425 'HLT Corrected Jets' : 'CorrJetPt',
0426 'Gen-level Jets' : 'GenPt',
0427 'Photon Multiplicity' : 'photonMultiplicity',
0428 'Neutral Multiplicity' : 'neutralMultiplicity',
0429 'Charged Multiplicity' : 'chargedMultiplicity',
0430 'Neutral Hadron Multiplicity' : 'neutralHadronMultiplicity',
0431 'Charged Hadron Multiplicity' : 'chargedHadronMultiplicity',
0432 }
0433
0434 for Label, Var in Var1DList.items():
0435 plotter = Plotter(args.sample_label)
0436 root_hist = CheckRootFile(f"{dqm_dir}/{Var}", rebin=None)
0437 nbins, bin_edges, bin_centers, bin_widths = define_bins(root_hist)
0438 values, errors = histo_values_errors(root_hist)
0439
0440 plt.errorbar(bin_centers, values, xerr=None, yerr=errors,
0441 fmt='s', color='black', label=Label, **errorbar_kwargs)
0442 plt.step(bin_edges[:-1], values, where="post", color='black')
0443 plotter.ax.text(0.03, 0.97, f"{JetType}", transform=plotter.ax.transAxes, fontsize=fontsize,
0444 verticalalignment='top', horizontalalignment='left')
0445
0446 if 'Multiplicity' in Var:
0447 plotter.limits(y=(5E-1, 2*max(values)), logY=True)
0448 plotter.labels(x=Label, y='# Jets')
0449 else:
0450 plotter.limits(y=(0, 1.2*max(values)), logY=False)
0451 plotter.labels(x=HLabels.pt_label('pt'), y='# Jets', legend_title='')
0452
0453 plotter.save( os.path.join(args.odir, Var) )
0454
0455
0456
0457
0458
0459 if not os.path.exists(histo2d_dir):
0460 os.makedirs(histo2d_dir)
0461
0462 Var2DList = ('h2d_PtRecoOverGen_nCost_B', 'h2d_PtRecoOverGen_nCost_E', 'h2d_PtRecoOverGen_nCost_F',
0463 'h2d_PtRecoOverGen_chHad_B', 'h2d_PtRecoOverGen_chHad_E', 'h2d_PtRecoOverGen_chHad_F',
0464 'h2d_PtRecoOverGen_neHad_B', 'h2d_PtRecoOverGen_neHad_E', 'h2d_PtRecoOverGen_neHad_F',
0465 'h2d_PtRecoOverGen_chEm_B', 'h2d_PtRecoOverGen_chEm_E', 'h2d_PtRecoOverGen_chEm_F',
0466 'h2d_PtRecoOverGen_neEm_B', 'h2d_PtRecoOverGen_neEm_E', 'h2d_PtRecoOverGen_neEm_F',
0467 )
0468
0469 for Var2D in Var2DList:
0470 plotter = Plotter(args.sample_label, fontsize=15)
0471 root_hist = CheckRootFile(f"{dqm_dir}/{Var2D}")
0472
0473 nbins_x, nbins_y, x_edges, y_edges = define_bins_2D(root_hist)
0474 values = histo_values_2D(root_hist)
0475
0476 ylabel = root_hist.GetYaxis().GetTitle().replace('#', '\\')
0477
0478
0479
0480 pcm = plotter.ax.pcolormesh(x_edges, y_edges, np.where(values==0, np.nan, values),
0481 cmap='viridis', shading='auto')
0482
0483 for lab2d in ('nCost', 'chHad', 'neHad', 'chEm', 'neEm'):
0484 if lab2d in Var2D:
0485 xlabel = HLabels.fraction_label(lab2d)
0486
0487 plotter.labels(x=xlabel, y=f"${ylabel}$")
0488 plotter.fig.colorbar(pcm, ax=plotter.ax, label='# Jets')
0489
0490 plotter.save( os.path.join(histo2d_dir, Var2D) )
0491
0492
0493 nsigmas = (1, 2)
0494 ideal_fraction = {1: 0.3173, 2: 0.0455, 3: 0.0027}
0495 pt_bins = np.array((20, 30, 40, 60, 90, 150, 250, 400, 650, 1000))
0496 Var2DList = ('h2d_PtRecoOverGen_GenPt', )
0497
0498 for Var2D in Var2DList:
0499 plotter = Plotter(args.sample_label, fontsize=15)
0500 root_hist = CheckRootFile(f"{dqm_dir}/{Var2D}")
0501
0502 tail_fracs, tail_low_errors, tail_high_errors = ({ns:[] for ns in nsigmas} for _ in range(3))
0503 for ipt, (low,high) in enumerate(zip(pt_bins[:-1],pt_bins[1:])):
0504 if root_hist.GetXaxis().FindBin(low) <= 0 and root_hist.GetXaxis().FindBin(high) >= root_hist.GetNbinsX():
0505 mess = f"Low bin: {low} | Low value: {root_hist.GetXaxis().FindBin(low)} | High bin: {high} | High value: {root_hist.GetXaxis().FindBin(high)}"
0506 raise RuntimeError(mess)
0507
0508 hproj = root_hist.ProjectionY(root_hist.GetName() + "_proj" + str(ipt),
0509 root_hist.GetXaxis().FindBin(low), root_hist.GetXaxis().FindBin(high), "e")
0510
0511 integr_start, integr_stop = 1, hproj.GetNbinsX()+1
0512 integr = hproj.Integral(integr_start, integr_stop)
0513 gausTF1 = findBestGaussianCoreFit(hproj)
0514
0515
0516 plotter_single = Plotter(args.sample_label, fontsize=15)
0517 nbins, bin_edges, bin_centers, bin_widths = define_bins(hproj)
0518 values, errors = histo_values_errors(hproj)
0519 pt_label = str(low) + " < $p_T$ < " + str(high) + " GeV"
0520 plotter_single.ax.errorbar(bin_centers, values, xerr=None, yerr=errors/2, fmt='o', markersize=3,
0521 color='black', label=pt_label, **errorbar_kwargs)
0522
0523 for ins, ns in enumerate(nsigmas):
0524 xfunc = np.linspace(gausTF1.GetParameter(1) - ns*abs(gausTF1.GetParameter(2)), gausTF1.GetParameter(1) + ns*abs(gausTF1.GetParameter(2)))
0525 yfunc = np.array([gausTF1.Eval(xi) for xi in xfunc])
0526 plotter_single.ax.plot(xfunc, yfunc, 'o', color=colors[ins], linewidth=2, linestyle='--', label=f'{ns}$\sigma$ gaussian coverage')
0527
0528 tail_low = hproj.Integral(integr_start, hproj.FindBin(gausTF1.GetParameter(1) - ns*abs(gausTF1.GetParameter(2))))
0529 tail_high = hproj.Integral(hproj.FindBin(gausTF1.GetParameter(1) + ns*abs(gausTF1.GetParameter(2))), integr_stop)
0530 tail_fracs[ns].append((tail_low + tail_high) / integr)
0531 tail_frac_errlo, tail_frac_errhi = tails_errors(integr, tail_low + tail_high)
0532 tail_low_errors[ns].append(tail_frac_errlo)
0533 tail_high_errors[ns].append(tail_frac_errhi)
0534 xlabel = hproj.GetXaxis().GetTitle().replace('#', '\\')
0535 plotter_single.labels(x=f"{xlabel}", y=f"# Jets", legend_title='')
0536 plotter_single.save( os.path.join(args.odir, Var2D + '_tails_fit' + str(ipt)) )
0537
0538 pt_centers = pt_bins[:-1] + (pt_bins[1:] - pt_bins[:-1])/2
0539 for ins, ns in enumerate(nsigmas):
0540 plotter.ax.errorbar(pt_centers, tail_fracs[ns], xerr=None, yerr=[tail_low_errors[ns],tail_low_errors[ns]],
0541 fmt='o', markersize=3, color=colors[ins], **errorbar_kwargs)
0542 plotter.ax.stairs(tail_fracs[ns], pt_bins, color=colors[ins], linewidth=2)
0543 plotter.ax.axhline(y=ideal_fraction[ns], color=colors[ins], linestyle='--', label=f'{ns}$\sigma$ coverage')
0544 plotter.ax.set_xscale('log')
0545 xlabel = root_hist.GetXaxis().GetTitle().replace('#', '\\')
0546 plotter.labels(x=f"{xlabel}", y="Resolution Tail Fraction", legend_title='Ideal Gaussian')
0547 plotter.save( os.path.join(args.odir, Var2D + '_tails') )
0548
0549
0550 Var2DList = (
0551 'h2d_chEm_pt_B', 'h2d_chEm_pt_E', 'h2d_chEm_pt_F',
0552 'h2d_neEm_pt_B', 'h2d_neEm_pt_E', 'h2d_neEm_pt_F',
0553 'h2d_chHad_pt_B', 'h2d_chHad_pt_E', 'h2d_chHad_pt_F',
0554 'h2d_neHad_pt_B', 'h2d_neHad_pt_E', 'h2d_neHad_pt_F',
0555 'h2d_chMult_pt_B', 'h2d_chMult_pt_E', 'h2d_chMult_pt_F',
0556 'h2d_neMult_pt_B', 'h2d_neMult_pt_E', 'h2d_neMult_pt_F',
0557 'h2d_chHadMult_pt_B', 'h2d_chHadMult_pt_E', 'h2d_chHadMult_pt_F',
0558 'h2d_neHadMult_pt_B', 'h2d_neHadMult_pt_E', 'h2d_neHadMult_pt_F',
0559 'h2d_phoMult_pt_B', 'h2d_phoMult_pt_E', 'h2d_phoMult_pt_F',
0560 )
0561 pt_bins = (20, 60, 100, 1000)
0562
0563 for Var2D in Var2DList:
0564 root_hist = CheckRootFile(f"{dqm_dir}/{Var2D}")
0565
0566 plotter = Plotter(args.sample_label, fontsize=15)
0567
0568 nbins_x, nbins_y, x_edges, y_edges = define_bins_2D(root_hist)
0569 values = histo_values_2D(root_hist)
0570
0571 xlabel = root_hist.GetXaxis().GetTitle().replace('#', '\\')
0572
0573
0574
0575 pcm = plotter.ax.pcolormesh(x_edges, y_edges, np.where(values==0, np.nan, values),
0576 cmap='viridis', shading='auto')
0577
0578 for lab2d in ('nCost', 'chHad', 'neHad', 'chEm', 'neEm', 'chMult', 'neMult'):
0579 if lab2d in Var2D:
0580 if 'Mult' in Var2D:
0581 ylabel = HLabels.multiplicity_label(lab2d)
0582 else:
0583 ylabel = HLabels.fraction_label(lab2d)
0584 break
0585
0586 plotter.labels(x=f"${xlabel}$", y=ylabel)
0587 plotter.fig.colorbar(pcm, ax=plotter.ax, label='# Jets')
0588 plotter.save(os.path.join(histo2d_dir, Var2D))
0589
0590 if 'Mult' not in Var2D:
0591 continue
0592
0593 for ptbin in pt_bins:
0594 notfound = True
0595 for ibin in range(root_hist.GetNbinsX()+1):
0596 if root_hist.GetXaxis().GetBinLowEdge(ibin+1) == ptbin:
0597 notfound = False
0598 break
0599 if notfound:
0600 raise RuntimeError(f"The specified pT bin '{ptbin}' could not be matched to histogram {root_hist.GetName()}.")
0601
0602 hproj = {}
0603 values_all, errors_all, labels_all = [], [], []
0604 nybins = root_hist.GetNbinsY()
0605 bin_edges = np.array(root_hist.GetXaxis().GetXbins())
0606 for ipt, (low,high) in enumerate(zip(pt_bins[:-1],pt_bins[1:])):
0607 plotter_single = Plotter(args.sample_label, fontsize=15)
0608 hname = root_hist.GetName() + '_proj' + str(ipt)
0609 htitle = root_hist.GetTitle() + ' Proj PtBin' + str(ipt)
0610 if bin_edges.size:
0611 hproj[ipt] = ROOT.TH1F(hname, htitle, nybins, bin_edges)
0612 else:
0613 hproj[ipt] = ROOT.TH1F(hname, htitle, nybins,
0614 root_hist.GetXaxis().GetBinLowEdge(1), root_hist.GetXaxis().GetBinLowEdge(nybins+1))
0615
0616 for ibin in range(root_hist.GetNbinsX()+1):
0617 xlow = root_hist.GetXaxis().GetBinLowEdge(ibin+1)
0618 xhigh = root_hist.GetXaxis().GetBinLowEdge(ibin+2)
0619 if low <= xlow and high >= xhigh:
0620 for jbin in range(root_hist.GetNbinsX()+1):
0621
0622 val_curr = hproj[ipt].GetBinContent(jbin+1)
0623 val_next = root_hist.GetBinContent(ibin+1, jbin+1)
0624 error_curr = val_curr*hproj[ipt].GetBinError(jbin+1)
0625 error_next = val_next*root_hist.GetBinError(ibin+1, jbin+1)
0626 error = 0. if val_curr+val_next==0 else (error_curr + error_next) / (val_curr + val_next)
0627 hproj[ipt].SetBinError(jbin+1, error)
0628 hproj[ipt].SetBinContent(jbin+1, val_curr + val_next)
0629
0630 nbins, bin_edges, bin_centers, bin_widths = define_bins(hproj[ipt])
0631 values, errors = histo_values_errors(hproj[ipt])
0632 pt_label = str(low) + " < $p_T$ < " + str(high) + " GeV"
0633 values_all.append(values)
0634 errors_all.append(errors)
0635 labels_all.append(pt_label)
0636 plotter_single.ax.errorbar(bin_centers, values, xerr=None, yerr=errors/2, fmt='o', markersize=3,
0637 color='black', label=pt_label, **errorbar_kwargs)
0638 plotter_single.ax.stairs(values, bin_edges, color='black', linewidth=2)
0639 plotter_single.ax.set_yscale('log')
0640 plotter_single.ax.text(0.97, 0.92, f"{EtaInfo.label(Var2D[-1])}", transform=plotter_single.ax.transAxes, fontsize=fontsize,
0641 verticalalignment='top', horizontalalignment='right')
0642 plotter_single.labels(x=ylabel, y="# Jets", legend_title='')
0643 plotter_single.save(os.path.join(args.odir, Var2D + '_PtBin' + str(ipt)))
0644
0645 plotter_all = Plotter(args.sample_label, fontsize=15)
0646 for idx, (vals, errs, lab) in enumerate(zip(values_all, errors_all, labels_all)):
0647 values = np.zeros_like(vals) if sum(vals)==0 else vals / sum(vals)
0648 errors = np.zeros_like(vals) if sum(vals)==0 else errs / sum(vals)
0649 plotter_all.ax.errorbar(bin_centers, values, xerr=None, yerr=errors/2, fmt='o', markersize=3,
0650 color=colors[idx], label=lab, **errorbar_kwargs)
0651 plotter_all.ax.stairs(values, bin_edges, color=colors[idx], linewidth=2)
0652 plotter_all.ax.set_yscale('log')
0653 plotter_all.ax.text(0.97, 0.79, f"{EtaInfo.label(Var2D[-1])}", transform=plotter_all.ax.transAxes, fontsize=fontsize,
0654 verticalalignment='top', horizontalalignment='right')
0655 plotter_all.labels(x=ylabel, y="Counts [a.u.]", legend_title='Jet $p_T$ range')
0656 plotter_all.save(os.path.join(args.odir, Var2D + '_PtBinAll'))
0657
0658
0659
0660
0661
0662 GroupedVarList = dotdict({
0663 'JetPt_EtaRegions': dotdict(
0664 histos={'Barrel': 'JetPt_B', 'Endcap': 'JetPt_E', 'Forward': 'JetPt_F'},
0665 xlabel=HLabels.pt_label('pt'),
0666 ),
0667 'GenPt_EtaRegions': dotdict(
0668 histos={'Barrel': 'GenPt_B', 'Endcap': 'GenPt_E', 'Forward': 'GenPt_F'},
0669 xlabel=HLabels.pt_label('gen'),
0670 ),
0671 'JetTypes_Pt': dotdict(
0672 histos={'Gen-level Jets': 'GenPt', 'HLT Jets': 'JetPt', 'HLT Corrected Jets': 'CorrJetPt'},
0673 xlabel=HLabels.pt_label('pt'),
0674 ),
0675 'JetChargedHadFrac': dotdict(
0676 histos={'HLT jets': 'chargedHadronEnergyFraction',
0677 'HLT jets matched': 'MatchedJetchHad'},
0678 xlabel=HLabels.fraction_label('chHad'),
0679 ),
0680 'JetNeutralHadFrac': dotdict(
0681 histos={'HLT jets': 'neutralHadronEnergyFraction',
0682 'HLT jets matched': 'MatchedJetneHad'},
0683 xlabel=HLabels.fraction_label('neHad')
0684 ),
0685 'JetChargedEmFrac': dotdict(
0686 histos={'HLT jets': 'chargedEmEnergyFraction',
0687 'HLT jets matched': 'MatchedJetchEm'},
0688 xlabel=HLabels.fraction_label('chEm')
0689 ),
0690 'JetNeutralEmFrac': dotdict(
0691 histos={'HLT jets': 'neutralEmEnergyFraction',
0692 'HLT jets matched': 'MatchedJetneEm',},
0693 xlabel=HLabels.fraction_label('neEm')
0694 ),
0695 'JetnConst': dotdict(
0696 histos={'HLT jets': 'JetConstituents',
0697 'HLT jets matched': 'MatchedJetnCost'},
0698 xlabel=HLabels.fraction_label('nCost')
0699 ),
0700 'photonMultiplicity': dotdict(
0701 histos={'Barrel': 'photonMultiplicity_B',
0702 'Endcap': 'photonMultiplicity_E',
0703 'Forward': 'photonMultiplicity_F'},
0704 xlabel="Photon Multiplicity"
0705 ),
0706 'neutralMultiplicity': dotdict(
0707 histos={'Barrel': 'neutralMultiplicity_B',
0708 'Endcap': 'neutralMultiplicity_E',
0709 'Forward': 'neutralMultiplicity_F'},
0710 xlabel="Neutral Multiplicity"
0711 ),
0712 'chargedMultiplicity': dotdict(
0713 histos={'Barrel': 'chargedMultiplicity_B',
0714 'Endcap': 'chargedMultiplicity_E',
0715 'Forward': 'chargedMultiplicity_F'},
0716 xlabel="Charged Multiplicity"
0717 ),
0718 'chargedHadronMultiplicity': dotdict(
0719 histos={'Barrel': 'chargedHadronMultiplicity_B',
0720 'Endcap': 'chargedHadronMultiplicity_E',
0721 'Forward': 'chargedHadronMultiplicity_F'},
0722 xlabel="Charged Hadron Multiplicity"
0723 ),
0724 'neutralHadronMultiplicity': dotdict(
0725 histos={'Barrel': 'neutralHadronMultiplicity_B',
0726 'Endcap': 'neutralHadronMultiplicity_E',
0727 'Forward': 'neutralHadronMultiplicity_F'},
0728 xlabel="Neutral Hadron Multiplicity"
0729 ),
0730 })
0731
0732 for GroupedVar in GroupedVarList:
0733 plotter = Plotter(args.sample_label)
0734
0735 for i_var, (Label, Var) in enumerate(GroupedVarList[GroupedVar]['histos'].items()):
0736 root_hist = CheckRootFile(f"{dqm_dir}/{Var}", rebin=None)
0737 nbins, bin_edges, bin_centers, bin_widths = define_bins(root_hist)
0738 values, errors = histo_values_errors(root_hist)
0739
0740 plt.errorbar(bin_centers, values, xerr=None, yerr=errors,
0741 label=Label, color=colors[i_var], fmt=markers[i_var], **errorbar_kwargs)
0742 plt.step(bin_edges[:-1], values, where="post", color=colors[i_var], linewidth=2)
0743
0744 plotter.labels(x=GroupedVarList[GroupedVar].xlabel, y='# Jets' if 'Multiplicity' in GroupedVar else '# Jets', legend_title='')
0745 plotter.save( os.path.join(args.odir, GroupedVar) )
0746
0747
0748
0749
0750
0751 for i_res, ResType in enumerate(('PtRecoOverGen', 'PtCorrOverGen', 'PtCorrOverReco')):
0752
0753 for EtaRegion in EtaInfo.regions():
0754
0755 MergedPtBins = {
0756 '20 < $p_T$ < 40 GeV' : [f'h_{ResType}_{EtaRegion}_Pt20_30', f'h_{ResType}_{EtaRegion}_Pt30_40'],
0757 '40 < $p_T$ < 100 GeV' : [f'h_{ResType}_{EtaRegion}_Pt40_100'],
0758 '100 < $p_T$ < 300 GeV' : [f'h_{ResType}_{EtaRegion}_Pt100_200', f'h_{ResType}_{EtaRegion}_Pt200_300'],
0759 '300 < $p_T$ < 600 GeV' : [f'h_{ResType}_{EtaRegion}_Pt300_600',],
0760 '600 < $p_T$ < 5000 GeV': [f'h_{ResType}_{EtaRegion}_Pt600_2000', f'h_{ResType}_{EtaRegion}_Pt2000_5000']
0761 }
0762
0763 v_stacked_histo = []
0764 v_labels = []
0765 for i, (pt_label, hist_names) in enumerate(MergedPtBins.items()):
0766 stacked_hist = None
0767
0768 for j, hist_name in enumerate(hist_names):
0769 rebin = 2 if EtaRegion != 'F' else 4
0770 root_hist = CheckRootFile(f"{dqm_dir}/{hist_name}", rebin=rebin)
0771
0772 if stacked_hist is None:
0773 stacked_hist = root_hist.Clone()
0774 else:
0775 stacked_hist.Add(root_hist)
0776
0777 if stacked_hist is None: continue
0778 v_stacked_histo.append(stacked_hist)
0779 v_labels.append(pt_label)
0780
0781 if len(v_stacked_histo) == 0: continue
0782
0783 plotter = Plotter(args.sample_label)
0784 for i, (stacked_histo, pt_label) in enumerate(zip(v_stacked_histo, v_labels)):
0785 if stacked_histo.Integral() == 0:
0786 print(f"WARNING: Skipping empty histogram for {pt_label}")
0787 continue
0788 if stacked_histo.Integral() < 2:
0789 print(f"WARNING: Skipping histogram with low stat {pt_label}")
0790 continue
0791
0792 stacked_histo.Scale(1.0 / stacked_histo.Integral())
0793 nbins, bin_edges, bin_centers, bin_widths = define_bins(stacked_histo)
0794 values, errors = histo_values_errors(stacked_histo)
0795
0796 plotter.ax.hist(bin_edges[:-1], bins=bin_edges, weights=values, histtype='stepfilled', color=colors[i], alpha=0.1)
0797 plotter.ax.hist(bin_edges[:-1], bins=bin_edges, weights=values, histtype='step', color=colors[i], linewidth=1.5)
0798 plotter.ax.errorbar(bin_centers, values, xerr=None, yerr=errors, fmt='o', markersize=3,
0799 color=colors[i], label=pt_label, **errorbar_kwargs)
0800
0801 plotter.ax.text(0.03, 0.97, f"{JetType}\n{EtaInfo.label(EtaRegion)}", transform=plotter.ax.transAxes, fontsize=fontsize,
0802 verticalalignment='top', horizontalalignment='left')
0803
0804 plotter.labels(x="${}$".format(v_stacked_histo[0].GetXaxis().GetTitle()),
0805 y="Counts [a.u.]",
0806 legend_title=r"Jet $p_T$ range")
0807
0808 plotter.save( os.path.join(args.odir, f"Response_{ResType}_{EtaRegion}") )
0809
0810
0811
0812
0813
0814 ResolOptions = {
0815 'Scale' : (0.0, 2.7, '#7a21dd'),
0816 'Resolution' : (0.0, 0.8, '#5790fc')
0817 }
0818
0819 for key in ResolOptions.keys():
0820 for resol_type in HLabels.resol_types():
0821 myResolLabel = HLabels(resol_type)
0822 for myXvar in myResolLabel.xvars:
0823
0824 plotter = Plotter(args.sample_label)
0825 tprofile_mean = CheckRootFile( myResolLabel.mhisto(myXvar, dqm_dir), rebin=None)
0826 tprofile_sigma = CheckRootFile( myResolLabel.shisto(myXvar, dqm_dir), rebin=None)
0827 nbins, bin_edges, bin_centers, bin_widths = define_bins(tprofile_mean)
0828 means, mean_errors = histo_values_errors(tprofile_mean)
0829 sigmas, sigma_errors = histo_values_errors(tprofile_sigma)
0830 if key == 'Scale':
0831 y, y_errors = means, mean_errors
0832 ylabel = myResolLabel.ytitle(resol=False)
0833 else:
0834 y = [s / m if m != 0 else np.nan for s, m in zip(sigmas, means)]
0835 y_errors = [np.sqrt((ds / m)**2 + ((s / m) * (dm / m))**2) if m != 0 else np.nan
0836 for s, ds, m, dm in zip(sigmas, sigma_errors, means, mean_errors)]
0837 ylabel = myResolLabel.ytitle(resol=True)
0838
0839 plotter.ax.errorbar(bin_centers, y, xerr=None, yerr=y_errors,
0840 fmt='o', color=ResolOptions[key][2], label=f'{key} {resol_type}', **errorbar_kwargs)
0841
0842 if 'Pt' not in myXvar:
0843 xlabel = fr'$\{myXvar.lower()}$'
0844 else:
0845 xlabel = HLabels.pt_label('gen') if 'Gen' in resol_type else HLabels.pt_label('reco')
0846
0847 plotter.labels(x=xlabel, y=ylabel, legend_title='')
0848 plotter.limits(y=(ResolOptions[key][0], ResolOptions[key][1]))
0849 Text = f"{JetType}\n{EtaInfo.label(myXvar[-1])}" if any(x in myXvar for x in ('_B', '_E', '_F')) else f"{JetType}"
0850 plotter.ax.text(0.03, 0.97, Text, transform=plotter.ax.transAxes,
0851 fontsize=fontsize, verticalalignment='top', horizontalalignment='left')
0852
0853 if key == 'Scale':
0854 plotter.ax.axhline(1.0, color='gray', linestyle='--', linewidth=2)
0855
0856 plotter.save( os.path.join(args.odir, f'{key}_{resol_type}_{myXvar}') )
0857
0858 plotter = Plotter(args.sample_label, fontsize=fontsize)
0859 for etareg in EtaInfo.regions():
0860 tprofile_mean = CheckRootFile( myResolLabel.mhisto(f"Pt_{etareg}", dqm_dir), rebin=tprofile_rebinning[etareg] )
0861 tprofile_sigma = CheckRootFile( myResolLabel.shisto(f"Pt_{etareg}", dqm_dir), rebin=tprofile_rebinning[etareg] )
0862 nbins, bin_edges, bin_centers, bin_widths = define_bins(tprofile_mean)
0863 means, mean_errors = histo_values_errors(tprofile_mean)
0864 sigmas, sigma_errors = histo_values_errors(tprofile_sigma)
0865 if key == 'Scale':
0866 y, y_errors = means, mean_errors
0867 ylabel = myResolLabel.ytitle(resol=False)
0868 else:
0869 y = [s / m if m != 0 else np.nan for s, m in zip(sigmas, means)]
0870 y_errors = [np.sqrt((ds / m)**2 + ((s / m) * (dm / m))**2) if m != 0 else np.nan
0871 for s, ds, m, dm in zip(sigmas, sigma_errors, means, mean_errors)]
0872 ylabel = myResolLabel.ytitle(resol=True)
0873
0874 plt.errorbar(bin_centers, y, xerr=None, yerr=y_errors,
0875 fmt=EtaInfo.marker(etareg), color=EtaInfo.color(etareg), label=EtaInfo.label(etareg),
0876 elinewidth=0.8, linewidth=2)
0877 plt.stairs(y, bin_edges, color=EtaInfo.color(etareg))
0878
0879 xlabel = HLabels.pt_label('gen') if 'Gen' in resol_type else HLabels.pt_label('reco')
0880 plotter.labels(x=xlabel, y=ylabel, legend_title='')
0881 plotter.limits(y=(ResolOptions[key][0], ResolOptions[key][1]))
0882 plotter.ax.text(0.03, 0.97, f"{JetType}", transform=plotter.ax.transAxes,
0883 fontsize=fontsize, verticalalignment='top', horizontalalignment='left')
0884
0885 if key == 'Scale':
0886 plotter.ax.axhline(1.0, color='gray', linestyle='--', linewidth=2)
0887
0888 plotter.save( os.path.join(args.odir, f'{key}_{resol_type}_Pt_EtaBins') )
0889
0890
0891
0892
0893
0894 for key in ResolOptions.keys():
0895 plotter = Plotter(args.sample_label, fontsize=fontsize)
0896 for i_res, resol_type in enumerate(['RecoOverGen', 'CorrOverGen']):
0897 myResolLabel = HLabels(resol_type)
0898 for etareg in EtaInfo.regions():
0899 tprofile_mean = CheckRootFile( myResolLabel.mhisto(f"Pt_{etareg}", dqm_dir), rebin=tprofile_rebinning[etareg] )
0900 tprofile_sigma = CheckRootFile( myResolLabel.shisto(f"Pt_{etareg}", dqm_dir), rebin=tprofile_rebinning[etareg] )
0901 nbins, bin_edges, bin_centers, bin_widths = define_bins(tprofile_mean)
0902 means, mean_errors = histo_values_errors(tprofile_mean)
0903 sigmas, sigma_errors = histo_values_errors(tprofile_sigma)
0904 if key == 'Scale':
0905 y, y_errors = means, mean_errors
0906 else:
0907 y = [s / m if m != 0 else np.nan for s, m in zip(sigmas, means)]
0908 y_errors = [np.sqrt((ds / m)**2 + ((s / m) * (dm / m))**2) if m != 0 else np.nan
0909 for s, ds, m, dm in zip(sigmas, sigma_errors, means, mean_errors)]
0910
0911
0912 mfc = 'white' if i_res == 1 else EtaInfo.color(etareg)
0913 eb = plotter.ax.errorbar(bin_centers, y, xerr=0.5 * bin_widths, yerr=y_errors, mfc=mfc,
0914 fmt=EtaInfo.marker(etareg), color=EtaInfo.color(etareg), label=EtaInfo.label(etareg),
0915 elinewidth=0.8, linewidth=2)
0916 eb[-1][0].set_linestyle('-' if i_res==0 else '--')
0917
0918 if key == 'Scale':
0919 ylabel = HLabels.pt_label('pt/gen', average=True)
0920 else:
0921 ylabel = HLabels.resol_label('pt/gen')
0922 plotter.labels(x=HLabels.pt_label('gen'), y=ylabel)
0923 plotter.limits(y=(ResolOptions[key][0], ResolOptions[key][1]))
0924 plotter.ax.text(0.03, 0.97, f"{JetType}", transform=plotter.ax.transAxes,
0925 fontsize=fontsize, verticalalignment='top', horizontalalignment='left')
0926
0927 from matplotlib.lines import Line2D
0928
0929 eta_legend_elements = [
0930 Line2D([0], [0], color=EtaInfo.color(x), marker=EtaInfo.marker(x), linestyle='-', label=EtaInfo.label(x))
0931 for x in EtaInfo.regions()
0932 ]
0933
0934
0935 res_legend_elements = [
0936 Line2D([0], [0], color='grey', marker='o', mfc='grey', linestyle='-', label='Reco'),
0937 Line2D([0], [0], color='grey', marker='o', mfc='white', linestyle='--', label='Corrected')
0938 ]
0939
0940 legend_eta = plotter.ax.legend(handles=eta_legend_elements, loc='upper right', fontsize=fontsize)
0941 legend_res = plotter.ax.legend(handles=res_legend_elements, loc='upper right', fontsize=fontsize, bbox_to_anchor=(0.73, 0.99))
0942 plotter.ax.add_artist(legend_eta)
0943 plotter.save( os.path.join(args.odir, f'Pt{key}_CorrVsReco') )
0944
0945
0946
0947
0948
0949 eff_color = '#bd1f01'
0950 for eff_type in HLabels.eff_types():
0951 myEffLabel = HLabels(eff_type)
0952 common_kwargs = dict(linestyle='', color=eff_color, label=eff_type)
0953 if any(x in myEffLabel.ytitle() for x in ('Fake', 'Duplicate')):
0954 if eff_type == 'Fake Rate':
0955 axmin = 1E-2
0956 else:
0957 axmin = 1E-4
0958 axmax = 2.4
0959 else:
0960 axmin, axmax = 0, 1.25
0961
0962 for myXvar in myEffLabel.xvars:
0963 root_hist_num = CheckRootFile( myEffLabel.nhisto(myXvar, dqm_dir), rebin=2 )
0964 root_hist_den = CheckRootFile( myEffLabel.dhisto(myXvar, dqm_dir), rebin=2 )
0965 root_ratio = CheckRootFile( myEffLabel.rhisto(myXvar, dqm_dir), rebin=2 )
0966
0967 nbins, bin_edges, bin_centers, bin_widths = define_bins(root_ratio)
0968 numerator_vals, _ = histo_values_errors(root_hist_num)
0969 denominator_vals, _ = histo_values_errors(root_hist_den)
0970 eff_values, eff_errors = histo_values_errors(root_ratio)
0971 if eff_type == 'Fake Rate':
0972 eff_values = np.array([1-i if i != 0 else np.nan for i in eff_values])
0973
0974 plotter = Plotter(args.sample_label, grid_color=None, fontsize=fontsize)
0975 stepkwargs = dict(where="post", linewidth=2)
0976 plotter.ax.step(bin_edges[:-1], denominator_vals, label=myEffLabel.leglabel(isNum=False), color="black", **stepkwargs)
0977 plotter.ax.step(bin_edges[:-1], numerator_vals, label=myEffLabel.leglabel(isNum=True), color="#9c9ca1", linestyle='-.', **stepkwargs)
0978 plotter.ax.fill_between(bin_edges[:-1], numerator_vals, step="post", alpha=0.3, color="#9c9ca1")
0979
0980 label = root_hist_num.GetXaxis().GetTitle().replace('#', '\\')
0981 plotter.labels(x=f"${label}$", y="# Jets", legend_title='', legend_loc='upper left')
0982 plotter.limits(y=(0, 1.2*max(denominator_vals)))
0983
0984 Text = f"{JetType}\n{EtaInfo.label(myXvar[-1])}" if any(x in myXvar for x in ('_B', '_E', '_F')) else f"{JetType}"
0985 plotter.ax.text(0.97, 0.97, Text, transform=plotter.ax.transAxes,
0986 fontsize=fontsize, verticalalignment='top', horizontalalignment='right')
0987
0988 ax2 = plotter.ax.twinx()
0989 ax2.set_ylabel(myEffLabel.ytitle(), color=eff_color)
0990
0991 if any(x in myEffLabel.ytitle() for x in ('Fake', 'Duplicate')):
0992 ax2.set_yscale('log')
0993 ax2.set_ylim(axmin, axmax)
0994
0995 eff_filt, (err_filt_lo, err_filt_hi), up_error, transform = rate_errorbar_declutter(eff_values, eff_errors, axmin)
0996 ax2.errorbar(bin_centers, eff_filt, xerr=0.5 * bin_widths, yerr=[err_filt_lo,err_filt_hi],
0997 fmt='o', capthick=2, linewidth=1, capsize=2, **common_kwargs)
0998 ax2.plot(bin_centers, up_error, 'v', transform=transform, **common_kwargs)
0999
1000 ax2.grid(color=eff_color, axis='y')
1001 ax2.tick_params(axis='y', labelcolor=eff_color)
1002 plotter.ax.grid(color=eff_color, axis='x')
1003
1004 plotter.save( os.path.join(args.odir, myEffLabel.savename + '_' + myXvar) )
1005
1006 plotter = Plotter(args.sample_label, fontsize=fontsize)
1007 for etareg in EtaInfo.regions():
1008 root_ratio = CheckRootFile( myEffLabel.rhisto(f"Pt_{etareg}", dqm_dir), rebin=tprofile_rebinning[etareg] )
1009 nbins, bin_edges, bin_centers, bin_widths = define_bins(root_ratio)
1010 eff_values, eff_errors = histo_values_errors(root_ratio)
1011 if eff_type == 'Fake Rate':
1012 eff_values = np.array([1-i if i != 0 else np.nan for i in eff_values])
1013
1014 eff_filt, (err_filt_lo, err_filt_hi), up_error, transform = rate_errorbar_declutter(eff_values, eff_errors, axmin)
1015 plotter.ax.errorbar(bin_centers, eff_filt, xerr=0.5 * bin_widths, yerr=[err_filt_lo,err_filt_hi],
1016 fmt=EtaInfo.marker(etareg), color=EtaInfo.color(etareg), label=EtaInfo.label(etareg),
1017 **errorbar_kwargs)
1018 plotter.ax.plot(bin_centers, up_error, 'v', linestyle='', color=EtaInfo.color(etareg), transform=transform)
1019
1020 plotter.ax.text(0.03, 0.97, f"{JetType}", transform=plotter.ax.transAxes,
1021 fontsize=fontsize, verticalalignment='top', horizontalalignment='left')
1022 label = root_ratio.GetXaxis().GetTitle()
1023 plotter.labels(x=f"${label}$", y=myEffLabel.ytitle(), legend_title='')
1024 if any(x in myEffLabel.ytitle() for x in ('Fake', 'Duplicate')):
1025 plotter.limits(y=(axmin, axmax), logY=True)
1026 else:
1027 plotter.limits(y=(0, axmax))
1028
1029 plotter.save( os.path.join(args.odir, myEffLabel.ytitle().replace(' ', '_') + '_Pt_EtaBins') )
