diff --git a/Hlt/RecoConf/options/hlt1_reco_pvchecker.py b/Hlt/RecoConf/options/hlt1_reco_pvchecker.py index 3064151a39ddc78a748364d5ea5498ef24d505b4..29bcedf2c455f305bd5ca4a9c305e2dac7278ca7 100644 --- a/Hlt/RecoConf/options/hlt1_reco_pvchecker.py +++ b/Hlt/RecoConf/options/hlt1_reco_pvchecker.py @@ -12,8 +12,6 @@ from Moore import options, run_reconstruction from Moore.config import Reconstruction from RecoConf.hlt1_tracking import require_gec, make_hlt1_tracks, make_pvs from RecoConf.mc_checking import get_pv_checkers -from Configurables import ApplicationMgr -from Configurables import NTupleSvc def hlt1_reco_pvchecker(): @@ -27,10 +25,5 @@ def hlt1_reco_pvchecker(): return Reconstruction('PVperformance', data, [require_gec()]) +options.ntuple_file = 'Hlt1_PVperformance.root' run_reconstruction(options, hlt1_reco_pvchecker) - -NTupleSvc().Output += [ - "FILE1 DATAFILE='Hlt1_PVperformance.root' TYPE='ROOT' OPT='NEW'" -] -ApplicationMgr().ExtSvc += [NTupleSvc()] -ApplicationMgr().HistogramPersistency = "ROOT" diff --git a/Hlt/RecoConf/python/RecoConf/mc_checking.py b/Hlt/RecoConf/python/RecoConf/mc_checking.py index dca8dc26cabb5b8b7afcc154e812141d276208e0..036a82a3098e2eaaa1168199ec9a5fb29f71c40d 100644 --- a/Hlt/RecoConf/python/RecoConf/mc_checking.py +++ b/Hlt/RecoConf/python/RecoConf/mc_checking.py @@ -25,7 +25,7 @@ from PyConf.Algorithms import ( PrLHCbID2MCParticle, PrLHCbID2MCParticleVPUTFTMU, PrTrackAssociator, PrTrackChecker, PrUTHitChecker, TrackListRefiner, TrackResChecker, PrMultiplicityChecker, TrackIPResolutionCheckerNT, MCParticle2MCHitAlg, - PrTrackerDumper, PVDumper, + PrTrackerDumper, PVDumper, PrimaryVertexChecker, LHCb__Converters__RecVertex__v2__fromVectorLHCbRecVertices as FromVectorLHCbRecVertex) @@ -396,6 +396,28 @@ def get_best_tracks_checkers( return efficiency_checkers +@configurable +def get_pv_checkers(pvs, tracks, produce_ntuple=False): + + assert isinstance( + pvs, DataHandle), "Please provide reconstructed primary verticies" + + pv_checkers = [] + + pvchecker = PrimaryVertexChecker( + produceNtuple=produce_ntuple, + inputVerticesName=pvs, + inputTracksName=tracks["v1"], + MCVertexInput=mc_unpackers()["MCVertices"], + MCParticleInput=mc_unpackers()["MCParticles"], + MCHeaderLocation=make_data_with_FetchDataFromFile("/Event/MC/Header"), + MCPropertyInput=make_data_with_FetchDataFromFile( + "/Event/MC/TrackInfo")) + + pv_checkers.append(pvchecker) + return pv_checkers + + def make_track_filter(InputTracks, code): selector = LoKiTrackSelector(Code=code) filtered_tracks = TrackListRefiner( diff --git a/Hlt/RecoConf/python/RecoConf/standalone.py b/Hlt/RecoConf/python/RecoConf/standalone.py index 1761dec8e72f0fe9e39477deb1865a8d316c62f4..fb1c4723a5a074a9f34453ffe126baa55dcc5197 100644 --- a/Hlt/RecoConf/python/RecoConf/standalone.py +++ b/Hlt/RecoConf/python/RecoConf/standalone.py @@ -17,7 +17,9 @@ from .hlt1_tracking import ( from RecoConf.hlt1_muonmatch import make_tracks_with_muonmatch_ipcut from .hlt1_muonid import make_muon_id, make_tracks_with_muon_id from .hlt2_muonid import make_muon_ids as make_muon_id_hlt2 -from .hlt2_tracking import make_hlt2_tracks, make_TrackBestTrackCreator_tracks, get_default_tracks_for_calo, get_default_out_track_types_for_light_reco +from .hlt2_tracking import ( + make_hlt2_tracks, make_TrackBestTrackCreator_tracks, + get_default_tracks_for_calo, get_default_out_track_types_for_light_reco) from .calorimeter_reconstruction import (make_calo, make_calo_resolution_gamma, make_calo_resolution_pi0) from .calorimeter_mc_checking import ( @@ -36,10 +38,11 @@ from PyConf.Algorithms import ( RawBankSizeMonitor, ) from .mc_checking import (get_track_checkers, get_fitted_tracks_checkers, - get_best_tracks_checkers, + get_pv_checkers, get_best_tracks_checkers, get_track_checkers_multiplicity) from .reconstruction_objects import reconstruction -from .protoparticles import make_charged_protoparticles, make_neutral_protoparticles +from .protoparticles import (make_charged_protoparticles, + make_neutral_protoparticles) from .calo_data_monitoring import monitor_calo_clusters from .track_data_monitoring import monitor_tracking @@ -82,6 +85,7 @@ def standalone_hlt1_reco(do_mc_checking=False): "Forward": hlt1_tracks["Forward"], } data += get_track_checkers(types_and_locations_for_checkers) + data += get_pv_checkers(pvs, hlt1_tracks["Velo"]) return Reconstruction('hlt1_reco', data, reco_prefilters()) diff --git a/Hlt/RecoConf/scripts/PrimaryVertexCheckerBasic.py b/Hlt/RecoConf/scripts/PrimaryVertexCheckerBasic.py new file mode 100644 index 0000000000000000000000000000000000000000..0ae38bc39ea67b142e6895841d265e00e8305b60 --- /dev/null +++ b/Hlt/RecoConf/scripts/PrimaryVertexCheckerBasic.py @@ -0,0 +1,167 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#!/usr/bin/python + +# The script for plotting PV efficinecy as the function +# of various distributions: nTracks, z, r. +# As input the NTuple created by hlt1_reco_pvchecker.py +# is needed. +# +# The efficency is calculated usig TGraphAsymmErrors +# and Bayesian error bars +# +# author: Agnieszka Dziurda (agnieszka.dziurda@cern.ch) +# date: 02/2020 +# +# Example of usage: +# ../../../run python PrimaryVertexCheckerPull.py +# --file file1.root file2.root --label name1 name2 +# + +import os, sys +import argparse + +from ROOT import gROOT, TLegend + +parser = argparse.ArgumentParser() +parser.add_argument( + '--file', dest='fileName', default="", nargs='+', help='filename to plot') +parser.add_argument( + '--label', dest='label', default="", nargs='+', help='labels for files') +parser.add_argument( + '--tree', + dest='treeName', + default="", + nargs='+', + help='tree name to plot', +) +parser.add_argument( + '--smog', + dest='smog', + default=False, + action='store_true', + help='set true for SMOG') +parser.add_argument( + '--multi', + dest='multi', + default=False, + action='store_true', + help='add multiplicity plots') +parser.add_argument( + '--isol', + dest='isol', + default=False, + action='store_true', + help='add isolated/closed plots') +parser.add_argument( + '--dist', + dest='dist', + default=False, + action='store_true', + help='plot distributions in the canvas') +parser.add_argument( + '--prefix', + dest='prefix', + default="pv_basic", + help='prefix for the plot name', +) +parser.add_argument( + '--dir', + dest='directory', + default=os.getcwd(), + help='tree name to plot', +) + +parser.add_argument( + '--offset', + dest='offset', + default=0, + help='offset for plot colors', +) + + +def get_labels(number_of_files): + label = [] + for i in range(0, number_of_files): + label.append("PV Checker {number}".format(number=str(i + 1))) + return label + + +if __name__ == '__main__': + args = parser.parse_args() + path = args.directory + "/utils/" + sys.path.append(os.path.abspath(path)) + offset = int(args.offset) + gROOT.SetBatch() + + from pvutils import get_files, get_trees + from pvutils import get_global, plot_comparison + + from pvconfig import get_variable_ranges + from pvconfig import set_legend_simple, get_style, get_categories + + ranges = get_variable_ranges(args.smog) + style = get_style() + + cat = get_categories(args.multi, args.isol, args.smog) + label = args.label + if args.label == "": + label = get_labels(len(args.fileName)) + + tr = {} + tf = {} + tf = get_files(tf, label, args.fileName) + tr = get_trees(tf, tr, label, args.treeName, True) + + hist_trueMCPVs = {} + hist_visMCPVs = {} + hist_zMC = {} + hist_xMC = {} + hist_yMC = {} + hist_part = {} + + norm = True #to-do + hist_trueMCPVs = get_global(hist_trueMCPVs, tr, "mtruemcpv", "All MC PVs", + "Candidates Normalized", style, ranges, cat, + label, offset, True) + hist_visMCPVs = get_global( + hist_visMCPVs, tr, "nmcpv", "Reconstructible MC PVs}", + "Candidates Normalized", style, ranges, cat, label, offset, True) + hist_zMC = get_global(hist_zMC, tr, "zMC", "PV z [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset, True) + hist_xMC = get_global(hist_xMC, tr, "xMC", "PV x [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset, True) + hist_yMC = get_global(hist_yMC, tr, "yMC", "PV y [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset, True) + hist_part = get_global(hist_part, tr, "nrectrmc", "MC Particles in MC PVs", + "Candidates Normalized", style, ranges, cat, label, + offset, True) + + legend = TLegend(0.15, 0.86, 0.88, 0.98) + legend = set_legend_simple(legend, label, hist_trueMCPVs) + + plot_comparison(hist_trueMCPVs, args.prefix, "mtruemcpv", cat, label, + style, norm, offset, True, legend) + plot_comparison(hist_visMCPVs, args.prefix, "nmcpv", cat, label, style, + norm, offset, True, legend) + + plot_comparison(hist_zMC, args.prefix, "zMC", cat, label, style, norm, + offset, True, legend) + plot_comparison(hist_xMC, args.prefix, "xMC", cat, label, style, norm, + offset, True, legend) + plot_comparison(hist_yMC, args.prefix, "yMC", cat, label, style, norm, + offset, True, legend) + + plot_comparison(hist_part, args.prefix, "nrectrmc", cat, label, style, + norm, offset, True, legend) diff --git a/Hlt/RecoConf/scripts/PrimaryVertexCheckerEfficiency.py b/Hlt/RecoConf/scripts/PrimaryVertexCheckerEfficiency.py new file mode 100644 index 0000000000000000000000000000000000000000..612d9176001b44daec814a13f9b13257db25a02c --- /dev/null +++ b/Hlt/RecoConf/scripts/PrimaryVertexCheckerEfficiency.py @@ -0,0 +1,156 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#!/usr/bin/python + +# The script for plotting PV efficinecy as the function +# of various distributions: nTracks, z, r. +# As input the NTuple created by hlt1_reco_pvchecker.py +# is needed. +# +# The efficency is calculated usig TGraphAsymmErrors +# and Bayesian error bars +# +# author: Agnieszka Dziurda (agnieszka.dziurda@cern.ch) +# date: 02/2020 +# +# Example of usage: +# ../../../run python PrimaryVertexCheckerEfficiency.py +# --file file1.root file2.root --label name1 name2 --dist +# + +import os, sys +import argparse + +from ROOT import gROOT, TLegend + +parser = argparse.ArgumentParser() +parser.add_argument( + '--file', dest='fileName', default="", nargs='+', help='filename to plot') +parser.add_argument( + '--label', dest='label', default="", nargs='+', help='labels for files') +parser.add_argument( + '--tree', + dest='treeName', + default="", + nargs='+', + help='tree name to plot', +) +parser.add_argument( + '--smog', + dest='smog', + default=False, + action='store_true', + help='set true for SMOG') +parser.add_argument( + '--multi', + dest='multi', + default=False, + action='store_true', + help='add multiplicity plots') +parser.add_argument( + '--isol', + dest='isol', + default=False, + action='store_true', + help='add isolated/closed plots') +parser.add_argument( + '--dist', + dest='dist', + default=False, + action='store_true', + help='plot distributions in the canvas') +parser.add_argument( + '--prefix', + dest='prefix', + default="pv_eff", + help='prefix for the plot name', +) +parser.add_argument( + '--dir', + dest='directory', + default=os.getcwd(), + help='tree name to plot', +) + +parser.add_argument( + '--offset', + dest='offset', + default=0, + help='offset for plot colors', +) + + +def get_labels(number_of_files): + label = [] + for i in range(0, number_of_files): + label.append("PV Checker {number}".format(number=str(i + 1))) + return label + + +if __name__ == '__main__': + args = parser.parse_args() + path = args.directory + "/utils/" + sys.path.append(os.path.abspath(path)) + offset = int(args.offset) + + gROOT.SetBatch() + + from pvutils import get_files, get_trees, get_eff, plot + + from pvconfig import get_variable_ranges + from pvconfig import get_style, get_categories + from pvconfig import set_legend + + ranges = get_variable_ranges(args.smog) + style = get_style() + + label = args.label + if args.label == "": + label = get_labels(len(args.fileName)) + + tr = {} + tf = {} + tf = get_files(tf, label, args.fileName) + tr = get_trees(tf, tr, label, args.treeName, True) + + eff = {} + eff["tracks"] = {} + eff["z"] = {} + eff["r"] = {} + + hist = {} + hist["tracks"] = {} + hist["z"] = {} + hist["r"] = {} + + cat = get_categories(args.multi, args.isol, args.smog) + + eff["tracks"], hist["tracks"] = get_eff(eff["tracks"], hist["tracks"], tr, + "nrectrmc", style, ranges, cat, + label, offset) + eff["z"], hist["z"] = get_eff(eff["z"], hist["z"], tr, "zMC", style, + ranges, cat, label, offset) + eff["r"], hist["r"] = get_eff(eff["r"], hist["r"], tr, "rMC", style, + ranges, cat, label, offset) + + if args.dist: + legend = TLegend(0.15, 0.82, 0.88, 0.98) + else: + legend = TLegend(0.15, 0.86, 0.88, 0.98) + legend = set_legend(legend, label, eff["tracks"], "eff", hist["z"], + args.dist) + + plot(eff["tracks"], "eff", args.prefix, "ntracks", cat, label, legend, + hist["tracks"], args.dist, (0.9, 0.65)) + plot(eff["z"], "eff", args.prefix, "z", cat, label, legend, hist["z"], + args.dist, (0.9, 0.65)) + plot(eff["r"], "eff", args.prefix, "r", cat, label, legend, hist["r"], + args.dist, (0.9, 0.3)) diff --git a/Hlt/RecoConf/scripts/PrimaryVertexCheckerPull.py b/Hlt/RecoConf/scripts/PrimaryVertexCheckerPull.py new file mode 100644 index 0000000000000000000000000000000000000000..80c7785c2699d4292948fd715fd432847a94e76b --- /dev/null +++ b/Hlt/RecoConf/scripts/PrimaryVertexCheckerPull.py @@ -0,0 +1,296 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#!/usr/bin/python + +# The script for plotting PV efficinecy as the function +# of various distributions: nTracks, z, r. +# As input the NTuple created by hlt1_reco_pvchecker.py +# is needed. +# +# The efficency is calculated usig TGraphAsymmErrors +# and Bayesian error bars +# +# author: Agnieszka Dziurda (agnieszka.dziurda@cern.ch) +# date: 02/2020 +# +# Example of usage: +# ../../../run python PrimaryVertexCheckerPull.py +# --file file1.root file2.root --label name1 name2 +# + +import os, sys +import argparse + +from ROOT import gROOT, TLegend + +parser = argparse.ArgumentParser() +parser.add_argument( + '--file', dest='fileName', default="", nargs='+', help='filename to plot') +parser.add_argument( + '--label', dest='label', default="", nargs='+', help='labels for files') +parser.add_argument( + '--tree', + dest='treeName', + default="", + nargs='+', + help='tree name to plot', +) +parser.add_argument( + '--smog', + dest='smog', + default=False, + action='store_true', + help='set true for SMOG') +parser.add_argument( + '--multi', + dest='multi', + default=False, + action='store_true', + help='add multiplicity plots') +parser.add_argument( + '--isol', + dest='isol', + default=False, + action='store_true', + help='add isolated/closed plots') +parser.add_argument( + '--dist', + dest='dist', + default=False, + action='store_true', + help='plot distributions in the canvas') +parser.add_argument( + '--prefix', + dest='prefix', + default="pv_pull", + help='prefix for the plot name', +) +parser.add_argument( + '--dir', + dest='directory', + default=os.getcwd(), + help='tree name to plot', +) + +parser.add_argument( + '--offset', + dest='offset', + default=0, + help='offset for plot colors', +) + + +def get_labels(number_of_files): + label = [] + for i in range(0, number_of_files): + label.append("PV Checker {number}".format(number=str(i + 1))) + return label + + +if __name__ == '__main__': + args = parser.parse_args() + path = args.directory + "/utils/" + sys.path.append(os.path.abspath(path)) + offset = int(args.offset) + gROOT.SetBatch() + + from pvutils import get_files, get_trees + from pvutils import get_global, plot_comparison + from pvutils import get_dependence + from pvutils import plot + + from pvconfig import get_variable_ranges + from pvconfig import set_legend, get_style, get_categories + + ranges = get_variable_ranges(args.smog) + style = get_style() + + cat = get_categories(args.multi, args.isol, args.smog) + label = args.label + if args.label == "": + label = get_labels(len(args.fileName)) + + tr = {} + tf = {} + tf = get_files(tf, label, args.fileName) + tr = get_trees(tf, tr, label, args.treeName, True) + + hist_x = {} + hist_y = {} + hist_z = {} + norm = True #to-do + hist_x = get_global(hist_x, tr, "pullx", "#Delta x / #sigma_{x}", + "Candidates Normalized", style, ranges, cat, label, + offset) + hist_y = get_global(hist_y, tr, "pully", "#Delta y / #sigma_{y}", + "Candidates Normalized", style, ranges, cat, label, + offset) + hist_z = get_global(hist_z, tr, "pullz", "#Delta z / #sigma_{z}", + "Candidates Normalized", style, ranges, cat, label, + offset) + + plot_comparison(hist_x, args.prefix, "pullx", cat, label, style, norm, + offset) + plot_comparison(hist_y, args.prefix, "pully", cat, label, style, norm, + offset) + plot_comparison(hist_z, args.prefix, "pullz", cat, label, style, norm, + offset) + + from ROOT import gEnv + gEnv.SetValue("Hist.Binning.1D.x", "100") + + graph = {} + graph["tracks"] = {} + graph["tracks"]["pullx"] = {} + graph["tracks"]["pully"] = {} + graph["tracks"]["pullz"] = {} + graph["tracks"]["pullx"] = get_dependence(graph["tracks"]["pullx"], tr, + "pullx", "nrectrmc", ranges, + style, cat, label, offset) + graph["tracks"]["pully"] = get_dependence(graph["tracks"]["pully"], tr, + "pully", "nrectrmc", ranges, + style, cat, label, offset) + graph["tracks"]["pullz"] = get_dependence(graph["tracks"]["pullz"], tr, + "pullz", "nrectrmc", ranges, + style, cat, label, offset) + + legend = TLegend(0.15, 0.86, 0.88, 0.98) + legend = set_legend(legend, label, graph["tracks"]["pullz"], "sigma") + + lhcbtextpos = (0.9, 0.75) + lhcbtextposdown = (0.9, 0.25) + plot( + graph["tracks"]["pullx"], + "mean", + args.prefix + "_ntracks_mean", + "pullx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["pullx"], + "sigma", + args.prefix + "_ntracks_sigma", + "pullx", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) + + plot( + graph["tracks"]["pully"], + "mean", + args.prefix + "_ntracks_mean", + "pully", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["pully"], + "sigma", + args.prefix + "_ntracks_sigma", + "pully", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) + + plot( + graph["tracks"]["pullz"], + "mean", + args.prefix + "_ntracks_mean", + "pullz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["pullz"], + "sigma", + args.prefix + "_ntracks_sigma", + "pullz", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) + + graph["z"] = {} + graph["z"]["pullx"] = {} + graph["z"]["pully"] = {} + graph["z"]["pullz"] = {} + graph["z"]["pullx"] = get_dependence(graph["z"]["pullx"], tr, "pullx", + "zMC", ranges, style, cat, label, + offset) + graph["z"]["pully"] = get_dependence(graph["z"]["pully"], tr, "pully", + "zMC", ranges, style, cat, label, + offset) + graph["z"]["pullz"] = get_dependence(graph["z"]["pullz"], tr, "pullz", + "zMC", ranges, style, cat, label, + offset) + + plot( + graph["z"]["pullx"], + "mean", + args.prefix + "_z_mean", + "pullx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["pullx"], + "sigma", + args.prefix + "_z_sigma", + "pullx", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) + + plot( + graph["z"]["pully"], + "mean", + args.prefix + "_z_mean", + "pully", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["pully"], + "sigma", + args.prefix + "_z_sigma", + "pully", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) + + plot( + graph["z"]["pullz"], + "mean", + args.prefix + "_z_mean", + "pullz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["pullz"], + "sigma", + args.prefix + "_z_sigma", + "pullz", + cat, + label, + legend, + lhcbtextpos=lhcbtextposdown) diff --git a/Hlt/RecoConf/scripts/PrimaryVertexCheckerResolution.py b/Hlt/RecoConf/scripts/PrimaryVertexCheckerResolution.py new file mode 100644 index 0000000000000000000000000000000000000000..dd5424bf0cdd6b3c8f912f4f2c430f83ee3d2caa --- /dev/null +++ b/Hlt/RecoConf/scripts/PrimaryVertexCheckerResolution.py @@ -0,0 +1,291 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#!/usr/bin/python + +# The script for plotting PV efficinecy as the function +# of various distributions: nTracks, z, r. +# As input the NTuple created by hlt1_reco_pvchecker.py +# is needed. +# +# The efficency is calculated usig TGraphAsymmErrors +# and Bayesian error bars +# +# author: Agnieszka Dziurda (agnieszka.dziurda@cern.ch) +# date: 02/2020 +# +# Example of usage: +# ../../../run python PrimaryVertexCheckerResolution.py +# --file file1.root file2.root --label name1 name2 +# + +import os, sys +import argparse + +from ROOT import gROOT, TLegend + +parser = argparse.ArgumentParser() +parser.add_argument( + '--file', dest='fileName', default="", nargs='+', help='filename to plot') +parser.add_argument( + '--label', dest='label', default="", nargs='+', help='labels for files') +parser.add_argument( + '--tree', + dest='treeName', + default="", + nargs='+', + help='tree name to plot', +) +parser.add_argument( + '--smog', + dest='smog', + default=False, + action='store_true', + help='set true for SMOG') +parser.add_argument( + '--multi', + dest='multi', + default=False, + action='store_true', + help='add multiplicity plots') +parser.add_argument( + '--isol', + dest='isol', + default=False, + action='store_true', + help='add isolated/closed plots') +parser.add_argument( + '--dist', + dest='dist', + default=False, + action='store_true', + help='plot distributions in the canvas') +parser.add_argument( + '--prefix', + dest='prefix', + default="pv_resol", + help='prefix for the plot name', +) +parser.add_argument( + '--dir', + dest='directory', + default=os.getcwd(), + help='tree name to plot', +) + +parser.add_argument( + '--offset', + dest='offset', + default=0, + help='offset for plot colors', +) + + +def get_labels(number_of_files): + label = [] + for i in range(0, number_of_files): + label.append("PVChecker{number}".format(number=str(i + 1))) + return label + + +if __name__ == '__main__': + args = parser.parse_args() + path = args.directory + "/utils/" + sys.path.append(os.path.abspath(path)) + offset = int(args.offset) + + gROOT.SetBatch() + + from pvutils import get_files, get_trees + from pvutils import get_global, plot_comparison + from pvutils import get_dependence + from pvutils import plot + + from pvconfig import get_variable_ranges + from pvconfig import set_legend, get_style, get_categories + + ranges = get_variable_ranges(args.smog) + style = get_style() + + cat = get_categories(args.multi, args.isol, args.smog) + label = args.label + if args.label == "": + label = get_labels(len(args.fileName)) + + tr = {} + tf = {} + tf = get_files(tf, label, args.fileName) + tr = get_trees(tf, tr, label, args.treeName, True) + + hist_x = {} + hist_y = {} + hist_z = {} + norm = True + + hist_x = get_global(hist_x, tr, "dx", "#Delta x [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset) + hist_y = get_global(hist_y, tr, "dy", "#Delta y [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset) + hist_z = get_global(hist_z, tr, "dz", "#Delta z [mm]", + "Candidates Normalized", style, ranges, cat, label, + offset) + + plot_comparison(hist_x, args.prefix, "dx", cat, label, style, norm, offset) + plot_comparison(hist_y, args.prefix, "dy", cat, label, style, norm, offset) + plot_comparison(hist_z, args.prefix, "dz", cat, label, style, norm, offset) + + from ROOT import gEnv + gEnv.SetValue("Hist.Binning.1D.x", "100") + + graph = {} + graph["tracks"] = {} + graph["tracks"]["dx"] = {} + graph["tracks"]["dy"] = {} + graph["tracks"]["dz"] = {} + graph["tracks"]["dx"] = get_dependence(graph["tracks"]["dx"], tr, "dx", + "nrectrmc", ranges, style, cat, + label, offset) + graph["tracks"]["dy"] = get_dependence(graph["tracks"]["dy"], tr, "dy", + "nrectrmc", ranges, style, cat, + label, offset) + graph["tracks"]["dz"] = get_dependence(graph["tracks"]["dz"], tr, "dz", + "nrectrmc", ranges, style, cat, + label, offset) + + legend = TLegend(0.15, 0.86, 0.88, 0.98) + legend = set_legend(legend, label, graph["tracks"]["dz"], "sigma") + + lhcbtextpos = (0.9, 0.75) + plot( + graph["tracks"]["dx"], + "mean", + args.prefix + "_ntracks_mean", + "dx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["dx"], + "sigma", + args.prefix + "_ntracks_sigma", + "dx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + + plot( + graph["tracks"]["dy"], + "mean", + args.prefix + "_ntracks_mean", + "dy", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["dy"], + "sigma", + args.prefix + "_ntracks_sigma", + "dy", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + + plot( + graph["tracks"]["dz"], + "mean", + args.prefix + "_ntracks_mean", + "dz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["tracks"]["dz"], + "sigma", + args.prefix + "_ntracks_sigma", + "dz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + + graph["z"] = {} + graph["z"]["dx"] = {} + graph["z"]["dy"] = {} + graph["z"]["dz"] = {} + graph["z"]["dx"] = get_dependence(graph["z"]["dx"], tr, "dx", "zMC", + ranges, style, cat, label, offset) + graph["z"]["dy"] = get_dependence(graph["z"]["dy"], tr, "dy", "zMC", + ranges, style, cat, label, offset) + graph["z"]["dz"] = get_dependence(graph["z"]["dz"], tr, "dz", "zMC", + ranges, style, cat, label, offset) + + plot( + graph["z"]["dx"], + "mean", + args.prefix + "_z_mean", + "dx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["dx"], + "sigma", + args.prefix + "_z_sigma", + "dx", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + + plot( + graph["z"]["dy"], + "mean", + args.prefix + "_z_mean", + "dy", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["dy"], + "sigma", + args.prefix + "_z_sigma", + "dy", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + + plot( + graph["z"]["dz"], + "mean", + args.prefix + "_z_mean", + "dz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) + plot( + graph["z"]["dz"], + "sigma", + args.prefix + "_z_sigma", + "dz", + cat, + label, + legend, + lhcbtextpos=lhcbtextpos) diff --git a/Hlt/RecoConf/scripts/utils/pvconfig.py b/Hlt/RecoConf/scripts/utils/pvconfig.py new file mode 100644 index 0000000000000000000000000000000000000000..c0cf8b2de682fd253997c2efd1dbad9e00ddca52 --- /dev/null +++ b/Hlt/RecoConf/scripts/utils/pvconfig.py @@ -0,0 +1,220 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### +#!/usr/bin/python + +from ROOT import TH1F +from ROOT import kRed, kBlue, kOrange, kMagenta, kGreen, kCyan + + +def get_categories(multi, isol, smog): + cut = {} + cut["all"] = {"cut": ""} + if isol: + cut["isolated"] = {"cut": "&&isol==1"} + if not smog: + cut["close"] = {"cut": "&&isol==0"} + if multi: + cut["1st"] = {"cut": "&&multimc==1"} + cut["2nd"] = {"cut": "&&multimc==2"} + cut["3rd"] = {"cut": "&&multimc==3"} + cut["4th"] = {"cut": "&&multimc==4"} + cut["5th"] = {"cut": "&&multimc==5"} + + return cut + + +def get_style(): + return { + "color": [kRed, kBlue, kOrange, kMagenta, kGreen, kCyan], + "marker": [21, 20, 22, 23, 24, 25] + } + + +def get_default_tree_name(is_checker): + if is_checker: + return "PrimaryVertexChecker/101" + else: + return "VertexCompare/102" + + +def transfer_variable(dep): + dictionary = { + "pullx": "dx/errx", + "pully": "dy/erry", + "pullz": "dz/errz", + } + return dictionary[dep] + + +def get_variable_ranges(smog): + dictionary = { + "zMC": { + "bins": 50, + "min": -200, + "max": 200 + }, + "rMC": { + "bins": 50, + "min": 0.0, + "max": 0.2 + }, + "xMC": { + "bins": 50, + "min": -0.2, + "max": 0.2 + }, + "yMC": { + "bins": 50, + "min": -0.2, + "max": 0.2 + }, + "dx": { + "bins": 50, + "min": -0.10, + "max": 0.10 + }, + "dy": { + "bins": 50, + "min": -0.10, + "max": 0.10 + }, + "dz": { + "bins": 50, + "min": -0.40, + "max": 0.40 + }, + "pullx": { + "bins": 50, + "min": -3.50, + "max": 3.50 + }, + "pully": { + "bins": 50, + "min": -3.50, + "max": 3.50 + }, + "pullz": { + "bins": 50, + "min": -3.50, + "max": 3.50 + }, + "nrectrmc": { + "bins": 66, + "min": 4.0, + "max": 70.0 + }, + "mtruemcpv": { + "bins": 20, + "min": 0.0, + "max": 20.0 + }, + "nmcpv": { + "bins": 20, + "min": 0.0, + "max": 20.0 + }, + } + if smog: + dictionary["zMC"] = {"bins": 100, "min": -500, "max": 200} + + return dictionary + + +def get_y_axis(dependence): + dictionary = { + "dx": "Resolution #Delta x [#mu m]", + "dy": "Resolution #Delta y [#mu m]", + "dz": "Resolution #Delta z [#mu m]", + "pullx": "Pull #Delta x/#sigma_{x}", + "pully": "Pull #Delta y/#sigma_{y}", + "pullz": "Pull #Delta z/#sigma_{z}" + } + return dictionary[dependence] + + +def get_x_axis(dependence): + dictionary = { + "zMC": "z [mm]", + "rMC": "radial distance [mm]", + "nrectrmc": "number of tracks in Primary Vertex" + } + return dictionary[dependence] + + +def set_style(graph, color, marker, xaxis, yaxis, title): + graph.SetTitle("") + graph.SetLineColor(color) + graph.SetMarkerColor(color) + graph.SetMarkerSize(1.3) + graph.SetMarkerStyle(marker) + graph.GetYaxis().SetTitleOffset(0.85) + if type(graph) == TH1F: + graph.SetFillColor(color) + graph.SetLineWidth(1) + graph.SetStats(False) + graph.GetYaxis().SetTitleOffset(1.1) + graph.GetYaxis().SetTitleSize(0.06) + graph.GetYaxis().SetLabelSize(0.06) + graph.GetXaxis().SetTitleSize(0.06) + graph.GetXaxis().SetLabelSize(0.06) + graph.GetXaxis().SetTitleFont(132) + graph.GetXaxis().SetLabelFont(132) + graph.GetYaxis().SetTitleFont(132) + graph.GetYaxis().SetLabelFont(132) + + if title != "": + graph.SetTitle(title) + if xaxis != "": + graph.GetXaxis().SetTitle(xaxis) + if yaxis != "": + graph.GetYaxis().SetTitle(yaxis) + + +def set_legend_simple(legend, label, hist): + legend.SetTextSize(0.04) + legend.SetTextFont(12) + legend.SetFillColor(4000) + legend.SetShadowColor(0) + legend.SetBorderSize(0) + legend.SetTextFont(132) + legend.SetNColumns(2) + for lab in label: + legend.AddEntry(hist["all"][lab], "{lab}".format(lab=lab), "lep") + + return legend + + +def set_legend(legend, label, gr, gr_type, hist=None, dist=False): + legend.SetTextSize(0.04) + legend.SetTextFont(12) + legend.SetFillColor(4000) + legend.SetShadowColor(0) + legend.SetBorderSize(0) + legend.SetTextFont(132) + legend.SetNColumns(2) + for lab in label: + legend.AddEntry(gr["all"][lab][gr_type], "{lab}".format(lab=lab), + "lep") + if dist: + legend.AddEntry(hist["all"][label[0]]["den"], "Distribution MC", "f") + for lab in label: + legend.AddEntry(hist["all"][lab]["nom"], + "Distribution {lab}".format(lab=lab), "lep") + + return legend + + +def get_text_cor(): + return {"x": [0.17, 0.65, 0.17, 0.65], "y": [0.92, 0.92, 0.75, 0.75]} + + +def basic_cut(): + return "nrectrmc>=4 && dz < 2.0" diff --git a/Hlt/RecoConf/scripts/utils/pvutils.py b/Hlt/RecoConf/scripts/utils/pvutils.py new file mode 100644 index 0000000000000000000000000000000000000000..8112fc72a91ab917d1990dad2f45f377cce373e2 --- /dev/null +++ b/Hlt/RecoConf/scripts/utils/pvutils.py @@ -0,0 +1,469 @@ +############################################################################### +# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration # +# # +# This software is distributed under the terms of the GNU General Public # +#1;95;0c Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". # +# # +# In applying this licence, CERN does not waive the privileges and immunities # +# granted to it by virtue of its status as an Intergovernmental Organization # +# or submit itself to any jurisdiction. # +############################################################################### + +# The utils for processing output of the PrimaryVertexChecker output +# +# author: Agnieszka Dziurda (agnieszka.dziurda@cern.ch) +# date: 02/2020 +# + +from ROOT import TFile, gDirectory, TGraphAsymmErrors, TPaveText, gStyle +from ROOT import TCanvas +from ROOT import gPad, kGray, TLatex + +from pvconfig import get_default_tree_name +from pvconfig import transfer_variable +from pvconfig import get_x_axis, get_y_axis +from pvconfig import set_style +from pvconfig import get_text_cor +from pvconfig import basic_cut + +from array import array + + +def get_files(tf, label, files): + i = 0 + for f in files: + tf[label[i]] = TFile(f) + i += 1 + return tf + + +def get_trees(tf, tr, label, trees, is_checker): + i = 0 + for lab in label: + if len(tf) == len(trees): + tr[lab] = tf[lab].Get(trees[i]) + else: + tr[lab] = tf[lab].Get(get_default_tree_name(is_checker)) + i += 1 + return tr + + +def get_eff(eff, hist, trees, dependence, style, ranges, categories, label, + offset): + + for cat in categories: + eff[cat] = {} + hist[cat] = {} + i = offset + for lab in label: + var_den = '{dependence}>>hist{code}_{dep}_{cat}_{lab}_denom{i}({bins},{mins},{maxs})'.format( + dependence=dependence, + code="den", + dep=dependence, + cat=cat, + lab=lab, + bins=ranges[dependence]["bins"], + mins=ranges[dependence]["min"], + maxs=ranges[dependence]["max"], + i=i) + var_nom = '{dependence}>>hist{code}_{dep}_{cat}_{lab}_num{i}({bins},{mins},{maxs})'.format( + dependence=dependence, + code="nom", + dep=dependence, + cat=cat, + lab=lab, + bins=ranges[dependence]["bins"], + mins=ranges[dependence]["min"], + maxs=ranges[dependence]["max"], + i=i) + + cut_den = 'nrectrmc>=4 && dz < 2.0 {cuts}'.format( + cuts=categories[cat]["cut"]) + cut_nom = cut_den + ' && reco == 1' + + trees[lab].Draw(var_nom, cut_nom) + trees[lab].Draw(var_den, cut_den) + + h_nom = gDirectory.Get( + 'hist{code}_{dep}_{cat}_{lab}_num{i}'.format( + code="nom", dep=dependence, cat=cat, lab=lab, i=i)) + h_den = gDirectory.Get( + 'hist{code}_{dep}_{cat}_{lab}_denom{i}'.format( + code="den", dep=dependence, cat=cat, lab=lab, i=i)) + + g_eff = TGraphAsymmErrors() + g_eff.Divide(h_nom, h_den, "cl=0.683 b(1,1) mode") + + set_style(h_nom, style["color"][i] - 7, style["marker"][i] + 4, + get_x_axis(dependence), "Efficiency", "") + set_style(h_den, kGray, style["marker"][i], get_x_axis(dependence), + "Efficiency", "") + hist[cat][lab] = {} + hist[cat][lab]["nom"] = h_nom + hist[cat][lab]["den"] = h_den + + set_style(g_eff, style["color"][i], style["marker"][i], + get_x_axis(dependence), "Efficiency", "") + eff[cat][lab] = {} + eff[cat][lab]["eff"] = g_eff + i += 1 + + return eff, hist + + +def find_max(gr, gr_type, label): + m = -999999.0 + for lab in label: + if (gr[lab][gr_type].GetYaxis().GetXmax() > m): + m = gr[lab][gr_type].GetYaxis().GetXmax() + + return m + + +def find_min(gr, gr_type, label): + m = 999999.0 + for lab in label: + if (gr[lab][gr_type].GetYaxis().GetXmin() < m): + m = gr[lab][gr_type].GetYaxis().GetXmin() + return m + + +def plot(gr, + gr_type, + prefix, + dependence, + categories, + label, + legend=None, + hist=None, + dist=False, + lhcbtextpos=(0.9, 0.6)): + + for cat in categories: + can = TCanvas( + 'canvas_{depen}_{prefix}_{gr_type}_{cat}'.format( + depen=dependence, prefix=prefix, gr_type=gr_type, cat=cat), + "cR", 1200, 800) + can.SetBottomMargin(0.15) + can.SetLeftMargin(0.12) + can.SetTopMargin(0.15) + if dist: + can.SetTopMargin(0.20) + can.SetRightMargin(0.05) + can.cd() + + maximum = find_max(gr[cat], gr_type, label) + minimum = find_min(gr[cat], gr_type, label) + if (abs(maximum) > abs(minimum) and minimum < 0): + minimum = -maximum + + gr[cat][label[0]][gr_type].GetYaxis().SetRangeUser( + minimum * 1.1, maximum * 1.1) + if (gr_type == "eff"): + gr[cat][label[0]][gr_type].GetYaxis().SetRangeUser(0.0, 1.1) + if (gr_type == "sigma"): + gr[cat][label[0]][gr_type].GetYaxis().SetRangeUser( + 0.0, maximum * 1.1) + gr[cat][label[0]][gr_type].Draw("AP") + + for lab in label: + gr[cat][lab][gr_type].Draw("SAME P") + + if dist: + histmax_den = 1.1 * hist[cat][label[0]]["den"].GetMaximum() + scale = gPad.GetUymax() / histmax_den + hist[cat][label[0]]["den"].Scale(scale * 0.75) + hist[cat][label[0]]["den"].Draw("hist SAME") + + for lab in label: + hist[cat][lab]["nom"].Scale(scale * 0.75) + hist[cat][lab]["nom"].Draw("ep SAME") + + gr[cat][lab][gr_type].Draw("SAME P") + + if legend: + legend.Draw("SAME") + + lhcbtext = TLatex() + lhcbtext.SetTextFont(132) + lhcbtext.SetTextColor(1) + lhcbtext.SetTextSize(0.07) + lhcbtext.SetTextAlign(132) + + lhcbtext.DrawTextNDC(lhcbtextpos[0], lhcbtextpos[1], "LHCb simulation") + + saveName = '{prefix}_{dependence}_{cat}.pdf'.format( + prefix=prefix, dependence=dependence, cat=cat) + can.SaveAs(saveName) + + +def get_global(hist, + trees, + dependence, + x_axis, + y_axis, + style, + ranges, + categories, + label, + offset, + simple=False): + + if dependence.find("pull") > -1: + dep = transfer_variable(dependence) + else: + dep = dependence + + for cat in categories: + hist[cat] = {} + i = offset + for lab in label: + var = '{dependence}>>hist_{dep}_{cat}_{lab}_{i}({bins},{mins},{maxs})'.format( + dependence=dep, + dep=dependence, + cat=cat, + lab=lab, + bins=ranges[dependence]["bins"], + mins=ranges[dependence]["min"], + maxs=ranges[dependence]["max"], + i=i) + + cut = 'nrectrmc>=4 && dz < 2.0 && reco == 1 {cuts}'.format( + cuts=categories[cat]["cut"]) + + trees[lab].Draw(var, cut) + h = gDirectory.Get('hist_{dep}_{cat}_{lab}_{i}'.format( + dep=dependence, cat=cat, lab=lab, i=i)) + if i == 0 and not simple: + col = style["color"][i] - 10 + else: + col = style["color"][i] + set_style(h, col, style["marker"][i], x_axis, y_axis, "") + + hist[cat][lab] = h + i += 1 + + return hist + + +def set_text(text, color, x, y, lab, mean, mean_err, rms, rms_err, scale, + units): + s = 1.0 + if scale: + s = 1000.0 + + ur = "" + um = "#times 10^{-3}" + if units: + ur = "[#mu m]" + um = "[#mu m]" + + text.SetNDC() + text.SetTextFont(132) + text.SetTextColor(color) + text.DrawLatex(x, y * 1.0, lab) + + text.DrawLatex( + x, y * 0.95, "#mu = ({0:0.2f} #pm {1:0.2f}) {unit}".format( + mean * 1000.0, mean_err * 1000.0, unit=um)) + text.DrawLatex( + x, y * 0.90, "#sigma = ({0:0.2f} #pm {1:0.2f}) {unit}".format( + rms * s, rms_err * s, unit=ur)) + return text + + +def plot_comparison(hist, + prefix, + dependence, + categories, + label, + style, + norm, + offset, + simple=False, + legend=None): + + for cat in categories: + can = TCanvas('canvas_{depen}_{cat}'.format(depen=dependence, cat=cat), + "cR", 1200, 800) + can.SetBottomMargin(0.15) + can.SetLeftMargin(0.15) + can.SetTopMargin(0.20) + can.SetRightMargin(0.05) + + can.cd() + cor = get_text_cor() + + hist[cat][label[0]].GetYaxis().SetRangeUser( + 0.0, hist[cat][label[0]].GetMaximum() * 1.1) + scale = True + units = True + if "pull" in dependence: + scale = False + units = False + if norm: + if not simple: + hist[cat][label[0]].DrawNormalized("hist") + else: + hist[cat][label[0]].DrawNormalized("PE") + i = offset + for lab in label: + hist[cat][lab].DrawNormalized("SAME PE") + if not simple: + text = TLatex() + text = set_text(text, style["color"][i], cor["x"][i], + cor["y"][i], lab, hist[cat][lab].GetMean(), + hist[cat][lab].GetMeanError(), + hist[cat][lab].GetRMS(), + hist[cat][lab].GetRMSError(), scale, units) + i += 1 + + else: + hist[cat][label[0]].Draw("hist") + for lab in label: + hist[cat][lab].Draw("SAME PE") + if legend: + legend.Draw("SAME") + + pavetext = TPaveText(0.72, 0.77 - gStyle.GetPadTopMargin(), 0.87, + 0.89 - gStyle.GetPadTopMargin(), "NBNDC") + pavetext.AddText("LHCb simulation") + pavetext.SetFillColor(0) + pavetext.SetTextSize(0.06) + pavetext.SetTextFont(132) + pavetext.Draw() + + can.Update() + saveName = '{prefix}_{dependence}_{cat}.pdf'.format( + prefix=prefix, dependence=dependence, cat=cat) + can.SaveAs(saveName) + + +def get_robust_sigma(hist, resol): + + y = array('d', [0.] * 3) + x = array('d', [0.] * 3) + x[0] = 0.25 + x[1] = 0.50 + x[2] = 0.75 + + hist.GetQuantiles(3, y, x) + + _median = y[1] + _approxstdev = (y[2] - y[0]) / 1.34898 + #factor gives correspondence between IQR and stdev for a Gaussian + mult = 4.0 + + histclone = hist.Clone() + nb = histclone.GetNbinsX() + + for i in range(1, nb + 1): + if ((histclone.GetBinCenter(i) < (_median - _approxstdev * mult)) + or (histclone.GetBinCenter(i) > + (_median + _approxstdev * mult))): + histclone.SetBinContent(i, 0) + + if resol: + robustsigma = { + "sigma": { + "var": histclone.GetRMS() * 1000.0, + "err": histclone.GetRMSError() * 1000.0 + }, + "mean": { + "var": histclone.GetMean() * 1000.0, + "err": histclone.GetMeanError() * 1000.0 + } + } + + else: + robustsigma = { + "sigma": { + "var": histclone.GetRMS() * 1.0, + "err": histclone.GetRMSError() * 1.0 + }, + "mean": { + "var": histclone.GetMean(), + "err": histclone.GetMeanError() + } + } + + return robustsigma + + +def get_dependence(graph, trees, variable, dependence, ranges, style, + categories, label, offset): + + if variable.find("pull") > -1: + var = transfer_variable(dependence) + else: + var = variable + bin_width = (ranges[dependence]["max"] - + ranges[dependence]["min"]) / ranges[dependence]["bins"] + + resol = True + mean = "Mean [#mu m]" + if "pull" in variable: + resol = False + mean = "Mean" + + for cat in categories: + graph[cat] = {} + i = offset + for lab in label: + graph[cat][lab] = {} + graph_rms = TGraphAsymmErrors(ranges[dependence]["bins"]) + graph_mean = TGraphAsymmErrors(ranges[dependence]["bins"]) + + bin_min = ranges[dependence]["min"] + for b in range(0, ranges[dependence]["bins"]): + + cut_bin = '{variable1}>={mins} && {variable2}<{maxs}'.format( + variable1=dependence, + mins=bin_min, + variable2=dependence, + maxs=bin_min + bin_width) + + hist_string = '{dependence}>>hist_{dep}_{cat}_{i}_{number}'.format( + dependence=var, + dep=variable, + cat=cat, + i=i, + number=int(b), + ) + + cut = '{basic} && reco == 1 && {bin_cut} {cuts}'.format( + basic=basic_cut(), + bin_cut=cut_bin, + cuts=categories[cat]["cut"]) + + trees[lab].Draw(hist_string, cut) + + bin_min += bin_width + + h = gDirectory.Get('hist_{dep}_{cat}_{i}_{number}'.format( + dep=variable, cat=cat, i=i, number=int(b))) + + robust_sigma = get_robust_sigma(h, resol) + + graph_mean.SetPoint(b, bin_min + bin_width / 2.0, + robust_sigma["mean"]["var"]) + graph_mean.SetPointError(b, bin_width / 2.0, bin_width / 2.0, + robust_sigma["mean"]["err"], + robust_sigma["mean"]["err"]) + + graph_rms.SetPoint(b, bin_min + bin_width / 2.0, + robust_sigma["sigma"]["var"]) + graph_rms.SetPointError(b, bin_width / 2.0, bin_width / 2.0, + robust_sigma["sigma"]["err"], + robust_sigma["sigma"]["err"]) + + set_style(graph_rms, style["color"][i], style["marker"][i], + get_x_axis(dependence), get_y_axis(variable), "") + set_style(graph_mean, style["color"][i], style["marker"][i], + get_x_axis(dependence), mean, "") + + graph[cat][lab]["sigma"] = graph_rms + graph[cat][lab]["mean"] = graph_mean + i += 1 + + return graph