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") diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test.py index ce598840a5fd4716323e1ba1d07406cbe4694a51..7d51513b54e2dd9c70ce1910b5e8639e9c19d5ba 100755 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test.py @@ -10,6 +10,7 @@ import random RUN2016 = 297730l RUN2015 = 252604l + def arange(xmin, xmax, delta): # just to don't inject dep from numpy result = [] @@ -35,7 +36,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): def setUpClass(cls): cls.event = ROOT.xAOD.TEvent() cls.factory = ROOT.EgammaFactory() - cls.tool_es2012c = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2012c") + cls.tool_es2012c = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2012c") def test_initialization(self): tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool") @@ -48,14 +50,19 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): self.assertTrue(tool.setProperty("ESModel", "es2010").isSuccess()) self.assertTrue(tool.initialize().isSuccess()) - def generator_kinematics(self, eta_range=None, e_range=None, phi_range=None): + def generator_kinematics(self, eta_range=None, + e_range=None, + phi_range=None): eta_range = eta_range or arange(-2.49, 2.49, 0.05) e_range = e_range or arange(5E3, 1000E3, 100E3) phi_range = phi_range or [0.] from itertools import product return product(e_range, eta_range, phi_range) - def generator_photon(self, eta_range=None, e_range=None, phi_range=None): + def generator_photon(self, + eta_range=None, + e_range=None, + phi_range=None): random.seed(10) factory = self.factory for e, eta, phi in self.generator_kinematics(eta_range, e_range, phi_range): @@ -66,21 +73,34 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): yield factory.create_unconverted_photon(eta, phi, e) factory.clear() - def generator_converted(self, eta_range=None, e_range=None, phi_range=None): + def generator_converted(self, eta_range=None, + e_range=None, + phi_range=None): factory = self.factory - for e, eta, phi in self.generator_kinematics(eta_range, e_range, phi_range): + for e, eta, phi in self.generator_kinematics(eta_range, + e_range, + phi_range): yield factory.create_converted_photon(eta, phi, e) factory.clear() - def generator_unconverted(self, eta_range=None, e_range=None, phi_range=None): + def generator_unconverted(self, eta_range=None, + e_range=None, + phi_range=None): factory = self.factory - for e, eta, phi in self.generator_kinematics(eta_range, e_range, phi_range): + for e, eta, phi in self.generator_kinematics(eta_range, + e_range, + phi_range): yield factory.create_unconverted_photon(eta, phi, e) factory.clear() - def generator_electron(self, eta_range=None, e_range=None, phi_range=None): + def generator_electron(self, + eta_range=None, + e_range=None, + phi_range=None): factory = self.factory - for e, eta, phi in self.generator_kinematics(eta_range, e_range, phi_range): + for e, eta, phi in self.generator_kinematics(eta_range, + e_range, + phi_range): yield factory.create_electron(eta, phi, e) factory.clear() @@ -96,7 +116,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): """ tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool") self.assertTrue(tool.setProperty("ESModel", "es2012c").isSuccess()) - self.assertTrue(tool.setProperty("int")("randomRunNumber", RUN2015).isSuccess()) + self.assertTrue(tool.setProperty("int")( + "randomRunNumber", RUN2015).isSuccess()) self.assertTrue(tool.setProperty("int")("doSmearing", 0).isSuccess()) tool.msg().setLevel(ROOT.MSG.WARNING) @@ -118,7 +139,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_2015PRE") self.assertTrue(tool.setProperty("ESModel", "es2015PRE").isSuccess()) self.assertTrue(tool.setProperty("int")("doSmearing", 0).isSuccess()) - self.assertTrue(tool.setProperty("decorrelationModel", "1NP_v1"). isSuccess()) + self.assertTrue(tool.setProperty( + "decorrelationModel", "1NP_v1"). isSuccess()) tool.msg().setLevel(ROOT.MSG.WARNING) self.assertTrue(tool.initialize().isSuccess()) ei = self.factory.create_eventinfo(True, 100000) # simulation @@ -141,12 +163,15 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): ei = self.factory.create_eventinfo(True, 100000) # simulation for el in self.generator_electron(): e_before = el.e() - e_cluster = el.caloCluster().energyBE(0) + el.caloCluster().energyBE(1) + el.caloCluster().energyBE(2) + el.caloCluster().energyBE(3) + e_cluster = el.caloCluster().energyBE(0) + el.caloCluster().energyBE(1) + \ + el.caloCluster().energyBE(2) + el.caloCluster().energyBE(3) tool.applyCorrection(el, ei) e_after = el.e() if e_before > 7E3: self.assertTrue(abs(e_before / e_after - 1) < 0.5, - msg="energy cluser/calibrated very different at eta = %f: %f/%f, e0=%.2f, e1=%.2f, e2=%.2f, e3=%.2f" % (el.eta(), e_cluster, e_after, el.caloCluster().energyBE(0), el.caloCluster().energyBE(1), el.caloCluster().energyBE(2), el.caloCluster().energyBE(3))) + msg="energy cluser/calibrated very different at eta = %f: %f/%f, e0=%.2f, e1=%.2f, e2=%.2f, e3=%.2f" % (el.eta(), + e_cluster, e_after, el.caloCluster().energyBE(0), el.caloCluster().energyBE(1), + el.caloCluster().energyBE(2), el.caloCluster().energyBE(3))) self.assertGreater(e_after, 0) def test_es2012XX(self): @@ -264,11 +289,14 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): import random import csv with open("testmva_%s_%s_%s.csv" % (esmodel, particle, "data" if isdata else "fullsim"), 'wb') as csvfile: - spamwriter = csv.writer(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL) + spamwriter = csv.writer( + csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL) if particle == "photon": - spamwriter.writerow(('eta', 'phi', 'e0', 'e1', 'e2', 'e3', 'e', 'rconv', 'zconv', 'output')) + spamwriter.writerow( + ('eta', 'phi', 'e0', 'e1', 'e2', 'e3', 'e', 'rconv', 'zconv', 'output')) elif particle == "electron": - spamwriter.writerow(('eta', 'phi', 'e0', 'e1', 'e2', 'e3', 'e', 'output')) + spamwriter.writerow( + ('eta', 'phi', 'e0', 'e1', 'e2', 'e3', 'e', 'output')) for eta in arange(-3.0, 3.0, 0.05): phi = 0 @@ -321,11 +349,13 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): import os basedir = os.path.dirname(os.path.abspath(__file__)) - filename = os.path.join(basedir, "testmva_%s_%s_%s.csv" % (esmodel, particle, "data" if isdata else "fullsim")) + filename = os.path.join(basedir, "testmva_%s_%s_%s.csv" % ( + esmodel, particle, "data" if isdata else "fullsim")) import csv with open(filename, 'rb') as csvfile: - reader = csv.reader(csvfile, delimiter=' ', quotechar='|', quoting=csv.QUOTE_MINIMAL) + reader = csv.reader(csvfile, delimiter=' ', + quotechar='|', quoting=csv.QUOTE_MINIMAL) next(reader, None) # skip header for row in reader: args = map(float, row[:-1]) @@ -333,7 +363,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): 'photon': self.factory.create_photon}[particle](*args) tool.applyCorrection(p, ei) calibrated_energy = p.e() - self.assertAlmostEqual(calibrated_energy, float(row[-1]), delta=0.1) + self.assertAlmostEqual( + calibrated_energy, float(row[-1]), delta=0.1) self.factory.clear() def test_list_syst(self): @@ -347,56 +378,84 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): """ tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool") tool.msg().setLevel(ROOT.MSG.WARNING) - self.assertTrue(tool.setProperty("int")("useMVACalibration", 0).isSuccess()) + self.assertTrue(tool.setProperty("int")( + "useMVACalibration", 0).isSuccess()) self.assertTrue(tool.setProperty("ESModel", model).isSuccess()) if decorrelation is not None: - self.assertTrue(tool.setProperty("decorrelationModel", decorrelation).isSuccess()) + self.assertTrue(tool.setProperty( + "decorrelationModel", decorrelation).isSuccess()) if decorrelation_scale is not None: - self.assertTrue(tool.setProperty("decorrelationModelScale", decorrelation_scale).isSuccess()) + self.assertTrue(tool.setProperty( + "decorrelationModelScale", decorrelation_scale).isSuccess()) if decorrelation_resolution is not None: - self.assertTrue(tool.setProperty("decorrelationModelResolution", decorrelation_resolution).isSuccess()) + self.assertTrue(tool.setProperty( + "decorrelationModelResolution", decorrelation_resolution).isSuccess()) tool.msg().setLevel(ROOT.MSG.WARNING) if success: - self.assertTrue(tool.initialize().isSuccess(), msg='cannot initialize tool with %s' % model) + self.assertTrue(tool.initialize().isSuccess(), + msg='cannot initialize tool with %s' % model) else: - self.assertFalse(tool.initialize().isSuccess(), msg='should not possible to initialize tool with %s' % model) + self.assertFalse(tool.initialize().isSuccess( + ), msg='should not possible to initialize tool with %s' % model) return sys_list = get_names_sys(tool) if type(allsyst) is int: - self.assertEqual(len(sys_list), allsyst, msg='actual (expected) number of sys %d(%d): %s' % (len(sys_list), allsyst, sys_list)) + self.assertEqual(len(sys_list), allsyst, msg='actual (expected) number of sys %d(%d): %s' % ( + len(sys_list), allsyst, sys_list)) else: self.assertItemsEqual(sys_list, allsyst) return sys_list - list_1NP_scale = ['EG_SCALE_ALL__1down', 'EG_SCALE_ALL__1up'] - list_1NP_resolution = ['EG_RESOLUTION_ALL__1down', 'EG_RESOLUTION_ALL__1up'] - list_FULL_resolution = ['EG_RESOLUTION_MATERIALCALO__1down', 'EG_RESOLUTION_MATERIALCALO__1up', - 'EG_RESOLUTION_MATERIALCRYO__1down', 'EG_RESOLUTION_MATERIALCRYO__1up', - 'EG_RESOLUTION_MATERIALGAP__1down', 'EG_RESOLUTION_MATERIALGAP__1up', - 'EG_RESOLUTION_MATERIALID__1down', 'EG_RESOLUTION_MATERIALID__1up', - 'EG_RESOLUTION_PILEUP__1down', 'EG_RESOLUTION_PILEUP__1up', - 'EG_RESOLUTION_SAMPLINGTERM__1down', 'EG_RESOLUTION_SAMPLINGTERM__1up', - 'EG_RESOLUTION_ZSMEARING__1down', 'EG_RESOLUTION_ZSMEARING__1up'] - list_1NPCOR_PLUS_UNCOR_scale = ['EG_SCALE_ALLCORR__1down', 'EG_SCALE_ALLCORR__1up', - 'EG_SCALE_E4SCINTILLATOR__1down', 'EG_SCALE_E4SCINTILLATOR__1up', - 'EG_SCALE_LARCALIB_EXTRA2015PRE__1down', 'EG_SCALE_LARCALIB_EXTRA2015PRE__1up', - 'EG_SCALE_LARTEMPERATURE_EXTRA2015PRE__1down', 'EG_SCALE_LARTEMPERATURE_EXTRA2015PRE__1up', - 'EG_SCALE_LARTEMPERATURE_EXTRA2016PRE__1down', 'EG_SCALE_LARTEMPERATURE_EXTRA2016PRE__1up'] - _test_list_syst("es2015PRE", "1NP_v1", None, None, list_1NP_scale + list_1NP_resolution) - _test_list_syst("es2012c", "1NP_v1", None, None, list_1NP_scale + list_1NP_resolution) - _test_list_syst("es2016PRE", None, "1NP_v1", "1NP_v1", list_1NP_scale + list_1NP_resolution) - _test_list_syst("es2016PRE", None, "1NP_v1", "FULL_v1", list_1NP_scale + list_FULL_resolution) - _test_list_syst("es2015PRE", "1NPCOR_PLUS_UNCOR", None, None, list_1NP_resolution + list_1NPCOR_PLUS_UNCOR_scale) - _test_list_syst("es2015PRE", "1NP_v1", "1NPCOR_PLUS_UNCOR", None, list_1NP_resolution + list_1NPCOR_PLUS_UNCOR_scale) - _test_list_syst("es2015c_summer", "1NP_v1", None, None, list_1NP_scale + list_1NP_resolution) + list_1NP_resolution = [ + 'EG_RESOLUTION_ALL__1down', 'EG_RESOLUTION_ALL__1up'] + list_FULL_resolution = ['EG_RESOLUTION_MATERIALCALO__1down', + 'EG_RESOLUTION_MATERIALCALO__1up', + 'EG_RESOLUTION_MATERIALCRYO__1down', + 'EG_RESOLUTION_MATERIALCRYO__1up', + 'EG_RESOLUTION_MATERIALGAP__1down', + 'EG_RESOLUTION_MATERIALGAP__1up', + 'EG_RESOLUTION_MATERIALID__1down', + 'EG_RESOLUTION_MATERIALID__1up', + 'EG_RESOLUTION_PILEUP__1down', + 'EG_RESOLUTION_PILEUP__1up', + 'EG_RESOLUTION_SAMPLINGTERM__1down', + 'EG_RESOLUTION_SAMPLINGTERM__1up', + 'EG_RESOLUTION_ZSMEARING__1down', + 'EG_RESOLUTION_ZSMEARING__1up'] + list_1NPCOR_PLUS_UNCOR_scale = ['EG_SCALE_ALLCORR__1down', + 'EG_SCALE_ALLCORR__1up', + 'EG_SCALE_E4SCINTILLATOR__1down', + 'EG_SCALE_E4SCINTILLATOR__1up', + 'EG_SCALE_LARCALIB_EXTRA2015PRE__1down', + 'EG_SCALE_LARCALIB_EXTRA2015PRE__1up', + 'EG_SCALE_LARTEMPERATURE_EXTRA2015PRE__1down', + 'EG_SCALE_LARTEMPERATURE_EXTRA2015PRE__1up', + 'EG_SCALE_LARTEMPERATURE_EXTRA2016PRE__1down', + 'EG_SCALE_LARTEMPERATURE_EXTRA2016PRE__1up'] + _test_list_syst("es2015PRE", "1NP_v1", None, None, + list_1NP_scale + list_1NP_resolution) + _test_list_syst("es2012c", "1NP_v1", None, None, + list_1NP_scale + list_1NP_resolution) + _test_list_syst("es2016PRE", None, "1NP_v1", "1NP_v1", + list_1NP_scale + list_1NP_resolution) + _test_list_syst("es2016PRE", None, "1NP_v1", "FULL_v1", + list_1NP_scale + list_FULL_resolution) + _test_list_syst("es2015PRE", "1NPCOR_PLUS_UNCOR", None, + None, list_1NP_resolution + list_1NPCOR_PLUS_UNCOR_scale) + _test_list_syst("es2015PRE", "1NP_v1", "1NPCOR_PLUS_UNCOR", + None, list_1NP_resolution + list_1NPCOR_PLUS_UNCOR_scale) + _test_list_syst("es2015c_summer", "1NP_v1", None, None, + list_1NP_scale + list_1NP_resolution) _test_list_syst("es2015PRE", "FULL_ETACORRELATED_v1", None, None, 58) _test_list_syst("es2012c", "FULL_ETACORRELATED_v1", None, None, 54) _test_list_syst("es2016PRE", "FULL_ETACORRELATED_v1", None, None, 62) - _test_list_syst("es2015c_summer", "FULL_ETACORRELATED_v1", None, None, 60) - _test_list_syst("es2016data_mc15c", "FULL_ETACORRELATED_v1", None, None, 68) + _test_list_syst("es2015c_summer", + "FULL_ETACORRELATED_v1", None, None, 60) + _test_list_syst("es2016data_mc15c", + "FULL_ETACORRELATED_v1", None, None, 68) _test_list_syst("es2012c", "FULL_v1", None, None, 148) _test_list_syst("es2012c", None, "FULL_v1", "FULL_v1", 148) _test_list_syst("es2015PRE", "FULL_v1", None, None, 158) @@ -404,9 +463,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): _test_list_syst("es2015PRE", None, None, None, 158) # these works, but generate FATALS, as expected - _test_list_syst("es2016PRE", "1NP_v1", "1NP_v1", "1NP_v1", [], success=False) - - + _test_list_syst("es2016PRE", "1NP_v1", "1NP_v1", + "1NP_v1", [], success=False) def test_same_smearing(self): tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool") @@ -462,7 +520,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): tool.applySystematicVariation(sys_set).ignore() e_zee_syst = tool.getEnergy(ph, ei) sys_set = ROOT.CP.SystematicSet() - sys_set.insert(ROOT.CP.SystematicVariation("EG_SCALE_LARCALIB__ETABIN0", 1.)) + sys_set.insert(ROOT.CP.SystematicVariation( + "EG_SCALE_LARCALIB__ETABIN0", 1.)) tool.applySystematicVariation(sys_set).ignore() e_larcalib__bin0 = tool.getEnergy(ph, ei) if abs(ph.eta()) < 2.47: @@ -478,9 +537,11 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): def test_intermodule_correction_working(self): tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool") self.assertTrue(tool.setProperty("ESModel", "es2012c").isSuccess()) - tool_no_correction = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_no_correction") + tool_no_correction = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_no_correction") tool_no_correction.setProperty("ESModel", "es2012c").ignore() - tool_no_correction.setProperty("int")("useIntermoduleCorrection", 0).ignore() + tool_no_correction.setProperty("int")( + "useIntermoduleCorrection", 0).ignore() self.assertTrue(tool.initialize().isSuccess()) self.assertTrue(tool_no_correction.initialize().isSuccess()) @@ -493,7 +554,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): energy.append(tool.getEnergy(ph, ei)) energy_no_correction.append(tool_no_correction.getEnergy(ph, ei)) - self.assertFalse(all([x == y for x, y in zip(energy, energy_no_correction)])) + self.assertFalse( + all([x == y for x, y in zip(energy, energy_no_correction)])) def test_idempotence(self): tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015PRE") @@ -516,23 +578,26 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): def test_1NP_vs_FULL_es2017(self): """ check that the 1NP model is the squared sum of the single systematics """ - tool_1NP = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2016data_mc15c_1NP") + tool_1NP = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2016data_mc15c_1NP") tool_1NP.setProperty("ESModel", "es2016data_mc15c").ignore() tool_1NP.setProperty("decorrelationModel", "1NP_v1").ignore() tool_1NP.setProperty("int")("randomRunNumber", RUN2015).ignore() - #tool_1NP.setProperty("int")("doSmearing", 0).ignore() # remove - #tool_1NP.msg().setLevel(ROOT.MSG.DEBUG) + # tool_1NP.setProperty("int")("doSmearing", 0).ignore() # remove + # tool_1NP.msg().setLevel(ROOT.MSG.DEBUG) tool_1NP.initialize().ignore() - tool_FULL = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2016data_mc15c_FULL") + tool_FULL = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2016data_mc15c_FULL") tool_FULL.setProperty("ESModel", "es2016data_mc15c").ignore() tool_FULL.setProperty("int")("randomRunNumber", RUN2015).ignore() # use ETACORRELATED to compare. FULL_v1 will differ (very small difference) since by default # FULL_v1 divide the ZEESTAT by the sqrt(#bins) - tool_FULL.setProperty("decorrelationModel", "FULL_ETACORRELATED_v1").ignore() - #tool_FULL.setProperty("int")("doSmearing", 0).ignore() # remove - #tool_FULL.msg().setLevel(ROOT.MSG.DEBUG) + tool_FULL.setProperty("decorrelationModel", + "FULL_ETACORRELATED_v1").ignore() + # tool_FULL.setProperty("int")("doSmearing", 0).ignore() # remove + # tool_FULL.msg().setLevel(ROOT.MSG.DEBUG) tool_FULL.initialize().ignore() ei = self.factory.create_eventinfo(True, 100000) # MC @@ -550,9 +615,11 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): all_syst_FULL = tool_FULL.recommendedSystematics() all_syst_1NP = tool_1NP.recommendedSystematics() - self.assertTrue('EG_SCALE_ALL__1up' in [s.name() for s in list(all_syst_1NP)]) + self.assertTrue('EG_SCALE_ALL__1up' in [ + s.name() for s in list(all_syst_1NP)]) - all_syst_FULL_scale_up = [s for s in all_syst_FULL if ('SCALE' in s.name() and '__1up' in s.name())] + all_syst_FULL_scale_up = [s for s in all_syst_FULL if ( + 'SCALE' in s.name() and '__1up' in s.name())] all_biases = [] for sys in all_syst_FULL_scale_up: @@ -561,12 +628,11 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): tool_FULL.applySystematicVariation(sys_set).ignore() e = tool_FULL.getEnergy(particle, ei) bias = e - e_full - #print "{:<40} {:10.5f} {:10.5f}".format(sys.name(), e, bias) + # print "{:<40} {:10.5f} {:10.5f}".format(sys.name(), e, bias) all_biases.append(bias) sum_quadrature = math.sqrt(sum([b * b for b in all_biases])) - #print "bias sum quadrature: ", sum_quadrature - + # print "bias sum quadrature: ", sum_quadrature bias_up, bias_down = None, None for sys in all_syst_1NP: @@ -581,12 +647,13 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): bias_up = bias elif sys.name() == 'EG_SCALE_ALL__1down': bias_down = bias - #print sys.name() + "\t" + str(bias) - + # print sys.name() + "\t" + str(bias) - self.assertAlmostEqual(sum_quadrature, bias_up, delta=1, msg="sum of errors not equal to 1NP (up) %.7f!=%.7f" % (sum_quadrature, bias_up)) + self.assertAlmostEqual(sum_quadrature, bias_up, delta=1, + msg="sum of errors not equal to 1NP (up) %.7f!=%.7f" % (sum_quadrature, bias_up)) # TODO: check why fails - #self.assertAlmostEqual(sum_quadrature, -bias_down, delta=1, msg="sum of errors not equal to 1NP (down) %.7f!=%.7f" % (sum_quadrature, bias_down)) + # self.assertAlmostEqual(sum_quadrature, -bias_down, delta=1, + # msg="sum of errors not equal to 1NP (down) %.7f!=%.7f" % (sum_quadrature, bias_down)) def test_1NP_100GeV_electron_es2015PRE(self): """ @@ -595,7 +662,8 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): """ tool = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015PRE") self.assertTrue(tool.setProperty("ESModel", "es2015PRE").isSuccess()) - self.assertTrue(tool.setProperty("decorrelationModel", "1NP_v1").isSuccess()) + self.assertTrue(tool.setProperty( + "decorrelationModel", "1NP_v1").isSuccess()) tool.msg().setLevel(ROOT.MSG.WARNING) tool.initialize().ignore() ei = self.factory.create_eventinfo(False, 100000) # data @@ -632,13 +700,15 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): check that es2015cPRE == es2015PRE for electron except for crack region [1.4-1.6] check that es2015cPRE == es2015PRE for photons """ - tool_es2015PRE = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015PRE") + tool_es2015PRE = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2015PRE") tool_es2015PRE.setProperty("ESModel", "es2015PRE").ignore() tool_es2015PRE.setProperty("decorrelationModel", "1NP_v1").ignore() tool_es2015PRE.setProperty("int")("doSmearing", 0).ignore() tool_es2015PRE.initialize().ignore() - tool_es2015cPRE = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015cPRE") + tool_es2015cPRE = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2015cPRE") tool_es2015cPRE.setProperty("ESModel", "es2015cPRE").ignore() tool_es2015cPRE.setProperty("decorrelationModel", "1NP_v1").ignore() tool_es2015cPRE.setProperty("int")("doSmearing", 0).ignore() @@ -651,7 +721,7 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): e_es2015PRE = tool_es2015PRE.getEnergy(el, ei) e_es2015cPRE = tool_es2015cPRE.getEnergy(el, ei) - #print el.eta(), el.e(), e_es2015PRE, e_es2015cPRE + # print el.eta(), el.e(), e_es2015PRE, e_es2015cPRE if 1.4 < abs(el.eta()) < 1.6: self.assertNotAlmostEqual(e_es2015PRE, e_es2015cPRE) else: @@ -670,16 +740,20 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): check that es2015c_summer == es2015cPRE for electrons != for photons """ - tool_es2015c_summer = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015c_summer") + tool_es2015c_summer = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2015c_summer") tool_es2015c_summer.setProperty("ESModel", "es2015c_summer").ignore() - tool_es2015c_summer.setProperty("decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() + tool_es2015c_summer.setProperty( + "decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() tool_es2015c_summer.setProperty("int")("doSmearing", 0).ignore() tool_es2015c_summer.msg().setLevel(ROOT.MSG.WARNING) tool_es2015c_summer.initialize().ignore() - tool_es2015cPRE = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015cPRE") + tool_es2015cPRE = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2015cPRE") tool_es2015cPRE.setProperty("ESModel", "es2015cPRE").ignore() - tool_es2015cPRE.setProperty("decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() + tool_es2015cPRE.setProperty( + "decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() tool_es2015cPRE.setProperty("int")("doSmearing", 0).ignore() tool_es2015cPRE.msg().setLevel(ROOT.MSG.WARNING) tool_es2015cPRE.initialize().ignore() @@ -701,15 +775,18 @@ class TestEgammaCalibrationAndSmearingTool(unittest.TestCase): e_es2015cPRE = tool_es2015cPRE.getEnergy(ph, ei) self.assertGreater(e_es2015cPRE, 0) self.assertNotAlmostEqual(e_es2015c_summer, e_es2015cPRE, - msg="e_es2015c_summer == e_es2015cPRE == %.2f at eta=%.2f for photons" % (e_es2015c_summer, ph.eta())) + msg="e_es2015c_summer == e_es2015cPRE == %.2f at eta=%.2f for photons" % (e_es2015c_summer, + ph.eta())) def test_es2015c_summer_data(self): """ check that energy > 0 """ - tool_es2015c_summer = ROOT.CP.EgammaCalibrationAndSmearingTool("tool_es2015c_summer") + tool_es2015c_summer = ROOT.CP.EgammaCalibrationAndSmearingTool( + "tool_es2015c_summer") tool_es2015c_summer.setProperty("ESModel", "es2015c_summer").ignore() - tool_es2015c_summer.setProperty("decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() + tool_es2015c_summer.setProperty( + "decorrelationModel", "1NPCOR_PLUS_UNCOR").ignore() tool_es2015c_summer.msg().setLevel(ROOT.MSG.WARNING) tool_es2015c_summer.initialize().ignore() @@ -728,5 +805,5 @@ if __name__ == '__main__': ROOT.gROOT.ProcessLine(".x $ROOTCOREDIR/scripts/load_packages.C") # from ROOT import EgammaCalibPeriodRunNumbersExample - #ROOT.xAOD.TReturnCode.enableFailure() + # ROOT.xAOD.TReturnCode.enableFailure() unittest.main() diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_factory.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_factory.py index ff2ae3998b8014c9efb8f7405a883bd8e030e215..644e69fe4f8ce8419ddff62f334f96649fc683cc 100755 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_factory.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_factory.py @@ -8,7 +8,7 @@ import unittest class TestEgammaFactory(unittest.TestCase): def test_eventinfo(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() ei1 = factory.create_eventinfo(True, 100000) self.assertTrue(ei1.eventType(ROOT.xAOD.EventInfo.IS_SIMULATION)) @@ -16,7 +16,7 @@ class TestEgammaFactory(unittest.TestCase): self.assertFalse(ei2.eventType(ROOT.xAOD.EventInfo.IS_SIMULATION)) def test_unconverted(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() runnumber = 10000 @@ -32,11 +32,13 @@ class TestEgammaFactory(unittest.TestCase): self.assertAlmostEqual(ph.caloCluster().eta(), eta) self.assertAlmostEqual(ph.phi(), phi) self.assertAlmostEqual(ph.caloCluster().phi(), phi) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("etaCalo"), eta) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("phiCalo"), phi) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("etaCalo"), eta) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("phiCalo"), phi) self.assertAlmostEqual(ph.e(), e, delta=0.01) self.assertEqual(ph.caloCluster().e(), e) - for i in xrange(3): + for i in range(3): self.assertGreater(ph.caloCluster().energyBE(i), 0) def test_converted(self): @@ -56,15 +58,17 @@ class TestEgammaFactory(unittest.TestCase): self.assertAlmostEqual(ph.caloCluster().eta(), eta) self.assertAlmostEqual(ph.phi(), phi) self.assertAlmostEqual(ph.caloCluster().phi(), phi) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("etaCalo"), eta) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("phiCalo"), phi) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("etaCalo"), eta) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("phiCalo"), phi) self.assertAlmostEqual(ph.e(), e, delta=0.01) self.assertEqual(ph.caloCluster().e(), e) - for i in xrange(3): + for i in range(3): self.assertGreater(ph.caloCluster().energyBE(i), 0) def test_photon(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() runnumber = 10000 @@ -89,8 +93,10 @@ class TestEgammaFactory(unittest.TestCase): self.assertAlmostEqual(ph.caloCluster().eta(), eta) self.assertAlmostEqual(ph.phi(), phi) self.assertAlmostEqual(ph.caloCluster().phi(), phi) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("etaCalo"), eta) - self.assertAlmostEqual(ph.caloCluster().auxdata("float")("phiCalo"), phi) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("etaCalo"), eta) + self.assertAlmostEqual( + ph.caloCluster().auxdata("float")("phiCalo"), phi) self.assertAlmostEqual(ph.e(), e, delta=0.01) self.assertEqual(ph.caloCluster().e(), e) @@ -101,7 +107,7 @@ class TestEgammaFactory(unittest.TestCase): tool.applyCorrection(ph, ei) def test_electron(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() runnumber = 10000 @@ -124,14 +130,16 @@ class TestEgammaFactory(unittest.TestCase): self.assertAlmostEqual(el.caloCluster().eta(), eta) self.assertAlmostEqual(el.phi(), phi) self.assertAlmostEqual(el.caloCluster().phi(), phi) - self.assertAlmostEqual(el.caloCluster().auxdata("float")("etaCalo"), eta) - self.assertAlmostEqual(el.caloCluster().auxdata("float")("phiCalo"), phi) + self.assertAlmostEqual( + el.caloCluster().auxdata("float")("etaCalo"), eta) + self.assertAlmostEqual( + el.caloCluster().auxdata("float")("phiCalo"), phi) self.assertAlmostEqual(el.e(), e, delta=0.01) self.assertEqual(el.caloCluster().e(), e) self.assertAlmostEqual(el.trackParticle().eta(), eta, delta=0.00001) def test_clean(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() runnumber = 10000 @@ -147,7 +155,7 @@ class TestEgammaFactory(unittest.TestCase): self.assertAlmostEqual(el.eta(), 1) def test_decoration(self): - event = ROOT.xAOD.TEvent() + ROOT.xAOD.TEvent() factory = ROOT.EgammaFactory() runnumber = 10000 @@ -158,7 +166,6 @@ class TestEgammaFactory(unittest.TestCase): el.auxdecor("double")("mydeco") el.auxdataConst("double")("mydeco") - #factory.clear() el = factory.create_electron(1, 2, 3, 4, 5, 6, 7) el.auxdataConst("double")("mydeco") diff --git a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_resolution.py b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_resolution.py index b4a330f56aceeb8ef213a825043a182c0c7a807f..72b99d4cbf6ffdc8e8b839d3f65b3d2bbdcb2351 100755 --- a/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_resolution.py +++ b/PhysicsAnalysis/ElectronPhotonID/ElectronPhotonFourMomentumCorrection/test/ut_test_resolution.py @@ -31,16 +31,16 @@ class TestEgammaResolution(unittest.TestCase): cls.factory = ROOT.EgammaFactory() def loop_combination(self, eta_min=-2.5, eta_max=2.5, eta_step=0.05, - energy_min=0, energy_max=200E3, energy_step=10E3, - ptypes=(0, 1, 2, 3), - resol_types=(0, 1, 2)): + energy_min=0, energy_max=200E3, energy_step=10E3, + ptypes=(0, 1, 2, 3), + resol_types=(0, 1, 2)): for ptype in ptypes: for resol_type in resol_types: for eta in arange(eta_min, eta_max, eta_step): for energy in arange(energy_min, energy_max, energy_step): yield ptype, energy, eta, resol_type - #@unittest.expectedFailure # return nan in the crack (1.37-1.52), inf at >= 2.4 + # @unittest.expectedFailure # return nan in the crack (1.37-1.52), inf at >= 2.4 def test_resolution_positive(self): for ptype, energy, eta, resol_type in self.loop_combination(): for tool in self.tool_run1, self.tool_run2: @@ -49,7 +49,8 @@ class TestEgammaResolution(unittest.TestCase): if ptype == 3: continue resolution = tool.getResolution(ptype, energy, eta, resol_type) - self.assertFalse(math.isnan(resolution), msg="resolution is nan for eta=%f" % eta) + self.assertFalse(math.isnan(resolution), + msg="resolution is nan for eta=%f" % eta) self.assertGreater(resolution, 0) def create_response_run1(self): @@ -57,9 +58,13 @@ class TestEgammaResolution(unittest.TestCase): this should not be executed, to don't invalidate reference file created on 2015-04-27 with ElectronPhotonFourMomentumCorrection-00-01-23 """ - fout = ROOT.TFile("$ROOTCOREBIN/data/ElectronPhotonFourMomentumCorrection/test_resolution_nonregression_run1_data.root", "recreate") - tree = ROOT.TTree("test_resolution_nonregression_data_run1", "test_resolution_nonregression_data") - data_ptype, data_resol_type, data_eta, data_energy, data_resolution = map(lambda t: array(t, [0]), ('i', 'i', 'f', 'f', 'f')) + fout = ROOT.TFile( + "$ROOTCOREBIN/data/ElectronPhotonFourMomentumCorrection/test_resolution_nonregression_run1_data.root", + "recreate") + tree = ROOT.TTree("test_resolution_nonregression_data_run1", + "test_resolution_nonregression_data") + data_ptype, data_resol_type, data_eta, data_energy, data_resolution = map( + lambda t: array(t, [0]), ('i', 'i', 'f', 'f', 'f')) tree.Branch("ptype", data_ptype, "ptype/I") tree.Branch("resol_type", data_resol_type, "resol_type/I") @@ -68,7 +73,8 @@ class TestEgammaResolution(unittest.TestCase): tree.Branch("resolution", data_resolution, "resolution/F") for ptype, energy, eta, resol_type in self.loop_combination_run1(): - resolution = self.tool_run1.getResolution(ptype, energy, eta, resol_type) + resolution = self.tool_run1.getResolution( + ptype, energy, eta, resol_type) data_ptype[0] = ptype data_energy[0] = energy data_eta[0] = eta @@ -86,18 +92,23 @@ class TestEgammaResolution(unittest.TestCase): for eta in arange(0, 1.37, 0.1): for e in arange(10E3, 1000E3, 100E3): for t in 0, 1, 2: - resolution_run1 = self.tool_run1.getResolution(ptype, e, eta, t) - resolution_run2 = tool_run2.getResolution(ptype, e, eta, t) - self.assertNotAlmostEqual(resolution_run1, resolution_run2, - msg="resolution should be different for particle=%d, eta=%f, e=%d, rtype=%d" % (ptype, eta, e, t)) + resolution_run1 = self.tool_run1.getResolution( + ptype, e, eta, t) + resolution_run2 = tool_run2.getResolution( + ptype, e, eta, t) + self.assertNotAlmostEqual(resolution_run1, + resolution_run2, + msg="resolution should be different for particle=%d, eta=%f, e=%d, rtype=%d" % (ptype, + eta, e, t)) def test_nonregression_run1(self): from PathResolver import PathResolver - rootfile = PathResolver.FindCalibFile("ElectronPhotonFourMomentumCorrection/v8/test_resolution_nonregression_run1_data.root") + rootfile = PathResolver.FindCalibFile( + "ElectronPhotonFourMomentumCorrection/v8/test_resolution_nonregression_run1_data.root") f = ROOT.TFile(rootfile) tree = f.Get("test_resolution_nonregression_data_run1") - for ievent in xrange(tree.GetEntries()): + for ievent in range(tree.GetEntries()): tree.GetEntry(ievent) resolution = self.tool_run1.getResolution(tree.ptype, @@ -105,12 +116,20 @@ class TestEgammaResolution(unittest.TestCase): tree.resol_type) expected_response = tree.resolution if math.isnan(resolution): - self.assertTrue(math.isnan(expected_response), msg="resolution is nan, but expected resonse is %f" % expected_response) + self.assertTrue(math.isnan( + expected_response), msg="resolution is nan, but expected resonse is %f" % expected_response) elif math.isnan(expected_response): - self.assertTrue(math.isnan(resolution), msg="expected response is nan, but resolution is %f" % resolution) + self.assertTrue(math.isnan( + resolution), msg="expected response is nan, but resolution is %f" % resolution) else: self.assertAlmostEqual(resolution, expected_response, - msg="resolution mismatch %f != %f at eta=%f, energy=%f, ptype=%d, resol_type=%d" % (resolution, expected_response, tree.eta, tree.energy, tree.ptype, tree.resol_type)) + msg="resolution mismatch %f != %f at eta=%f, energy=%f, ptype=%d, resol_type=%d" % ( + resolution, + expected_response, + tree.eta, + tree.energy, + tree.ptype, + tree.resol_type)) @unittest.skip("CHECK") def test_resolution_interface(self): @@ -125,22 +144,29 @@ class TestEgammaResolution(unittest.TestCase): for energy in energies: if particle == 'unconverted': p = self.factory.create_photon(0., 0.1, energy, 0) - r1 = tool.resolution(energy, 0., 0., ROOT.PATCore.ParticleType.UnconvertedPhoton) + r1 = tool.resolution( + energy, 0., 0., + ROOT.PATCore.ParticleType.UnconvertedPhoton) r2 = tool.getResolution(p) self.assertAlmostEqual(r1, r2) elif particle == 'converted': p = self.factory.create_photon(0., 0.1, energy, 100) - r1 = tool.resolution(energy, 0., 0., ROOT.PATCore.ParticleType.ConvertedPhoton) + r1 = tool.resolution( + energy, 0., 0., + ROOT.PATCore.ParticleType.ConvertedPhoton) r2 = tool.getResolution(p) self.assertAlmostEqual(r1, r2) elif particle == 'electron': p = self.factory.create_electron(0., 0.1, energy) - r1 = tool.resolution(energy, 0., 0., ROOT.PATCore.ParticleType.Electron) + r1 = tool.resolution( + energy, 0., 0., + ROOT.PATCore.ParticleType.Electron) r2 = tool.getResolution(p) self.assertAlmostEqual(r1, r2) else: raise ValueError() + if __name__ == '__main__': ROOT.gROOT.ProcessLine(".x $ROOTCOREDIR/scripts/load_packages.C") unittest.main()