diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/create_input.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/create_input.py index f1b161dacf85c9de9a8ee08ee542c60e02e86e3a..7114879c35211ff0efefaa1ca94da2af69283060 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/create_input.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/create_input.py @@ -2,6 +2,7 @@ # this file do not work out of the box +import shutil import ROOT from array import array import numpy as np @@ -27,7 +28,8 @@ assert old_ct_sys file_christophe = ROOT.TFile("~/ElectronEnergyScaleFactor.root") scales_christophe = file_christophe.Get("alpha68") ct_christophe = file_christophe.Get("sigma24") -scales_sys_christophe = file_christophe.Get("systAlpha") # sum of 8 / 13 TeV diff + 64 / 32 bins diff +# sum of 8 / 13 TeV diff + 64 / 32 bins diff +scales_sys_christophe = file_christophe.Get("systAlpha") ct_sys_christophe = file_christophe.Get("systSigma") # 8 / 13 TeV diff assert scales_christophe @@ -39,7 +41,7 @@ assert ct_sys_christophe def qsum_histograms(histo1, histo2): new_histo = histo1.Clone() new_histo.Reset() - for ibin in xrange(histo1.GetNbinsX() + 2): + for ibin in range(histo1.GetNbinsX() + 2): value1 = histo1.GetBinContent(ibin) central_value = histo1.GetBinCenter(ibin) ibin2 = histo2.FindBin(central_value) @@ -127,23 +129,23 @@ def merge_histograms(old, new, merge_error=True): UNDERFLOW = 0 OVERFLOW = new.GetNbinsX() + 1 - for iold in xrange(1, old.GetNbinsX()): - l = old.GetBinLowEdge(iold) - r = l + old.GetBinWidth(iold) + for iold in range(1, old.GetNbinsX()): + low = old.GetBinLowEdge(iold) + r = low + old.GetBinWidth(iold) - il_new = new.FindFixBin(l) + il_new = new.FindFixBin(low) ir_new = new.FindFixBin(r) remainer = None if il_new == UNDERFLOW and ir_new == UNDERFLOW: - new_binning.append((l, r)) + new_binning.append((low, r)) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) elif il_new == UNDERFLOW and ir_new > UNDERFLOW and ir_new < OVERFLOW: - if abs(new.GetBinLowEdge(1) - l) < 1E-100: + if abs(new.GetBinLowEdge(1) - low) < 1E-100: continue - new_binning.append((l, new.GetBinLowEdge(1))) + new_binning.append((low, new.GetBinLowEdge(1))) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) if ir_new == OVERFLOW: @@ -151,27 +153,29 @@ def merge_histograms(old, new, merge_error=True): break last_old = iold - for inew in xrange(1, new.GetNbinsX() + 1): - l = new.GetBinLowEdge(inew) - r = l + new.GetBinWidth(inew) - new_binning.append((l, r)) + for inew in range(1, new.GetNbinsX() + 1): + low = new.GetBinLowEdge(inew) + r = low + new.GetBinWidth(inew) + new_binning.append((low, r)) new_values.append(new.GetBinContent(inew)) new_errors.append(new.GetBinError(inew)) if remainer is not None: - new_binning.append((new.GetBinLowEdge(new.GetNbinsX()), old.GetBinLowEdge(remainer) + old.GetBinWidth(remainer))) + new_binning.append((new.GetBinLowEdge(new.GetNbinsX()), + old.GetBinLowEdge(remainer) + + old.GetBinWidth(remainer))) new_values.append(old.GetBinContent(remainer)) new_errors.append(old.GetBinError(remainer)) - for iold in xrange(last_old, old.GetNbinsX() + 1): - l = old.GetBinLowEdge(iold) - r = l + old.GetBinWidth(iold) + for iold in range(last_old, old.GetNbinsX() + 1): + low = old.GetBinLowEdge(iold) + r = low + old.GetBinWidth(iold) - il_new = new.FindFixBin(l) + il_new = new.FindFixBin(low) ir_new = new.FindFixBin(r) if il_new == OVERFLOW and ir_new == OVERFLOW: - new_binning.append((l, r)) + new_binning.append((low, r)) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) elif il_new < OVERFLOW and ir_new == OVERFLOW: @@ -183,7 +187,8 @@ def merge_histograms(old, new, merge_error=True): new_edges = array('f', [x[0] for x in new_binning] + [new_binning[-1][1]]) histo_type = type(new) - result = histo_type(new.GetName(), new.GetTitle(), len(new_edges) - 1, new_edges) + result = histo_type(new.GetName(), new.GetTitle(), + len(new_edges) - 1, new_edges) for i, (v, e) in enumerate(zip(new_values, new_errors), 1): result.SetBinContent(i, v) if merge_error: @@ -223,15 +228,16 @@ histo_scale.SetName("alphaZee_errStat") histo_ct.SetName("ctZee_errStat") # this created a file structured as the official one, with empty directories -#output_file = create_structured_file("calibration_constants_run2.root") -import shutil +# output_file = create_structured_file("calibration_constants_run2.root") print(old_filename) shutil.copy2(old_filename, "xxx.root") output_file = ROOT.TFile("xxx.root", "update") create_new_directories(output_file) -output_file.GetDirectory("Scales").GetDirectory("es2015PRE").WriteObject(histo_scale, "alphaZee_errStat") -output_file.GetDirectory("Resolution").GetDirectory("es2015PRE").WriteObject(histo_ct, "ctZee_errStat") +output_file.GetDirectory("Scales").GetDirectory( + "es2015PRE").WriteObject(histo_scale, "alphaZee_errStat") +output_file.GetDirectory("Resolution").GetDirectory( + "es2015PRE").WriteObject(histo_ct, "ctZee_errStat") # systematics @@ -239,10 +245,13 @@ new_scale_sys = qsum_histograms(scales_sys_christophe, old_scale_sys) new_ct_sys = qsum_histograms(ct_sys_christophe, old_ct_sys) new_scale_sys = merge_histograms(old_scale_sys, new_scale_sys, False) new_ct_sys = merge_histograms(old_ct_sys, new_ct_sys, False) -new_scale_sys.SetTitle("es2015PRE scales sys = es2012c sys + 7/13 TeV diff + 34/68 bin diff") +new_scale_sys.SetTitle( + "es2015PRE scales sys = es2012c sys + 7/13 TeV diff + 34/68 bin diff") new_ct_sys.SetTitle("es2015 ct sys = es2012c sys + 7/13 TeV diff") -output_file.GetDirectory("Scales").GetDirectory("es2015PRE").WriteObject(new_scale_sys, "alphaZee_errSyst") -output_file.GetDirectory("Resolution").GetDirectory("es2015PRE").WriteObject(new_ct_sys, "ctZee_errSyst") +output_file.GetDirectory("Scales").GetDirectory( + "es2015PRE").WriteObject(new_scale_sys, "alphaZee_errSyst") +output_file.GetDirectory("Resolution").GetDirectory( + "es2015PRE").WriteObject(new_ct_sys, "ctZee_errSyst") legend3 = ROOT.TLegend(0.6, 0.7, 0.9, 0.9) legend3.SetBorderSize(0) @@ -277,8 +286,10 @@ legend4.Draw() # stefano input for sensitivity stefano_file = ROOT.TFile("egammaEnergyCorrectionDataMC15.root") -copy_dir(stefano_file.GetDirectory("PSRecalibration"), output_file.GetDirectory("PSRecalibration").GetDirectory("es2015PRE")) -copy_dir(stefano_file.GetDirectory("S1Recalibration"), output_file.GetDirectory("S1Recalibration").GetDirectory("es2015PRE")) +copy_dir(stefano_file.GetDirectory("PSRecalibration"), + output_file.GetDirectory("PSRecalibration").GetDirectory("es2015PRE")) +copy_dir(stefano_file.GetDirectory("S1Recalibration"), + output_file.GetDirectory("S1Recalibration").GetDirectory("es2015PRE")) # correction for uA2MeV first week 2015 @@ -286,4 +297,4 @@ f_ua2mev = ROOT.TFile.Open("~/uA2MeV.root") histo_ua2mev = f_ua2mev.Get("histo_uA2MeV_week12") output_file.Get("Scales").Get("es2015PRE").WriteTObject(histo_ua2mev) -raw_input() +input("Press Key") diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/dump_layer.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/dump_layer.py index 42f5aa849b457a45faf34b41d7af5b542d8a581c..53d0d86e27080abcba14a29a6ce742c492aad32b 100644 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/dump_layer.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/dump_layer.py @@ -31,26 +31,34 @@ def main(path, particle): layers = "rawcl_Es0", "rawcl_Es1", "rawcl_Es2", "rawcl_Es3" layers = [prefix + "_" + v for v in layers] - title = ["ratio_Es0_true_E", "ratio_Es1_true_E", "ratio_Es2_true_E", "ratio_Es3_true_E"] + title = ["ratio_Es0_true_E", + "ratio_Es1_true_E", + "ratio_Es2_true_E", + "ratio_Es3_true_E"] aetaCalo = "abs(" + prefix + "_" + "cl_etaCalo" + ")" truth_E = prefix + "_truth_E" - ratio_layers = ["%s / %s" % (l, truth_E) for l in layers] + ratio_layers = ["%s / %s" % (layer, truth_E) for layer in layers] truth_eta = prefix + "_truth_eta" truth_pt = "(%s / cosh(%s))" % (truth_E, truth_eta) - matched = "abs(({0}_rawcl_Es0 + {0}_rawcl_Es1 + {0}_rawcl_Es2 + {0}_rawcl_Es3)/{0}_truth_E - 1)<1".format(prefix) + matched = ("abs(({0}_rawcl_Es0 + " + "{0}_rawcl_Es1 + {0}_rawcl_Es2 + " + "{0}_rawcl_Es3)/{0}_truth_E - 1)<1").format(prefix) vars_to_plot = ratio_layers selection += " && " + "(" + matched + ")" pt_binning = [0, 5E3, 10E3, 15E3, 20E3, - 30E3, 35E3, 40E3, 45E3, 50E3, 55E3, 60E3, 65E3, 70E3, 75E3, 80E3, 85E3, 90E3, + 30E3, 35E3, 40E3, 45E3, 50E3, 55E3, 60E3, 65E3, + 70E3, 75E3, 80E3, 85E3, 90E3, 100E3, 110E3, 120E3, 130E3, 140E3, 160E3, 180E3, 200E3, 220E3, 250E3, 300E3, 350E3, 400E3, 450E3, 500E3, 550E3, 600E3, 700E3, 800E3, 900E3, 1000E3, 1200E3, 1400E3, 1600E3, 1800E3, 2000E3, 2500E3] - aeta_binning = [0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 1.0, - 1.1, 1.2, 1.3, 1.425, 1.550, 1.6, 1.65, 1.7, 1.75, 1.8, 1.9, 2.0, 2.1, + aeta_binning = [0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, + 0.7, 0.75, 0.8, 0.85, 0.9, 1.0, + 1.1, 1.2, 1.3, 1.425, 1.550, 1.6, + 1.65, 1.7, 1.75, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.35, 2.4, 2.5] pt_binning = array('d', pt_binning) aeta_binning = array('d', aeta_binning) @@ -62,7 +70,8 @@ def main(path, particle): print("plotting counting") histo_name = "histo_count_%s" % particle - histo_count = ROOT.TH2F(histo_name, "count %s" % particle, len(pt_binning) - 1, pt_binning, len(aeta_binning) - 1, aeta_binning) + histo_count = ROOT.TH2F(histo_name, "count %s" % particle, len( + pt_binning) - 1, pt_binning, len(aeta_binning) - 1, aeta_binning) chain.Draw(aetaCalo + ":" + truth_pt + ">>" + histo_name, selection) result = [] @@ -73,7 +82,8 @@ def main(path, particle): sel = selection + " && abs(ph_zconv) < 5000" else: sel = selection - histo = ROOT.TProfile2D(histo_name, t, len(pt_binning) - 1, pt_binning, len(aeta_binning) - 1, aeta_binning) + histo = ROOT.TProfile2D(histo_name, t, len( + pt_binning) - 1, pt_binning, len(aeta_binning) - 1, aeta_binning) arg = v + ":" + aetaCalo + ":" + truth_pt + ">>" + histo_name chain.Draw(arg, sel, "profcolz") histo.GetXaxis().SetTitle("true pT [MeV]") @@ -82,10 +92,18 @@ def main(path, particle): result.append(histo_count) return result + if __name__ == "__main__": output_file = ROOT.TFile.Open("average_layers.root", "recreate") - result = main("/storage_tmp/atlas/MVA2015/inputs_MC15/photon_grid_v5/data-output/*/*.root", "unconv") - result += main("/storage_tmp/atlas/MVA2015/inputs_MC15/photon_grid_v5/data-output/*/*.root", "conv") - result += main("/storage_tmp/atlas/MVA2015/inputs_MC15/electron_local_v5.1/data-output/mc15_electron.root", "electron") + result = main( + "/storage_tmp/atlas/MVA2015/inputs_MC15/" + "photon_grid_v5/data-output/*/*.root", + "unconv") + result += main("/storage_tmp/atlas/MVA2015/inputs_MC15/" + "photon_grid_v5/data-output/*/*.root", + "conv") + result += main("/storage_tmp/atlas/MVA2015/inputs_MC15/" + "electron_local_v5.1/data-output/mc15_electron.root", + "electron") for r in result: r.Write() diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/merge_scale_histograms.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/merge_scale_histograms.py index 12548deec608a8b987000f970a951dd89968082d..0e461885105ac35fee0526c238e4feb98b8417b9 100755 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/merge_scale_histograms.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/python/merge_scale_histograms.py @@ -3,21 +3,26 @@ # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +from array import array +import logging +import ROOT doc = """ This is an utility to merge histograms. The common case is when you have old scale factors for the whole calorimeter and new scale factor only for a subpart. """ -import ROOT ROOT.PyConfig.IgnoreCommandLineOptions = True -import logging logging.basicConfig(level=logging.INFO) -from array import array + def merge_histograms(old, new, merge_error=True): - print("old binning: " + ", ".join(("%.3f" % old.GetBinLowEdge(ibin)) for ibin in xrange(1, old.GetNbinsX() + 2))) - print("new binning: " + ", ".join(("%.3f" % new.GetBinLowEdge(ibin)) for ibin in xrange(1, new.GetNbinsX() + 2))) + print("old binning: " + ", ".join(("%.3f" % old.GetBinLowEdge(ibin)) + for ibin in range(1, old.GetNbinsX() + + 2))) + print("new binning: " + ", ".join(("%.3f" % new.GetBinLowEdge(ibin)) + for ibin in range(1, new.GetNbinsX() + + 2))) new_binning = [] new_values = [] @@ -25,24 +30,24 @@ def merge_histograms(old, new, merge_error=True): UNDERFLOW = 0 OVERFLOW = new.GetNbinsX() + 1 - for iold in xrange(1, old.GetNbinsX()): - l = old.GetBinLowEdge(iold) - r = l + old.GetBinWidth(iold) + for iold in range(1, old.GetNbinsX()): + low = old.GetBinLowEdge(iold) + r = low + old.GetBinWidth(iold) - il_new = new.FindFixBin(l) + il_new = new.FindFixBin(low) ir_new = new.FindFixBin(r) remainer = None if il_new == UNDERFLOW and ir_new == UNDERFLOW: - print("1. adding %.3f - %.3f from old" % (l, r)) - new_binning.append((l, r)) + print("1. adding %.3f - %.3f from old" % (low, r)) + new_binning.append((low, r)) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) elif il_new == UNDERFLOW and ir_new > UNDERFLOW: - if abs(new.GetBinLowEdge(1) - l) < 1E-100: + if abs(new.GetBinLowEdge(1) - low) < 1E-100: continue - new_binning.append((l, new.GetBinLowEdge(1))) + new_binning.append((low, new.GetBinLowEdge(1))) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) if ir_new == OVERFLOW: @@ -51,30 +56,24 @@ def merge_histograms(old, new, merge_error=True): break last_old = iold - for inew in xrange(1, new.GetNbinsX() + 1): - l = new.GetBinLowEdge(inew) - r = l + new.GetBinWidth(inew) - print("2. adding %.3f - %.3f from new" % (l, r)) - new_binning.append((l, r)) + for inew in range(1, new.GetNbinsX() + 1): + low = new.GetBinLowEdge(inew) + r = low + new.GetBinWidth(inew) + print("2. adding %.3f - %.3f from new" % (low, r)) + new_binning.append((low, r)) new_values.append(new.GetBinContent(inew)) new_errors.append(new.GetBinError(inew)) - """ - if remainer is not None: - print("3. adding %.3f - %.3f from old" % (new.GetBinLowEdge(new.GetNbinsX()), old.GetBinLowEdge(remainer) + old.GetBinWidth(remainer))) - new_binning.append((new.GetBinLowEdge(new.GetNbinsX()), old.GetBinLowEdge(remainer) + old.GetBinWidth(remainer))) - new_values.append(old.GetBinContent(remainer)) - new_errors.append(old.GetBinError(remainer)) - """ - for iold in xrange(last_old, old.GetNbinsX() + 1): - l = old.GetBinLowEdge(iold) - r = l + old.GetBinWidth(iold) - - il_new = new.FindFixBin(l) + + for iold in range(last_old, old.GetNbinsX() + 1): + low = old.GetBinLowEdge(iold) + r = low + old.GetBinWidth(iold) + + il_new = new.FindFixBin(low) ir_new = new.FindFixBin(r) if il_new == OVERFLOW and ir_new == OVERFLOW: - print("4. adding %.3f - %.3f from old" % (l, r)) - new_binning.append((l, r)) + print("4. adding %.3f - %.3f from old" % (low, r)) + new_binning.append((low, r)) new_values.append(old.GetBinContent(iold)) new_errors.append(old.GetBinError(iold)) elif il_new < OVERFLOW and ir_new == OVERFLOW: @@ -87,19 +86,20 @@ def merge_histograms(old, new, merge_error=True): print(new_binning) new_edges = array('f', [x[0] for x in new_binning] + [new_binning[-1][1]]) histo_type = type(new) - result = histo_type(new.GetName(), new.GetTitle(), len(new_edges) - 1, new_edges) + result = histo_type(new.GetName(), new.GetTitle(), + len(new_edges) - 1, new_edges) for i, (v, e) in enumerate(zip(new_values, new_errors), 1): result.SetBinContent(i, v) if merge_error: result.SetBinError(i, e) - print("merged binning: " + ", ".join(("%.3f" % result.GetBinLowEdge(ibin)) for ibin in xrange(1, result.GetNbinsX() + 1))) - + print("merged binning: " + ", ".join(("%.3f" % result.GetBinLowEdge(ibin)) + for ibin in range(1, result.GetNbinsX() + + 1))) return result - if __name__ == "__main__": try: ROOT.gROOT.ProcessLine(".x ~/atlasstyle/AtlasStyle.C") @@ -110,20 +110,13 @@ if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description=doc, - formatter_class=argparse.RawTextHelpFormatter, - epilog='''example (merge scales 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../RootCoreBin/download/ElectronPhotonFourMomentumCorrection/v7/egammaEnergyCorrectionData.root:Scales/es2015PRE/alphaZee_errStat ~/Scaricati/EnergyScaleFactors.root:centVal_alpha --title="2015 summer" --name alphaZee_errStat - -example (merge ct 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../RootCoreBin/download/ElectronPhotonFourMomentumCorrection/v7/egammaEnergyCorrectionData.root:Resolution/es2015PRE/ctZee_errStat ~/Scaricati/EnergyScaleFactors.root:centVal_c --title "2015 summer" --name ctZee_errStat - -example (merge scales sys 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../RootCoreBin/download/ElectronPhotonFourMomentumCorrection/v7/egammaEnergyCorrectionData.root:Scales/es2015PRE/alphaZee_errSyst ~/Scaricati/EnergyScaleFactors.root:totSyst_alpha --title="2015 summer" --name alphaZee_errSyst - -example (merge ct sys 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../RootCoreBin/download/ElectronPhotonFourMomentumCorrection/v7/egammaEnergyCorrectionData.root:Resolution/es2015PRE_res_improved/ctZee_errSyst ~/Scaricati/EnergyScaleFactors.root:totSyst_c --title="2015 summer" --name ctZee_errSyst -''') + formatter_class=argparse.RawTextHelpFormatter) parser.add_argument('histo_old') parser.add_argument('histo_new') parser.add_argument('--title', help='title of the new histogram') parser.add_argument('--name', help='name of the new histogram') - parser.add_argument('--recreate', help='create a new file', action='store_true') + parser.add_argument( + '--recreate', help='create a new file', action='store_true') parser.add_argument('--output-filename', default='output.root') args = parser.parse_args() @@ -139,7 +132,8 @@ example (merge ct sys 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../ if not histo_new: raise IOError("cannot find histogram %s" % args.histo_new) - logging.info("going to merge %s with %s", histo_old.GetName(), histo_new.GetName()) + logging.info("going to merge %s with %s", + histo_old.GetName(), histo_new.GetName()) histo_merged = merge_histograms(histo_old, histo_new) canvas = ROOT.TCanvas() @@ -159,11 +153,12 @@ example (merge ct sys 2015PRE + 2015 summer): ./merge_scale_histograms.py ../../ legend.SetBorderSize(0) legend.Draw() - fout = ROOT.TFile.Open(args.output_filename, "recreate" if args.recreate else "update") + fout = ROOT.TFile.Open(args.output_filename, + "recreate" if args.recreate else "update") if args.title is not None: histo_merged.SetTitle(args.title) name = args.name or histo_old.GetName() histo_merged.SetName(name) histo_merged.Write() - raw_input() + input("press a key")