From cbd6e7f5a5dc49fd490e124a9e7face84e66d782 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 19 Mar 2019 12:15:25 +0100 Subject: [PATCH 01/13] added option to set label size in plot --- CombineTools/scripts/plotImpacts.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CombineTools/scripts/plotImpacts.py b/CombineTools/scripts/plotImpacts.py index 9479ea43205..4e25c0dea5c 100755 --- a/CombineTools/scripts/plotImpacts.py +++ b/CombineTools/scripts/plotImpacts.py @@ -23,6 +23,7 @@ parser.add_argument('--color-groups', default=None, help='Comma separated list of GROUP=COLOR') parser.add_argument("--pullDef", default=None, help="Choose the definition of the pull, see HiggsAnalysis/CombinedLimit/python/calculate_pulls.py for options") parser.add_argument('--POI', default=None, help='Specify a POI to draw') +parser.add_argument("--set-label-size", dest = "labelsize",default = 0.022, help = "Set size of labels on the left hand side (default = 0.022", type = float) args = parser.parse_args() if args.transparent: @@ -219,7 +220,7 @@ def GetRounded(nom, e_hi, e_lo): else: plot.Set(h_pulls.GetXaxis(), TitleSize=0.04, LabelSize=0.03, Title='(#hat{#theta}-#theta_{0})/#Delta#theta') - plot.Set(h_pulls.GetYaxis(), LabelSize=0.021, TickLength=0.0) + plot.Set(h_pulls.GetYaxis(), LabelSize=args.labelsize, TickLength=0.0) h_pulls.GetYaxis().LabelsOption('v') h_pulls.Draw() @@ -324,4 +325,4 @@ def GetRounded(nom, e_hi, e_lo): extra = '(' if page == n - 1: extra = ')' - canv.Print('.pdf%s' % extra) + canv.Print('.pdf%s' % extra) \ No newline at end of file From 5f57af8a01cf68571463b919472d6c13cf987bc4 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 18 Jan 2022 10:23:09 +0100 Subject: [PATCH 02/13] implemented option to merge processes at runtime --- .../bin/PostFitShapesFromWorkspace.cpp | 58 ++++++++++++++++++- 1 file changed, 57 insertions(+), 1 deletion(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index a7b626d1a25..9e42236b3ac 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -47,6 +47,9 @@ int main(int argc, char* argv[]) { bool skip_proc_errs = false; bool total_shapes = false; std::vector reverse_bins_; + // Containers to parse processes that are to be merged at runtime + std::vector input_merge_procs_; + std::map merged_procs; po::options_description help_config("Help"); help_config.add_options() @@ -103,7 +106,9 @@ int main(int argc, char* argv[]) { ("total-shapes", po::value(&total_shapes)->default_value(total_shapes)->implicit_value(true), "Save signal- and background shapes added for all channels/categories") - ("reverse-bins", po::value>(&reverse_bins_)->multitoken(), "List of bins to reverse the order for"); + ("reverse-bins", po::value>(&reverse_bins_)->multitoken(), "List of bins to reverse the order for") + ("merge-procs,p", po::value>(&input_merge_procs_)->multitoken(), + "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'"); po::variables_map vm; @@ -140,6 +145,8 @@ int main(int argc, char* argv[]) { // Create CH instance and parse the workspace ch::CombineHarvester cmb; cmb.SetFlag("workspaces-use-clone", true); + // Allow regex expressions to combine processes on the fly + cmb.SetFlag("filters-use-regex", true); ch::ParseCombineWorkspace(cmb, *ws, "ModelConfig", data, false); // Only evaluate in case parameters to freeze are provided @@ -170,6 +177,13 @@ int main(int argc, char* argv[]) { } // cmb.GetParameter("r")->set_frozen(true); + // parse processes that are to be merged + for (auto& in: input_merge_procs_){ + vector parts; + boost::split(parts, in, boost::is_any_of("=")); + merged_procs[parts[0]] = parts[1]; + } + ch::CombineHarvester cmb_card; cmb_card.SetFlag("workspaces-use-clone",true); if (datacard != "") { @@ -248,6 +262,28 @@ int main(int argc, char* argv[]) { cmb_bin.cp().process({proc}).GetShapeWithUncertainty(); } } + // Create prefit shapes for merged processes + for (auto iter: merged_procs){ + // First element of the iterator is the name of the merged process + auto proc=iter.first; + std::cout << ">> Doing prefit: " << bin << "," << proc << std::endl; + // Second element is the regex expression for the processes + // that are to be merged + auto proc_regex = iter.second; + auto cmb_proc = cmb_bin.cp().process({proc_regex}); + // First check for matches + if (cmb_proc.process_set().size() == 0){ + std::cout << ">> WARNING: found no processes matching " << proc << std::endl; + continue; + } + if (skip_proc_errs) { + pre_shapes[bin][proc] = + cmb_proc.GetShape(); + } else { + pre_shapes[bin][proc] = + cmb_proc.GetShapeWithUncertainty(); + } + } // The fill total signal and total bkg hists std::cout << ">> Doing prefit: " << bin << "," << "TotalBkg" << std::endl; @@ -362,6 +398,26 @@ int main(int argc, char* argv[]) { : cmb_proc.GetShapeWithUncertainty(res, samples); } } + + // Generate postfit distributions for merged processes + for (auto iter: merged_procs){ + auto proc=iter.first; + std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; + auto proc_regex = iter.second; + auto cmb_proc = cmb_bin.cp().process({proc_regex}); + if (cmb_proc.process_set().size() == 0){ + std::cout << ">> WARNING: found no processes matching " << proc << std::endl; + continue; + } + if (skip_proc_errs) { + post_shapes[bin][proc] = cmb_proc.GetShape(); + } else { + post_shapes[bin][proc] = + sampling ? cmb_proc.GetShapeWithUncertainty(res, samples) + : cmb_proc.GetShapeWithUncertainty(); + } + } + if (!no_sampling && covariance) { post_yield_cov[bin] = cmb_bin.GetRateCovariance(res, samples); post_yield_cor[bin] = cmb_bin.GetRateCorrelation(res, samples); From 5ce7aa8294af1b841d0c0d87b9332df19556c750 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 18 Jan 2022 10:52:02 +0100 Subject: [PATCH 03/13] implemented verbose mode in generation of uncertainty band for debugging --- CombineTools/interface/CombineHarvester.h | 4 +- CombineTools/src/CombineHarvester_Evaluate.cc | 51 ++++++++++++++++--- CombineTools/src/CombineHarvester_Python.cc | 2 +- 3 files changed, 46 insertions(+), 11 deletions(-) diff --git a/CombineTools/interface/CombineHarvester.h b/CombineTools/interface/CombineHarvester.h index 1f6403bc687..5ba6f29e370 100644 --- a/CombineTools/interface/CombineHarvester.h +++ b/CombineTools/interface/CombineHarvester.h @@ -362,8 +362,8 @@ class CombineHarvester { * version of this method instead - this method will be removed in an * upcoming release */ - TH1F GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples); - TH1F GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples); + TH1F GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples, const bool& verbose = false); + TH1F GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples, const bool& verbose = false); TH1F GetObservedShape(); TH2F GetRateCovariance(RooFitResult const& fit, unsigned n_samples); diff --git a/CombineTools/src/CombineHarvester_Evaluate.cc b/CombineTools/src/CombineHarvester_Evaluate.cc index 6a6dea751a4..f260a58c4df 100644 --- a/CombineTools/src/CombineHarvester_Evaluate.cc +++ b/CombineTools/src/CombineHarvester_Evaluate.cc @@ -14,6 +14,7 @@ #include "TDirectory.h" #include "TH1.h" #include "TH2.h" +#include "TString.h" #include "CombineHarvester/CombineTools/interface/Observation.h" #include "CombineHarvester/CombineTools/interface/Process.h" #include "CombineHarvester/CombineTools/interface/Systematic.h" @@ -126,16 +127,21 @@ TH1F CombineHarvester::GetShapeWithUncertainty() { } TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const* fit, - unsigned n_samples) { - return GetShapeWithUncertainty(*fit, n_samples); + unsigned n_samples, + const bool& verbose) { + return GetShapeWithUncertainty(*fit, n_samples, verbose); } TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, - unsigned n_samples) { + unsigned n_samples, + const bool& verbose) +{ auto lookup = GenerateProcSystMap(); TH1F shape = GetShapeInternal(lookup); + std::vector full_rand_shape_summary; for (int i = 1; i <= shape.GetNbinsX(); ++i) { shape.SetBinError(i, 0.0); + full_rand_shape_summary.push_back(""); } // Create a backup copy of the current parameter values auto backup = GetParameters(); @@ -152,22 +158,51 @@ TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, for (unsigned n = 0; n < p_vec.size(); ++n) { r_vec[n] = dynamic_cast(rands.at(n)); p_vec[n] = GetParameter(r_vec[n]->GetName()); + // if(! p_vec[n]) std::cout << "could not initialize parameter '" << r_vec[n]->GetName() << "'" << std::endl; } + TString tmp; + // Main loop through n_samples for (unsigned i = 0; i < n_samples; ++i) { // Randomise and update values fit.randomizePars(); + std::string rand_shape_summary = ""; for (int n = 0; n < n_pars; ++n) { - if (p_vec[n]) p_vec[n]->set_val(r_vec[n]->getVal()); + if (p_vec[n] && !p_vec[n]->frozen()) { + p_vec[n]->set_val(r_vec[n]->getVal()); + tmp.Form("setting parameter '%s' to %f\n", p_vec[n]->name().c_str(), r_vec[n]->getVal()); + // std::cout << tmp.Data(); + rand_shape_summary = tmp.Data(); + // std::cout << rand_shape_summary.c_str() << std::endl; + TH1F rand_shape = this->GetShapeInternal(lookup); + for (int j = 1; j <= rand_shape.GetNbinsX(); ++j) { + tmp.Form("\tBin Content %i: %f\n", j, rand_shape.GetBinContent(j)); + // std::cout << tmp.Data(); + full_rand_shape_summary.at(j-1) = rand_shape_summary; + full_rand_shape_summary.at(j-1) += tmp.Data(); + } + } + // else if (!p_vec[n]) { + // std::cout << "could not find parameter '" << r_vec[n]->GetName() << "'"<< std::endl; + // } } - TH1F rand_shape = this->GetShapeInternal(lookup); - for (int i = 1; i <= shape.GetNbinsX(); ++i) { + for (int j = 1; j <= shape.GetNbinsX(); ++j) { + double err = - std::fabs(rand_shape.GetBinContent(i) - shape.GetBinContent(i)); - shape.SetBinError(i, err*err + shape.GetBinError(i)); + std::fabs(rand_shape.GetBinContent(j) - shape.GetBinContent(j)); + if (std::fabs(err/shape.GetBinContent(j)) >= 0.5 and verbose){ + std::cout << "rel. error > 0.5 detected in bin " << j << " for toy " << i << std::endl; + std::cout << "nominal bin content (" << j << "):\t" << shape.GetBinContent(j) <UpdateParameters(backup); } for (int i = 1; i <= shape.GetNbinsX(); ++i) { shape.SetBinError(i, std::sqrt(shape.GetBinError(i)/double(n_samples))); diff --git a/CombineTools/src/CombineHarvester_Python.cc b/CombineTools/src/CombineHarvester_Python.cc index 93002b02a8d..70d9c48d7c1 100644 --- a/CombineTools/src/CombineHarvester_Python.cc +++ b/CombineTools/src/CombineHarvester_Python.cc @@ -158,7 +158,7 @@ TH1F (CombineHarvester::*Overload1_GetShapeWithUncertainty)( void) = &CombineHarvester::GetShapeWithUncertainty; TH1F (CombineHarvester::*Overload2_GetShapeWithUncertainty)( - RooFitResult const&, unsigned) = &CombineHarvester::GetShapeWithUncertainty; + RooFitResult const&, unsigned, bool const&) = &CombineHarvester::GetShapeWithUncertainty; void (CombineHarvester::*Overload1_UpdateParameters)( RooFitResult const&) = &CombineHarvester::UpdateParameters; From 31c1e5cd127890a2b4802e2c214c025cd77384cf Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 18 Jan 2022 10:59:53 +0100 Subject: [PATCH 04/13] implemented option for verbose generation of postfit uncertainty bands --- .../bin/PostFitShapesFromWorkspace.cpp | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index 9e42236b3ac..16b7a00e17b 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -46,6 +46,7 @@ int main(int argc, char* argv[]) { bool skip_prefit = false; bool skip_proc_errs = false; bool total_shapes = false; + bool verbose = false; std::vector reverse_bins_; // Containers to parse processes that are to be merged at runtime std::vector input_merge_procs_; @@ -106,6 +107,9 @@ int main(int argc, char* argv[]) { ("total-shapes", po::value(&total_shapes)->default_value(total_shapes)->implicit_value(true), "Save signal- and background shapes added for all channels/categories") + ("verbose,v", + po::value(&verbose)->default_value(verbose)->implicit_value(true), + "Genererate additional output if uncertainties in a given bin are large") ("reverse-bins", po::value>(&reverse_bins_)->multitoken(), "List of bins to reverse the order for") ("merge-procs,p", po::value>(&input_merge_procs_)->multitoken(), "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'"); @@ -356,15 +360,15 @@ int main(int argc, char* argv[]) { std::cout << ">> Doing postfit: TotalBkg" << std::endl; post_shapes_tot["TotalBkg"] = no_sampling ? cmb_bkgs.GetShapeWithUncertainty() - : cmb_bkgs.GetShapeWithUncertainty(res, samples); + : cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose); std::cout << ">> Doing postfit: TotalSig" << std::endl; post_shapes_tot["TotalSig"] = no_sampling ? cmb_sigs.GetShapeWithUncertainty() - : cmb_sigs.GetShapeWithUncertainty(res, samples); + : cmb_sigs.GetShapeWithUncertainty(res, samples, verbose); std::cout << ">> Doing postfit: TotalProcs" << std::endl; post_shapes_tot["TotalProcs"] = no_sampling ? cmb.cp().GetShapeWithUncertainty() - : cmb.cp().GetShapeWithUncertainty(res, samples); + : cmb.cp().GetShapeWithUncertainty(res, samples, verbose); if (datacard != "") { TH1F ref = cmb_card.cp().GetObservedShape(); @@ -395,7 +399,7 @@ int main(int argc, char* argv[]) { } else { post_shapes[bin][proc] = no_sampling ? cmb_proc.GetShapeWithUncertainty() - : cmb_proc.GetShapeWithUncertainty(res, samples); + : cmb_proc.GetShapeWithUncertainty(res, samples, verbose); } } @@ -413,7 +417,7 @@ int main(int argc, char* argv[]) { post_shapes[bin][proc] = cmb_proc.GetShape(); } else { post_shapes[bin][proc] = - sampling ? cmb_proc.GetShapeWithUncertainty(res, samples) + sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) : cmb_proc.GetShapeWithUncertainty(); } } @@ -428,15 +432,15 @@ int main(int argc, char* argv[]) { std::cout << ">> Doing postfit: " << bin << "," << "TotalBkg" << std::endl; post_shapes[bin]["TotalBkg"] = no_sampling ? cmb_bkgs.GetShapeWithUncertainty() - : cmb_bkgs.GetShapeWithUncertainty(res, samples); + : cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose); std::cout << ">> Doing postfit: " << bin << "," << "TotalSig" << std::endl; post_shapes[bin]["TotalSig"] = no_sampling ? cmb_sigs.GetShapeWithUncertainty() - : cmb_sigs.GetShapeWithUncertainty(res, samples); + : cmb_sigs.GetShapeWithUncertainty(res, samples, verbose); std::cout << ">> Doing postfit: " << bin << "," << "TotalProcs" << std::endl; post_shapes[bin]["TotalProcs"] = no_sampling ? cmb_bin.cp().GetShapeWithUncertainty() - : cmb_bin.cp().GetShapeWithUncertainty(res, samples); + : cmb_bin.cp().GetShapeWithUncertainty(res, samples, verbose); if (datacard != "") { TH1F ref = cmb_card.cp().bin({bin}).GetObservedShape(); From a88c57745662ac489b1d88a801e5e48f8de8df39 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 7 Mar 2023 17:23:48 +0100 Subject: [PATCH 05/13] restructure code, implement useful additional options --- .../bin/PostFitShapesFromWorkspace.cpp | 746 ++++++++++++------ 1 file changed, 486 insertions(+), 260 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index 16b7a00e17b..e022a66005e 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -1,8 +1,14 @@ #include +#include + #include "boost/program_options.hpp" #include "boost/format.hpp" #include "TSystem.h" #include "TH2F.h" +#include "RooRealVar.h" +#include "TMath.h" +#include "TString.h" +#include "RooMsgService.h" #include "CombineHarvester/CombineTools/interface/CombineHarvester.h" #include "CombineHarvester/CombineTools/interface/ParseCombineWorkspace.h" #include "CombineHarvester/CombineTools/interface/TFileIO.h" @@ -26,6 +32,129 @@ void ReverseBins(TH1F & h) { // return h; } +void load_shapes_and_yields( + ch::CombineHarvester& cmb, + map& shape_container, + vector& yield_container, + const std::string& proc, + const bool& skip_proc_errs=true, + const bool& sampling=false, + const int& samples=0, + RooFitResult* res=NULL, + const bool& verbose=false){ + TString helper; + if (skip_proc_errs) { + shape_container[proc] = + cmb.GetShape(); + } else { + shape_container[proc] = + sampling ? cmb.GetShapeWithUncertainty(*res, samples, verbose) + : cmb.GetShapeWithUncertainty(); + } + helper.Form("%s_%s", "yield" , proc.c_str()); + yield_container.push_back(RooRealVar(helper, helper, + cmb.GetRate())); + yield_container.back().setError(sampling ? cmb.GetUncertainty(*res, samples) + : cmb.GetUncertainty()); +} + +void do_complete_set( + ch::CombineHarvester& cmb_bin, + TFile& outfile, + const std::string& output_prefix, + map& shape_container, + vector& yield_container, + const std::map& merged_procs, + const std::set skip_procs, + ch::CombineHarvester& cmb_card, + const std::string& datacard="", + const bool& skip_proc_errs=true, + const bool& factors=false, + const bool& reverse_bins=false, + const bool& sampling=false, + const int& samples=0, + RooFitResult* res=NULL, + const bool& verbose=false + ){ + TString helper; + // ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); + // This next line is a temporary fix for models with parameteric RooFit pdfs + // - we try and set the number of bins to evaluate the pdf to be the same as + // the number of bins in data + // cmb_bin.SetPdfBins(cmb_bin.GetObservedShape().GetNbinsX()); + + // Fill the data and process histograms + shape_container["data_obs"] = cmb_bin.GetObservedShape(); + yield_container.push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb_bin.GetObservedRate())); + yield_container.back().setError(TMath::Power(cmb_bin.GetObservedRate(), 0.5)); + for (auto proc : cmb_bin.process_set()) { + if(skip_procs.find(proc) == skip_procs.end()){ + std::cout << ">> Doing " << output_prefix.c_str() << "," << proc << std::endl; + load_shapes_and_yields(cmb_bin.cp().process({proc}), shape_container, yield_container, proc, skip_proc_errs, sampling, samples, res, verbose); + } + else{ + std::cout << ">> Skipping processes " << proc << std::endl; + } + } + for (auto iter: merged_procs){ + auto proc=iter.first; + std::cout << ">> Doing " << output_prefix.c_str() << "," << proc << std::endl; + auto proc_regex = iter.second; + auto cmb_proc = cmb_bin.cp().process({proc_regex}); + if (cmb_proc.process_set().size() == 0){ + std::cout << ">> WARNING: found no processes matching " << proc << std::endl; + continue; + } + load_shapes_and_yields(cmb_proc, shape_container, yield_container, proc, skip_proc_errs, sampling, samples, res, verbose); + } + + // The fill total signal and total bkg hists + std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalBkg" << std::endl; + load_shapes_and_yields(cmb_bin.cp().backgrounds(), shape_container, yield_container, "TotalBkg", false, sampling, samples, res, verbose); + + // Print out the relative uncert. on the bkg + if (factors) { + cout << boost::format("%-25s %-32s\n") % "Bin" % + "Total relative bkg uncert."; + cout << string(58, '-') << "\n"; + + double rate = yield_container.back().getVal(); + double err = yield_container.back().getError(); + cout << boost::format("%-25s %-10.5f") % output_prefix.c_str() % + (rate > 0. ? (err / rate) : 0.) << std::endl; + } + + std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalSig" << std::endl; + load_shapes_and_yields(cmb_bin.cp().signals(), shape_container, yield_container, "TotalSig", false, sampling, samples, res, verbose); + + std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalProcs" << std::endl; + load_shapes_and_yields(cmb_bin, shape_container, yield_container, "TotalProcs", false, sampling, samples, res, verbose); + + + if (datacard != "") { + TH1F ref = cmb_card.GetObservedShape(); + for (auto & it : shape_container) { + it.second = ch::RestoreBinning(it.second, ref); + } + } + + if(reverse_bins){ + for (auto it = shape_container.begin(); it != shape_container.end(); ++it) { + ReverseBins(it->second); + } + } + // Can write these straight into the output file + outfile.cd(); + for (auto& iter : shape_container) { + ch::WriteToTFile(&(iter.second), &outfile, output_prefix + "/" + iter.first); + } + for (auto& yield: yield_container){ + helper.Form("%s/%s", output_prefix.c_str() , yield.GetName()); + ch::WriteToTFile(&(yield), &outfile, helper.Data()); + } +} + + int main(int argc, char* argv[]) { // Need this to read combine workspaces gSystem->Load("libHiggsAnalysisCombinedLimit"); @@ -36,7 +165,6 @@ int main(int argc, char* argv[]) { string mass = ""; bool postfit = false; bool sampling = false; - bool no_sampling = false; string output = ""; bool factors = false; unsigned samples = 500; @@ -48,10 +176,13 @@ int main(int argc, char* argv[]) { bool total_shapes = false; bool verbose = false; std::vector reverse_bins_; - // Containers to parse processes that are to be merged at runtime + std::vector bins_; + std::set skip_procs; + std::vector skip_procs_; std::vector input_merge_procs_; std::map merged_procs; + po::options_description help_config("Help"); help_config.add_options() ("help,h", "produce help message"); @@ -82,10 +213,7 @@ int main(int argc, char* argv[]) { "Create post-fit histograms in addition to pre-fit") ("sampling", po::value(&sampling)->default_value(sampling)->implicit_value(true), - "Use the cov. matrix sampling method for the post-fit uncertainty (deprecated, this is the default)") - ("no-sampling", - po::value(&no_sampling)->default_value(no_sampling)->implicit_value(true), - "Do not use the cov. matrix sampling method for the post-fit uncertainty") + "Use the cov. matrix sampling method for the post-fit uncertainty") ("samples", po::value(&samples)->default_value(samples), "Number of samples to make in each evaluate call") @@ -104,16 +232,23 @@ int main(int argc, char* argv[]) { ("skip-proc-errs", po::value(&skip_proc_errs)->default_value(skip_proc_errs)->implicit_value(true), "Skip evaluation of errors on individual processes") - ("total-shapes", - po::value(&total_shapes)->default_value(total_shapes)->implicit_value(true), - "Save signal- and background shapes added for all channels/categories") ("verbose,v", po::value(&verbose)->default_value(verbose)->implicit_value(true), "Genererate additional output if uncertainties in a given bin are large") + ("total-shapes", + po::value(&total_shapes)->default_value(total_shapes)->implicit_value(true), + "Save signal- and background shapes added for all channels/categories") ("reverse-bins", po::value>(&reverse_bins_)->multitoken(), "List of bins to reverse the order for") + ("bins", po::value>(&bins_)->multitoken(), "List of bins to produce shapes for (default: all bins") ("merge-procs,p", po::value>(&input_merge_procs_)->multitoken(), - "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'"); + "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'") + ("skip-procs,s", po::value>(&skip_procs_)->multitoken(), + "Skip these processes. Regex expression allowed. Format: 'expression'. Can be called multiple times"); + if (sampling && !postfit) { + throw logic_error( + "Can't sample the fit covariance matrix for pre-fit!"); + } po::variables_map vm; @@ -125,7 +260,7 @@ int main(int argc, char* argv[]) { if (vm.count("help")) { cout << config << "\nExample usage:\n"; cout << "PostFitShapesFromWorkspace.root -d htt_mt_125.txt -w htt_mt_125.root -o htt_mt_125_shapes.root -m 125 " - "-f mlfit.root:fit_s --postfit --print\n"; + "-f mlfit.root:fit_s --postfit --sampling --print\n"; return 1; } @@ -133,10 +268,6 @@ int main(int argc, char* argv[]) { po::store(po::command_line_parser(argc, argv).options(config).run(), vm); po::notify(vm); - if (sampling) { - std::cout<<"WARNING: the default behaviour of PostFitShapesFromWorkspace is to use the covariance matrix sampling method for the post-fit uncertainty. The option --sampling is deprecated and will be removed in future versions of CombineHarvester"<(gDirectory->Get("w")); @@ -149,33 +280,61 @@ int main(int argc, char* argv[]) { // Create CH instance and parse the workspace ch::CombineHarvester cmb; cmb.SetFlag("workspaces-use-clone", true); - // Allow regex expressions to combine processes on the fly cmb.SetFlag("filters-use-regex", true); ch::ParseCombineWorkspace(cmb, *ws, "ModelConfig", data, false); - + RooMsgService::instance().setSilentMode(true); // Only evaluate in case parameters to freeze are provided if(! freeze_arg.empty()) { vector freeze_vec; boost::split(freeze_vec, freeze_arg, boost::is_any_of(",")); + vector parameters; + for (auto& par : cmb.GetParameters()) + { + parameters.push_back(par.name()); + } for (auto const& item : freeze_vec) { vector parts; boost::split(parts, item, boost::is_any_of("=")); - if (parts.size() == 1) { - ch::Parameter *par = cmb.GetParameter(parts[0]); - if (par) par->set_frozen(true); - else throw std::runtime_error( + auto current_expr = parts[0]; + + // check for regex syntax: rgx{regex} + if (boost::starts_with(current_expr, "rgx{") && boost::ends_with(current_expr, "}")) { + bool matched = false; + + std::string reg_esp = current_expr.substr(4, current_expr.size()-5); + std::cout<<"interpreting "<set_val(boost::lexical_cast(parts[1])); + } + par->set_frozen(true); + } + } + // if not match is found, throw runtime error + if(!matched){ + throw std::runtime_error( FNERROR("Requested variable to freeze does not exist in workspace")); + } } else { - if (parts.size() == 2) { - ch::Parameter *par = cmb.GetParameter(parts[0]); - if (par) { + ch::Parameter *par = cmb.GetParameter(current_expr); + if (par) { + std::cout << "freezing parameter '" << par->name() << "'" << std::endl; + if (parts.size() == 2) { par->set_val(boost::lexical_cast(parts[1])); - par->set_frozen(true); } - else throw std::runtime_error( - FNERROR("Requested variable to freeze does not exist in workspace")); + par->set_frozen(true); } + else throw std::runtime_error( + FNERROR("Requested variable to freeze does not exist in workspace")); } } } @@ -203,137 +362,94 @@ int main(int argc, char* argv[]) { } return no_shape; }); + vector bins; + if (bins_.size() == 0) + { + auto bin_set = cmb.cp().bin_set(); + std::copy(bin_set.begin(), bin_set.end(), std::back_inserter(bins)); + } + else{ + bins = bins_; + } - auto bins = cmb.cp().bin_set(); + // parse processes + std::set process_set; + for(auto& expr: skip_procs_){ + process_set = cmb.cp().process({expr}).process_set(); + skip_procs.insert(process_set.begin(), process_set.end()); + } TFile outfile(output.c_str(), "RECREATE"); TH1::AddDirectory(false); - // Create a map of maps for storing histograms in the form: - // pre_shapes[][] - map> pre_shapes; - - // Also create a simple map for storing total histograms, summed - // over all bins, in the form: - // pre_shapes_tot[] - map pre_shapes_tot; - // We can always do the prefit version, // Loop through the bins writing the shapes to the output file + map> pre_shapes; + map> pre_yields; + TString helper; if (!skip_prefit) { if(total_shapes){ - pre_shapes_tot["data_obs"] = cmb.GetObservedShape(); - // Then fill total signal and total bkg hists - std::cout << ">> Doing prefit: TotalBkg" << std::endl; - pre_shapes_tot["TotalBkg"] = - cmb.cp().backgrounds().GetShapeWithUncertainty(); - std::cout << ">> Doing prefit: TotalSig" << std::endl; - pre_shapes_tot["TotalSig"] = - cmb.cp().signals().GetShapeWithUncertainty(); - std::cout << ">> Doing prefit: TotalProcs" << std::endl; - pre_shapes_tot["TotalProcs"] = - cmb.cp().GetShapeWithUncertainty(); - - if (datacard != "") { - TH1F ref = cmb_card.cp().GetObservedShape(); - for (auto & it : pre_shapes_tot) { - it.second = ch::RestoreBinning(it.second, ref); - } - } - - // Can write these straight into the output file - outfile.cd(); - for (auto& iter : pre_shapes_tot) { - ch::WriteToTFile(&(iter.second), &outfile, "prefit/" + iter.first); - } + pre_shapes["total"] = map(); + pre_yields["total"] = vector(); + do_complete_set(cmb, outfile, "total_prefit", pre_shapes["total"], pre_yields["total"], merged_procs, skip_procs, cmb_card, datacard, skip_proc_errs, factors); + // pre_yields["total"] = vector(); + // pre_shapes_tot["data_obs"] = cmb.GetObservedShape(); + // pre_yields["total"].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb.GetObservedRate())); + // pre_yields["total"].back().setError(TMath::Power(cmb.GetObservedRate(), 0.5)); + // // Then fill total signal and total bkg hists + // std::cout << ">> Doing prefit: TotalBkg" << std::endl; + // pre_shapes_tot["TotalBkg"] = + // cmb.cp().backgrounds().GetShapeWithUncertainty(); + // pre_yields["total"].push_back(RooRealVar("yield_TotalBkg", "yield_TotalBkg", + // cmb.cp().backgrounds().GetRate())); + // pre_yields["total"].back().setError(cmb.cp().backgrounds().GetUncertainty()); + // std::cout << ">> Doing prefit: TotalSig" << std::endl; + // pre_shapes_tot["TotalSig"] = + // cmb.cp().signals().GetShapeWithUncertainty(); + // pre_yields["total"].push_back(RooRealVar("yield_TotalSig", "yield_TotalSig", + // cmb.cp().signals().GetRate())); + // pre_yields["total"].back().setError(cmb.cp().signals().GetUncertainty()); + // std::cout << ">> Doing prefit: TotalProcs" << std::endl; + // pre_shapes_tot["TotalProcs"] = + // cmb.cp().GetShapeWithUncertainty(); + // pre_yields["total"].push_back(RooRealVar("yield_TotalProcs", "yield_TotalProcs", + // cmb.cp().GetRate())); + // pre_yields["total"].back().setError(cmb.cp().GetUncertainty()); + + // if (datacard != "") { + // TH1F ref = cmb_card.cp().GetObservedShape(); + // for (auto & it : pre_shapes_tot) { + // it.second = ch::RestoreBinning(it.second, ref); + // } + // } + + // // Can write these straight into the output file + // outfile.cd(); + // for (auto& iter : pre_shapes_tot) { + // ch::WriteToTFile(&(iter.second), &outfile, "prefit/" + iter.first); + // } + // for (auto& yield: pre_yields["total"]){ + // helper.Form("%s/%s", "prefit", yield.GetName()); + // ch::WriteToTFile(&(yield), &outfile, helper.Data()); + // } + // map>().swap(pre_yields); } for (auto bin : bins) { - ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); - // This next line is a temporary fix for models with parameteric RooFit pdfs - // - we try and set the number of bins to evaluate the pdf to be the same as - // the number of bins in data - // cmb_bin.SetPdfBins(cmb_bin.GetObservedShape().GetNbinsX()); - - // Fill the data and process histograms - pre_shapes[bin]["data_obs"] = cmb_bin.GetObservedShape(); - for (auto proc : cmb_bin.process_set()) { - std::cout << ">> Doing prefit: " << bin << "," << proc << std::endl; - if (skip_proc_errs) { - pre_shapes[bin][proc] = - cmb_bin.cp().process({proc}).GetShape(); - } else { - pre_shapes[bin][proc] = - cmb_bin.cp().process({proc}).GetShapeWithUncertainty(); - } - } - // Create prefit shapes for merged processes - for (auto iter: merged_procs){ - // First element of the iterator is the name of the merged process - auto proc=iter.first; - std::cout << ">> Doing prefit: " << bin << "," << proc << std::endl; - // Second element is the regex expression for the processes - // that are to be merged - auto proc_regex = iter.second; - auto cmb_proc = cmb_bin.cp().process({proc_regex}); - // First check for matches - if (cmb_proc.process_set().size() == 0){ - std::cout << ">> WARNING: found no processes matching " << proc << std::endl; - continue; - } - if (skip_proc_errs) { - pre_shapes[bin][proc] = - cmb_proc.GetShape(); - } else { - pre_shapes[bin][proc] = - cmb_proc.GetShapeWithUncertainty(); - } - } - - // The fill total signal and total bkg hists - std::cout << ">> Doing prefit: " << bin << "," << "TotalBkg" << std::endl; - pre_shapes[bin]["TotalBkg"] = - cmb_bin.cp().backgrounds().GetShapeWithUncertainty(); - std::cout << ">> Doing prefit: " << bin << "," << "TotalSig" << std::endl; - pre_shapes[bin]["TotalSig"] = - cmb_bin.cp().signals().GetShapeWithUncertainty(); - std::cout << ">> Doing prefit: " << bin << "," << "TotalProcs" << std::endl; - pre_shapes[bin]["TotalProcs"] = - cmb_bin.cp().GetShapeWithUncertainty(); - - - if (datacard != "") { - TH1F ref = cmb_card.cp().bin({bin}).GetObservedShape(); - for (auto & it : pre_shapes[bin]) { - it.second = ch::RestoreBinning(it.second, ref); - } - } - - for (auto const& rbin : reverse_bins_) { - if (rbin != bin) continue; - auto & hists = pre_shapes[bin]; - for (auto it = hists.begin(); it != hists.end(); ++it) { - ReverseBins(it->second); - } - } - // Can write these straight into the output file - outfile.cd(); - for (auto& iter : pre_shapes[bin]) { - ch::WriteToTFile(&(iter.second), &outfile, bin + "_prefit/" + iter.first); - } - } - - // Print out the relative uncert. on the bkg - if (factors) { - cout << boost::format("%-25s %-32s\n") % "Bin" % - "Total relative bkg uncert. (prefit)"; - cout << string(58, '-') << "\n"; - for (auto bin : bins) { - ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); - double rate = cmb_bin.cp().backgrounds().GetRate(); - double err = cmb_bin.cp().backgrounds().GetUncertainty(); - cout << boost::format("%-25s %-10.5f") % bin % - (rate > 0. ? (err / rate) : 0.) << std::endl; - } + pre_shapes[bin] = map(); + pre_yields[bin] = vector(); + do_complete_set( + cmb.cp().bin({bin}), + outfile, + bin+"_prefit", + pre_shapes[bin], + pre_yields[bin], + merged_procs, + skip_procs, + cmb_card.cp().bin({bin}), + datacard, + skip_proc_errs, + factors, + std::find(reverse_bins_.begin(), reverse_bins_.end(), bin) != reverse_bins_.end()); } } @@ -347,124 +463,234 @@ int main(int argc, char* argv[]) { // Calculate the post-fit fractional background uncertainty in each bin map> post_shapes; + map> post_yields; map post_yield_cov; map post_yield_cor; map post_shapes_tot; + - if(total_shapes){ - post_shapes_tot["data_obs"] = cmb.GetObservedShape(); - // Fill the total sig. and total bkg. hists - auto cmb_bkgs = cmb.cp().backgrounds(); - auto cmb_sigs = cmb.cp().signals(); - std::cout << ">> Doing postfit: TotalBkg" << std::endl; - post_shapes_tot["TotalBkg"] = - no_sampling ? cmb_bkgs.GetShapeWithUncertainty() - : cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose); - std::cout << ">> Doing postfit: TotalSig" << std::endl; - post_shapes_tot["TotalSig"] = - no_sampling ? cmb_sigs.GetShapeWithUncertainty() - : cmb_sigs.GetShapeWithUncertainty(res, samples, verbose); - std::cout << ">> Doing postfit: TotalProcs" << std::endl; - post_shapes_tot["TotalProcs"] = - no_sampling ? cmb.cp().GetShapeWithUncertainty() - : cmb.cp().GetShapeWithUncertainty(res, samples, verbose); - - if (datacard != "") { - TH1F ref = cmb_card.cp().GetObservedShape(); - for (auto & it : post_shapes_tot) { - it.second = ch::RestoreBinning(it.second, ref); - } - } - outfile.cd(); - // Write the post-fit histograms - for (auto & iter : post_shapes_tot) { - ch::WriteToTFile(&(iter.second), &outfile, - "postfit/" + iter.first); - } + if(total_shapes){ + post_yields["total"] = vector(); + post_shapes["total"] = map(); + do_complete_set( + cmb, + outfile, + "total_postfit", + post_shapes["total"], + post_yields["total"], + merged_procs, + skip_procs, + cmb_card, + datacard, + skip_proc_errs, + factors, + false, + sampling, + samples, + &res, + verbose); + // post_shapes["data_obs"] = cmb.GetObservedShape(); + // post_yields["total"].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb.GetObservedRate())); + // post_yields["total"].back().setError(TMath::Power(cmb.GetObservedRate(), 0.5)); + // // Fill the total sig. and total bkg. hists + // auto cmb_bkgs = cmb.cp().backgrounds(); + // auto cmb_sigs = cmb.cp().signals(); + // std::cout << ">> Doing postfit: TotalBkg" << std::endl; + // post_shapes_tot["TotalBkg"] = + // sampling ? cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_bkgs.GetShapeWithUncertainty(); + + // post_yields["total"].push_back(RooRealVar("yield_TotalBkg", "yield_TotalBkg", + // cmb_bkgs.cp().GetRate())); + // post_yields["total"].back().setError( + // sampling ? cmb_bkgs.GetUncertainty(res, samples) + // : cmb_bkgs.GetUncertainty()); + // std::cout << ">> Doing postfit: TotalSig" << std::endl; + // post_shapes_tot["TotalSig"] = + // sampling ? cmb_sigs.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_sigs.GetShapeWithUncertainty(); + // post_yields["total"].push_back(RooRealVar("yield_TotalSig", "yield_TotalSig", + // cmb_sigs.cp().GetRate())); + // post_yields["total"].back().setError( + // sampling ? cmb_sigs.GetUncertainty(res, samples) + // : cmb_sigs.GetUncertainty()); + // std::cout << ">> Doing postfit: TotalProcs" << std::endl; + // post_shapes_tot["TotalProcs"] = + // sampling ? cmb.cp().GetShapeWithUncertainty(res, samples, verbose) + // : cmb.cp().GetShapeWithUncertainty(); + + // post_yields["total"].push_back(RooRealVar("yield_TotalProcs", "yield_TotalProcs", + // cmb.cp().GetRate())); + // post_yields["total"].back().setError( + // sampling ? cmb.GetUncertainty(res, samples) + // : cmb.GetUncertainty()); + + // if (datacard != "") { + // TH1F ref = cmb_card.cp().GetObservedShape(); + // for (auto & it : post_shapes_tot) { + // it.second = ch::RestoreBinning(it.second, ref); + // } + // } + + // outfile.cd(); + // // Write the post-fit histograms + // for (auto & iter : post_shapes_tot) { + // ch::WriteToTFile(&(iter.second), &outfile, + // "postfit/" + iter.first); + // } + // for (auto& yield: post_yields["total"]){ + // helper.Form("%s/%s", "postfit", yield.GetName()); + // ch::WriteToTFile(&(yield), &outfile, helper.Data()); + // } + // map>().swap(post_yields); } for (auto bin : bins) { - ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); - post_shapes[bin]["data_obs"] = cmb_bin.GetObservedShape(); - for (auto proc : cmb_bin.process_set()) { - auto cmb_proc = cmb_bin.cp().process({proc}); - // Method to get the shape uncertainty depends on whether we are using - // the sampling method or the "wrong" method (assumes no correlations) - std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; - if (skip_proc_errs) { - post_shapes[bin][proc] = cmb_proc.GetShape(); - } else { - post_shapes[bin][proc] = - no_sampling ? cmb_proc.GetShapeWithUncertainty() - : cmb_proc.GetShapeWithUncertainty(res, samples, verbose); - } - } - - // Generate postfit distributions for merged processes - for (auto iter: merged_procs){ - auto proc=iter.first; - std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; - auto proc_regex = iter.second; - auto cmb_proc = cmb_bin.cp().process({proc_regex}); - if (cmb_proc.process_set().size() == 0){ - std::cout << ">> WARNING: found no processes matching " << proc << std::endl; - continue; - } - if (skip_proc_errs) { - post_shapes[bin][proc] = cmb_proc.GetShape(); - } else { - post_shapes[bin][proc] = - sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) - : cmb_proc.GetShapeWithUncertainty(); - } - } - - if (!no_sampling && covariance) { + // ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); + // post_yields[bin] = vector(); + // post_shapes[bin]["data_obs"] = cmb_bin.GetObservedShape(); + // post_yields[bin].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb_bin.GetObservedRate())); + // post_yields[bin].back().setError(TMath::Power(cmb_bin.GetObservedRate(), 0.5)); + // for (auto proc : cmb_bin.process_set()) { + // auto cmb_proc = cmb_bin.cp().process({proc}); + // // Method to get the shape uncertainty depends on whether we are using + // // the sampling method or the "wrong" method (assumes no correlations) + // std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; + // if (skip_proc_errs) { + // post_shapes[bin][proc] = cmb_proc.GetShape(); + // } else { + // post_shapes[bin][proc] = + // sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_proc.GetShapeWithUncertainty(); + // } + // helper.Form("%s_%s", "yield", proc.c_str()); + // post_yields[bin].push_back(RooRealVar(helper, helper, + // cmb_proc.cp().GetRate())); + // post_yields[bin].back().setError( + // sampling ? cmb_proc.GetUncertainty(res, samples) + // : cmb_proc.GetUncertainty()); + // } + // for (auto iter: merged_procs){ + // auto proc=iter.first; + // std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; + // auto proc_regex = iter.second; + // auto cmb_proc = cmb_bin.cp().process({proc_regex}); + // if (cmb_proc.process_set().size() == 0){ + // std::cout << ">> WARNING: found no processes matching " << proc << std::endl; + // continue; + // } + // if (skip_proc_errs) { + // post_shapes[bin][proc] = cmb_proc.GetShape(); + // } else { + // post_shapes[bin][proc] = + // sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_proc.GetShapeWithUncertainty(); + // } + // helper.Form("%s_%s", "yield", proc.c_str()); + // post_yields[bin].push_back(RooRealVar(helper, helper, + // cmb_proc.cp().GetRate())); + // post_yields[bin].back().setError( + // sampling ? cmb_proc.GetUncertainty(res, samples) + // : cmb_proc.GetUncertainty()); + // } + post_yields[bin] = vector(); + post_shapes[bin] = map(); + auto cmb_bin = cmb.cp().bin({bin}); + do_complete_set( + cmb_bin, + outfile, + bin+"_postfit", + post_shapes[bin], + post_yields[bin], + merged_procs, + skip_procs, + cmb_card.cp().bin({bin}), + datacard, + skip_proc_errs, + factors, + std::find(reverse_bins_.begin(), reverse_bins_.end(), bin) != reverse_bins_.end(), + sampling, + samples, + &res, + verbose); + if (sampling && covariance) { post_yield_cov[bin] = cmb_bin.GetRateCovariance(res, samples); post_yield_cor[bin] = cmb_bin.GetRateCorrelation(res, samples); } // Fill the total sig. and total bkg. hists - auto cmb_bkgs = cmb_bin.cp().backgrounds(); - auto cmb_sigs = cmb_bin.cp().signals(); - std::cout << ">> Doing postfit: " << bin << "," << "TotalBkg" << std::endl; - post_shapes[bin]["TotalBkg"] = - no_sampling ? cmb_bkgs.GetShapeWithUncertainty() - : cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose); - std::cout << ">> Doing postfit: " << bin << "," << "TotalSig" << std::endl; - post_shapes[bin]["TotalSig"] = - no_sampling ? cmb_sigs.GetShapeWithUncertainty() - : cmb_sigs.GetShapeWithUncertainty(res, samples, verbose); - std::cout << ">> Doing postfit: " << bin << "," << "TotalProcs" << std::endl; - post_shapes[bin]["TotalProcs"] = - no_sampling ? cmb_bin.cp().GetShapeWithUncertainty() - : cmb_bin.cp().GetShapeWithUncertainty(res, samples, verbose); - - if (datacard != "") { - TH1F ref = cmb_card.cp().bin({bin}).GetObservedShape(); - for (auto & it : post_shapes[bin]) { - it.second = ch::RestoreBinning(it.second, ref); - } - } + // auto cmb_bkgs = cmb_bin.cp().backgrounds(); + // auto cmb_sigs = cmb_bin.cp().signals(); + // std::cout << ">> Doing postfit: " << bin << "," << "TotalBkg" << std::endl; + // post_shapes[bin]["TotalBkg"] = + // sampling ? cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_bkgs.GetShapeWithUncertainty(); + // helper.Form("%s_%s", "yield", "TotalBkg"); + // std::cout << "calculating yields for TotalBkg\n"; + // post_yields[bin].push_back(RooRealVar(helper, helper, + // cmb_bkgs.cp().GetRate())); + // post_yields[bin].back().setError( + // sampling ? cmb_bkgs.GetUncertainty(res, samples) + // : cmb_bkgs.GetUncertainty()); + // std::cout << ">> Doing postfit: " << bin << "," << "TotalSig" << std::endl; + // auto signal_names = cmb_sigs.process_set(); + // for(auto name : signal_names){ + // std::cout << " " << name; + // } + // std::cout << std::endl; + // if(sampling) std::cout << "will generate " << samples << " toys" << std::endl; + // post_shapes[bin]["TotalSig"] = + // sampling ? cmb_sigs.GetShapeWithUncertainty(res, samples, verbose) + // : cmb_sigs.GetShapeWithUncertainty(); + // helper.Form("%s_%s", "yield", "TotalSig"); + // post_yields[bin].push_back(RooRealVar(helper, helper, + // cmb_sigs.cp().GetRate())); + // post_yields[bin].back().setError( + // sampling ? cmb_sigs.GetUncertainty(res, samples) + // : cmb_sigs.GetUncertainty()); + + // std::cout << ">> Doing postfit: " << bin << "," << "TotalProcs" << std::endl; + // post_shapes[bin]["TotalProcs"] = + // sampling ? cmb_bin.cp().GetShapeWithUncertainty(res, samples, verbose) + // : cmb_bin.cp().GetShapeWithUncertainty(); + + // helper.Form("%s_%s", "yield", "TotalProcs"); + // post_yields[bin].push_back(RooRealVar(helper, helper, + // cmb_bin.cp().GetRate())); + // post_yields[bin].back().setError( + // sampling ? cmb_bin.GetUncertainty(res, samples) + // : cmb_bin.GetUncertainty()); + + // if (datacard != "") { + // TH1F ref = cmb_card.cp().bin({bin}).GetObservedShape(); + // for (auto & it : post_shapes[bin]) { + // it.second = ch::RestoreBinning(it.second, ref); + // } + // } - outfile.cd(); // Write the post-fit histograms - for (auto const& rbin : reverse_bins_) { - if (rbin != bin) continue; - std::cout << ">> reversing hists in bin " << bin << "\n"; - auto & hists = post_shapes[bin]; - for (auto it = hists.begin(); it != hists.end(); ++it) { - ReverseBins(it->second); - } - } + // for (auto const& rbin : reverse_bins_) { + // if (rbin != bin) continue; + // std::cout << ">> reversing hists in bin " << bin << "\n"; + // auto & hists = post_shapes[bin]; + // for (auto it = hists.begin(); it != hists.end(); ++it) { + // ReverseBins(it->second); + // } + // } - for (auto & iter : post_shapes[bin]) { - ch::WriteToTFile(&(iter.second), &outfile, - bin + "_postfit/" + iter.first); - } + outfile.cd(); + + // for (auto & iter : post_shapes[bin]) { + // ch::WriteToTFile(&(iter.second), &outfile, + // bin + "_postfit/" + iter.first); + // } + // for (auto& yield: post_yields[bin]){ + // helper.Form("%s_%s/%s", bin.c_str(), "postfit" , yield.GetName()); + // ch::WriteToTFile(&(yield), &outfile, helper.Data()); + // } for (auto & iter : post_yield_cov) { ch::WriteToTFile(&(iter.second), &outfile, iter.first+"_cov"); @@ -476,19 +702,19 @@ int main(int argc, char* argv[]) { } - if (factors) { - cout << boost::format("\n%-25s %-32s\n") % "Bin" % - "Total relative bkg uncert. (postfit)"; - cout << string(58, '-') << "\n"; - for (auto bin : bins) { - ch::CombineHarvester cmb_bkgs = cmb.cp().bin({bin}).backgrounds(); - double rate = cmb_bkgs.GetRate(); - double err = no_sampling ? cmb_bkgs.GetUncertainty() - : cmb_bkgs.GetUncertainty(res, samples); - cout << boost::format("%-25s %-10.5f") % bin % - (rate > 0. ? (err / rate) : 0.) << std::endl; - } - } + // if (factors) { + // cout << boost::format("\n%-25s %-32s\n") % "Bin" % + // "Total relative bkg uncert. (postfit)"; + // cout << string(58, '-') << "\n"; + // for (auto bin : bins) { + // ch::CombineHarvester cmb_bkgs = cmb.cp().bin({bin}).backgrounds(); + // double rate = cmb_bkgs.GetRate(); + // double err = sampling ? cmb_bkgs.GetUncertainty(res, samples) + // : cmb_bkgs.GetUncertainty(); + // cout << boost::format("%-25s %-10.5f") % bin % + // (rate > 0. ? (err / rate) : 0.) << std::endl; + // } + // } // As we calculate the post-fit yields can also print out the post/pre scale // factors From 9e29eec044aa96bf8154d0074a910ab370d723df Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 7 Mar 2023 18:15:26 +0100 Subject: [PATCH 06/13] included previously overwritten older changes in PostFitShapesFromWorkspace.cpp --- .../bin/PostFitShapesFromWorkspace.cpp | 274 ++---------------- 1 file changed, 29 insertions(+), 245 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index e022a66005e..46687e9af7f 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -38,7 +38,7 @@ void load_shapes_and_yields( vector& yield_container, const std::string& proc, const bool& skip_proc_errs=true, - const bool& sampling=false, + const bool& no_sampling=true, const int& samples=0, RooFitResult* res=NULL, const bool& verbose=false){ @@ -48,14 +48,14 @@ void load_shapes_and_yields( cmb.GetShape(); } else { shape_container[proc] = - sampling ? cmb.GetShapeWithUncertainty(*res, samples, verbose) - : cmb.GetShapeWithUncertainty(); + no_sampling ? cmb.GetShapeWithUncertainty() + : cmb.GetShapeWithUncertainty(*res, samples); } helper.Form("%s_%s", "yield" , proc.c_str()); yield_container.push_back(RooRealVar(helper, helper, cmb.GetRate())); - yield_container.back().setError(sampling ? cmb.GetUncertainty(*res, samples) - : cmb.GetUncertainty()); + yield_container.back().setError(no_sampling ? cmb.GetUncertainty() + : cmb.GetUncertainty(*res, samples)); } void do_complete_set( @@ -71,7 +71,7 @@ void do_complete_set( const bool& skip_proc_errs=true, const bool& factors=false, const bool& reverse_bins=false, - const bool& sampling=false, + const bool& no_sampling=true, const int& samples=0, RooFitResult* res=NULL, const bool& verbose=false @@ -90,7 +90,7 @@ void do_complete_set( for (auto proc : cmb_bin.process_set()) { if(skip_procs.find(proc) == skip_procs.end()){ std::cout << ">> Doing " << output_prefix.c_str() << "," << proc << std::endl; - load_shapes_and_yields(cmb_bin.cp().process({proc}), shape_container, yield_container, proc, skip_proc_errs, sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().process({proc}), shape_container, yield_container, proc, skip_proc_errs, no_sampling, samples, res, verbose); } else{ std::cout << ">> Skipping processes " << proc << std::endl; @@ -105,12 +105,12 @@ void do_complete_set( std::cout << ">> WARNING: found no processes matching " << proc << std::endl; continue; } - load_shapes_and_yields(cmb_proc, shape_container, yield_container, proc, skip_proc_errs, sampling, samples, res, verbose); + load_shapes_and_yields(cmb_proc, shape_container, yield_container, proc, skip_proc_errs, no_sampling, samples, res, verbose); } // The fill total signal and total bkg hists std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalBkg" << std::endl; - load_shapes_and_yields(cmb_bin.cp().backgrounds(), shape_container, yield_container, "TotalBkg", false, sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().backgrounds(), shape_container, yield_container, "TotalBkg", false, no_sampling, samples, res, verbose); // Print out the relative uncert. on the bkg if (factors) { @@ -125,10 +125,10 @@ void do_complete_set( } std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalSig" << std::endl; - load_shapes_and_yields(cmb_bin.cp().signals(), shape_container, yield_container, "TotalSig", false, sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().signals(), shape_container, yield_container, "TotalSig", false, no_sampling, samples, res, verbose); std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalProcs" << std::endl; - load_shapes_and_yields(cmb_bin, shape_container, yield_container, "TotalProcs", false, sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin, shape_container, yield_container, "TotalProcs", false, no_sampling, samples, res, verbose); if (datacard != "") { @@ -165,6 +165,7 @@ int main(int argc, char* argv[]) { string mass = ""; bool postfit = false; bool sampling = false; + bool no_sampling = false; string output = ""; bool factors = false; unsigned samples = 500; @@ -213,7 +214,10 @@ int main(int argc, char* argv[]) { "Create post-fit histograms in addition to pre-fit") ("sampling", po::value(&sampling)->default_value(sampling)->implicit_value(true), - "Use the cov. matrix sampling method for the post-fit uncertainty") + "Use the cov. matrix sampling method for the post-fit uncertainty (deprecated, this is the default)") + ("no-sampling", + po::value(&no_sampling)->default_value(no_sampling)->implicit_value(true), + "Do not use the cov. matrix sampling method for the post-fit uncertainty") ("samples", po::value(&samples)->default_value(samples), "Number of samples to make in each evaluate call") @@ -245,22 +249,22 @@ int main(int argc, char* argv[]) { ("skip-procs,s", po::value>(&skip_procs_)->multitoken(), "Skip these processes. Regex expression allowed. Format: 'expression'. Can be called multiple times"); - if (sampling && !postfit) { - throw logic_error( - "Can't sample the fit covariance matrix for pre-fit!"); - } + // if (!no_sampling && !postfit) { + // throw logic_error( + // "Can't sample the fit covariance matrix for pre-fit!"); + // } po::variables_map vm; // First check if the user has set the "--help" or "-h" option, and if so - // just prin the usage information and quit + // just print the usage information and quit po::store(po::command_line_parser(argc, argv) .options(help_config).allow_unregistered().run(), vm); po::notify(vm); if (vm.count("help")) { cout << config << "\nExample usage:\n"; cout << "PostFitShapesFromWorkspace.root -d htt_mt_125.txt -w htt_mt_125.root -o htt_mt_125_shapes.root -m 125 " - "-f mlfit.root:fit_s --postfit --sampling --print\n"; + "-f mlfit.root:fit_s --postfit --print\n"; return 1; } @@ -268,6 +272,10 @@ int main(int argc, char* argv[]) { po::store(po::command_line_parser(argc, argv).options(config).run(), vm); po::notify(vm); + if (sampling) { + std::cout<<"WARNING: the default behaviour of PostFitShapesFromWorkspace is to use the covariance matrix sampling method for the post-fit uncertainty. The option --sampling is deprecated and will be removed in future versions of CombineHarvester"<(gDirectory->Get("w")); @@ -392,47 +400,6 @@ int main(int argc, char* argv[]) { pre_shapes["total"] = map(); pre_yields["total"] = vector(); do_complete_set(cmb, outfile, "total_prefit", pre_shapes["total"], pre_yields["total"], merged_procs, skip_procs, cmb_card, datacard, skip_proc_errs, factors); - // pre_yields["total"] = vector(); - // pre_shapes_tot["data_obs"] = cmb.GetObservedShape(); - // pre_yields["total"].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb.GetObservedRate())); - // pre_yields["total"].back().setError(TMath::Power(cmb.GetObservedRate(), 0.5)); - // // Then fill total signal and total bkg hists - // std::cout << ">> Doing prefit: TotalBkg" << std::endl; - // pre_shapes_tot["TotalBkg"] = - // cmb.cp().backgrounds().GetShapeWithUncertainty(); - // pre_yields["total"].push_back(RooRealVar("yield_TotalBkg", "yield_TotalBkg", - // cmb.cp().backgrounds().GetRate())); - // pre_yields["total"].back().setError(cmb.cp().backgrounds().GetUncertainty()); - // std::cout << ">> Doing prefit: TotalSig" << std::endl; - // pre_shapes_tot["TotalSig"] = - // cmb.cp().signals().GetShapeWithUncertainty(); - // pre_yields["total"].push_back(RooRealVar("yield_TotalSig", "yield_TotalSig", - // cmb.cp().signals().GetRate())); - // pre_yields["total"].back().setError(cmb.cp().signals().GetUncertainty()); - // std::cout << ">> Doing prefit: TotalProcs" << std::endl; - // pre_shapes_tot["TotalProcs"] = - // cmb.cp().GetShapeWithUncertainty(); - // pre_yields["total"].push_back(RooRealVar("yield_TotalProcs", "yield_TotalProcs", - // cmb.cp().GetRate())); - // pre_yields["total"].back().setError(cmb.cp().GetUncertainty()); - - // if (datacard != "") { - // TH1F ref = cmb_card.cp().GetObservedShape(); - // for (auto & it : pre_shapes_tot) { - // it.second = ch::RestoreBinning(it.second, ref); - // } - // } - - // // Can write these straight into the output file - // outfile.cd(); - // for (auto& iter : pre_shapes_tot) { - // ch::WriteToTFile(&(iter.second), &outfile, "prefit/" + iter.first); - // } - // for (auto& yield: pre_yields["total"]){ - // helper.Form("%s/%s", "prefit", yield.GetName()); - // ch::WriteToTFile(&(yield), &outfile, helper.Data()); - // } - // map>().swap(pre_yields); } for (auto bin : bins) { pre_shapes[bin] = map(); @@ -487,115 +454,14 @@ int main(int argc, char* argv[]) { skip_proc_errs, factors, false, - sampling, + no_sampling, samples, &res, verbose); - // post_shapes["data_obs"] = cmb.GetObservedShape(); - // post_yields["total"].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb.GetObservedRate())); - // post_yields["total"].back().setError(TMath::Power(cmb.GetObservedRate(), 0.5)); - // // Fill the total sig. and total bkg. hists - // auto cmb_bkgs = cmb.cp().backgrounds(); - // auto cmb_sigs = cmb.cp().signals(); - // std::cout << ">> Doing postfit: TotalBkg" << std::endl; - // post_shapes_tot["TotalBkg"] = - // sampling ? cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_bkgs.GetShapeWithUncertainty(); - - // post_yields["total"].push_back(RooRealVar("yield_TotalBkg", "yield_TotalBkg", - // cmb_bkgs.cp().GetRate())); - // post_yields["total"].back().setError( - // sampling ? cmb_bkgs.GetUncertainty(res, samples) - // : cmb_bkgs.GetUncertainty()); - // std::cout << ">> Doing postfit: TotalSig" << std::endl; - // post_shapes_tot["TotalSig"] = - // sampling ? cmb_sigs.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_sigs.GetShapeWithUncertainty(); - // post_yields["total"].push_back(RooRealVar("yield_TotalSig", "yield_TotalSig", - // cmb_sigs.cp().GetRate())); - // post_yields["total"].back().setError( - // sampling ? cmb_sigs.GetUncertainty(res, samples) - // : cmb_sigs.GetUncertainty()); - // std::cout << ">> Doing postfit: TotalProcs" << std::endl; - // post_shapes_tot["TotalProcs"] = - // sampling ? cmb.cp().GetShapeWithUncertainty(res, samples, verbose) - // : cmb.cp().GetShapeWithUncertainty(); - - // post_yields["total"].push_back(RooRealVar("yield_TotalProcs", "yield_TotalProcs", - // cmb.cp().GetRate())); - // post_yields["total"].back().setError( - // sampling ? cmb.GetUncertainty(res, samples) - // : cmb.GetUncertainty()); - - // if (datacard != "") { - // TH1F ref = cmb_card.cp().GetObservedShape(); - // for (auto & it : post_shapes_tot) { - // it.second = ch::RestoreBinning(it.second, ref); - // } - // } - - // outfile.cd(); - // // Write the post-fit histograms - // for (auto & iter : post_shapes_tot) { - // ch::WriteToTFile(&(iter.second), &outfile, - // "postfit/" + iter.first); - // } - // for (auto& yield: post_yields["total"]){ - // helper.Form("%s/%s", "postfit", yield.GetName()); - // ch::WriteToTFile(&(yield), &outfile, helper.Data()); - // } - // map>().swap(post_yields); } for (auto bin : bins) { - // ch::CombineHarvester cmb_bin = cmb.cp().bin({bin}); - // post_yields[bin] = vector(); - // post_shapes[bin]["data_obs"] = cmb_bin.GetObservedShape(); - // post_yields[bin].push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb_bin.GetObservedRate())); - // post_yields[bin].back().setError(TMath::Power(cmb_bin.GetObservedRate(), 0.5)); - // for (auto proc : cmb_bin.process_set()) { - // auto cmb_proc = cmb_bin.cp().process({proc}); - // // Method to get the shape uncertainty depends on whether we are using - // // the sampling method or the "wrong" method (assumes no correlations) - // std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; - // if (skip_proc_errs) { - // post_shapes[bin][proc] = cmb_proc.GetShape(); - // } else { - // post_shapes[bin][proc] = - // sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_proc.GetShapeWithUncertainty(); - // } - // helper.Form("%s_%s", "yield", proc.c_str()); - // post_yields[bin].push_back(RooRealVar(helper, helper, - // cmb_proc.cp().GetRate())); - // post_yields[bin].back().setError( - // sampling ? cmb_proc.GetUncertainty(res, samples) - // : cmb_proc.GetUncertainty()); - // } - // for (auto iter: merged_procs){ - // auto proc=iter.first; - // std::cout << ">> Doing postfit: " << bin << "," << proc << std::endl; - // auto proc_regex = iter.second; - // auto cmb_proc = cmb_bin.cp().process({proc_regex}); - // if (cmb_proc.process_set().size() == 0){ - // std::cout << ">> WARNING: found no processes matching " << proc << std::endl; - // continue; - // } - // if (skip_proc_errs) { - // post_shapes[bin][proc] = cmb_proc.GetShape(); - // } else { - // post_shapes[bin][proc] = - // sampling ? cmb_proc.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_proc.GetShapeWithUncertainty(); - // } - // helper.Form("%s_%s", "yield", proc.c_str()); - // post_yields[bin].push_back(RooRealVar(helper, helper, - // cmb_proc.cp().GetRate())); - // post_yields[bin].back().setError( - // sampling ? cmb_proc.GetUncertainty(res, samples) - // : cmb_proc.GetUncertainty()); - // } post_yields[bin] = vector(); post_shapes[bin] = map(); auto cmb_bin = cmb.cp().bin({bin}); @@ -612,85 +478,17 @@ int main(int argc, char* argv[]) { skip_proc_errs, factors, std::find(reverse_bins_.begin(), reverse_bins_.end(), bin) != reverse_bins_.end(), - sampling, + no_sampling, samples, &res, verbose); - if (sampling && covariance) { + if (!no_sampling && covariance) { post_yield_cov[bin] = cmb_bin.GetRateCovariance(res, samples); post_yield_cor[bin] = cmb_bin.GetRateCorrelation(res, samples); } - // Fill the total sig. and total bkg. hists - // auto cmb_bkgs = cmb_bin.cp().backgrounds(); - // auto cmb_sigs = cmb_bin.cp().signals(); - // std::cout << ">> Doing postfit: " << bin << "," << "TotalBkg" << std::endl; - // post_shapes[bin]["TotalBkg"] = - // sampling ? cmb_bkgs.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_bkgs.GetShapeWithUncertainty(); - // helper.Form("%s_%s", "yield", "TotalBkg"); - // std::cout << "calculating yields for TotalBkg\n"; - // post_yields[bin].push_back(RooRealVar(helper, helper, - // cmb_bkgs.cp().GetRate())); - // post_yields[bin].back().setError( - // sampling ? cmb_bkgs.GetUncertainty(res, samples) - // : cmb_bkgs.GetUncertainty()); - // std::cout << ">> Doing postfit: " << bin << "," << "TotalSig" << std::endl; - // auto signal_names = cmb_sigs.process_set(); - // for(auto name : signal_names){ - // std::cout << " " << name; - // } - // std::cout << std::endl; - // if(sampling) std::cout << "will generate " << samples << " toys" << std::endl; - // post_shapes[bin]["TotalSig"] = - // sampling ? cmb_sigs.GetShapeWithUncertainty(res, samples, verbose) - // : cmb_sigs.GetShapeWithUncertainty(); - // helper.Form("%s_%s", "yield", "TotalSig"); - // post_yields[bin].push_back(RooRealVar(helper, helper, - // cmb_sigs.cp().GetRate())); - // post_yields[bin].back().setError( - // sampling ? cmb_sigs.GetUncertainty(res, samples) - // : cmb_sigs.GetUncertainty()); - - // std::cout << ">> Doing postfit: " << bin << "," << "TotalProcs" << std::endl; - // post_shapes[bin]["TotalProcs"] = - // sampling ? cmb_bin.cp().GetShapeWithUncertainty(res, samples, verbose) - // : cmb_bin.cp().GetShapeWithUncertainty(); - - // helper.Form("%s_%s", "yield", "TotalProcs"); - // post_yields[bin].push_back(RooRealVar(helper, helper, - // cmb_bin.cp().GetRate())); - // post_yields[bin].back().setError( - // sampling ? cmb_bin.GetUncertainty(res, samples) - // : cmb_bin.GetUncertainty()); - - // if (datacard != "") { - // TH1F ref = cmb_card.cp().bin({bin}).GetObservedShape(); - // for (auto & it : post_shapes[bin]) { - // it.second = ch::RestoreBinning(it.second, ref); - // } - // } - - // Write the post-fit histograms - - // for (auto const& rbin : reverse_bins_) { - // if (rbin != bin) continue; - // std::cout << ">> reversing hists in bin " << bin << "\n"; - // auto & hists = post_shapes[bin]; - // for (auto it = hists.begin(); it != hists.end(); ++it) { - // ReverseBins(it->second); - // } - // } outfile.cd(); - // for (auto & iter : post_shapes[bin]) { - // ch::WriteToTFile(&(iter.second), &outfile, - // bin + "_postfit/" + iter.first); - // } - // for (auto& yield: post_yields[bin]){ - // helper.Form("%s_%s/%s", bin.c_str(), "postfit" , yield.GetName()); - // ch::WriteToTFile(&(yield), &outfile, helper.Data()); - // } for (auto & iter : post_yield_cov) { ch::WriteToTFile(&(iter.second), &outfile, iter.first+"_cov"); @@ -702,20 +500,6 @@ int main(int argc, char* argv[]) { } - // if (factors) { - // cout << boost::format("\n%-25s %-32s\n") % "Bin" % - // "Total relative bkg uncert. (postfit)"; - // cout << string(58, '-') << "\n"; - // for (auto bin : bins) { - // ch::CombineHarvester cmb_bkgs = cmb.cp().bin({bin}).backgrounds(); - // double rate = cmb_bkgs.GetRate(); - // double err = sampling ? cmb_bkgs.GetUncertainty(res, samples) - // : cmb_bkgs.GetUncertainty(); - // cout << boost::format("%-25s %-10.5f") % bin % - // (rate > 0. ? (err / rate) : 0.) << std::endl; - // } - // } - // As we calculate the post-fit yields can also print out the post/pre scale // factors if (factors && postfit) { From 31354ddc03b686b5d9da83228af2e2417ce03af0 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Tue, 7 Mar 2023 19:47:41 +0100 Subject: [PATCH 07/13] make yield generation optional --- .../bin/PostFitShapesFromWorkspace.cpp | 72 ++++++++++--------- 1 file changed, 40 insertions(+), 32 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index 46687e9af7f..e0802db8325 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -35,12 +35,12 @@ void ReverseBins(TH1F & h) { void load_shapes_and_yields( ch::CombineHarvester& cmb, map& shape_container, - vector& yield_container, const std::string& proc, + vector* yield_container=nullptr, const bool& skip_proc_errs=true, const bool& no_sampling=true, const int& samples=0, - RooFitResult* res=NULL, + RooFitResult* res=nullptr, const bool& verbose=false){ TString helper; if (skip_proc_errs) { @@ -51,11 +51,14 @@ void load_shapes_and_yields( no_sampling ? cmb.GetShapeWithUncertainty() : cmb.GetShapeWithUncertainty(*res, samples); } - helper.Form("%s_%s", "yield" , proc.c_str()); - yield_container.push_back(RooRealVar(helper, helper, - cmb.GetRate())); - yield_container.back().setError(no_sampling ? cmb.GetUncertainty() - : cmb.GetUncertainty(*res, samples)); + // fill the yield container if there is one + if(yield_container != nullptr){ + helper.Form("%s_%s", "yield" , proc.c_str()); + yield_container->push_back(RooRealVar(helper, helper, + cmb.GetRate())); + yield_container->back().setError(no_sampling ? cmb.GetUncertainty() + : cmb.GetUncertainty(*res, samples)); + } } void do_complete_set( @@ -63,10 +66,10 @@ void do_complete_set( TFile& outfile, const std::string& output_prefix, map& shape_container, - vector& yield_container, const std::map& merged_procs, const std::set skip_procs, ch::CombineHarvester& cmb_card, + bool save_yields=true, const std::string& datacard="", const bool& skip_proc_errs=true, const bool& factors=false, @@ -83,14 +86,19 @@ void do_complete_set( // the number of bins in data // cmb_bin.SetPdfBins(cmb_bin.GetObservedShape().GetNbinsX()); + // initialize container for yields if necessary + std::vector* yield_container = nullptr; + if(save_yields) yield_container = new std::vector(); // Fill the data and process histograms shape_container["data_obs"] = cmb_bin.GetObservedShape(); - yield_container.push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb_bin.GetObservedRate())); - yield_container.back().setError(TMath::Power(cmb_bin.GetObservedRate(), 0.5)); + if(yield_container != nullptr){ + yield_container->push_back(RooRealVar("yield_data_obs", "yield_data_obs", cmb_bin.GetObservedRate())); + yield_container->back().setError(TMath::Power(cmb_bin.GetObservedRate(), 0.5)); + } for (auto proc : cmb_bin.process_set()) { if(skip_procs.find(proc) == skip_procs.end()){ std::cout << ">> Doing " << output_prefix.c_str() << "," << proc << std::endl; - load_shapes_and_yields(cmb_bin.cp().process({proc}), shape_container, yield_container, proc, skip_proc_errs, no_sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().process({proc}), shape_container, proc, yield_container, skip_proc_errs, no_sampling, samples, res, verbose); } else{ std::cout << ">> Skipping processes " << proc << std::endl; @@ -105,30 +113,30 @@ void do_complete_set( std::cout << ">> WARNING: found no processes matching " << proc << std::endl; continue; } - load_shapes_and_yields(cmb_proc, shape_container, yield_container, proc, skip_proc_errs, no_sampling, samples, res, verbose); + load_shapes_and_yields(cmb_proc, shape_container, proc, yield_container, skip_proc_errs, no_sampling, samples, res, verbose); } // The fill total signal and total bkg hists std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalBkg" << std::endl; - load_shapes_and_yields(cmb_bin.cp().backgrounds(), shape_container, yield_container, "TotalBkg", false, no_sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().backgrounds(), shape_container, "TotalBkg", yield_container, false, no_sampling, samples, res, verbose); // Print out the relative uncert. on the bkg - if (factors) { + if (factors && yield_container != nullptr) { cout << boost::format("%-25s %-32s\n") % "Bin" % "Total relative bkg uncert."; cout << string(58, '-') << "\n"; - double rate = yield_container.back().getVal(); - double err = yield_container.back().getError(); + double rate = yield_container->back().getVal(); + double err = yield_container->back().getError(); cout << boost::format("%-25s %-10.5f") % output_prefix.c_str() % (rate > 0. ? (err / rate) : 0.) << std::endl; } std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalSig" << std::endl; - load_shapes_and_yields(cmb_bin.cp().signals(), shape_container, yield_container, "TotalSig", false, no_sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin.cp().signals(), shape_container, "TotalSig", yield_container, false, no_sampling, samples, res, verbose); std::cout << ">> Doing " << output_prefix.c_str() << "," << "TotalProcs" << std::endl; - load_shapes_and_yields(cmb_bin, shape_container, yield_container, "TotalProcs", false, no_sampling, samples, res, verbose); + load_shapes_and_yields(cmb_bin, shape_container, "TotalProcs", yield_container, false, no_sampling, samples, res, verbose); if (datacard != "") { @@ -148,9 +156,11 @@ void do_complete_set( for (auto& iter : shape_container) { ch::WriteToTFile(&(iter.second), &outfile, output_prefix + "/" + iter.first); } - for (auto& yield: yield_container){ - helper.Form("%s/%s", output_prefix.c_str() , yield.GetName()); - ch::WriteToTFile(&(yield), &outfile, helper.Data()); + if(yield_container != nullptr){ + for (auto& yield: *yield_container){ + helper.Form("%s/%s", output_prefix.c_str() , yield.GetName()); + ch::WriteToTFile(&(yield), &outfile, helper.Data()); + } } } @@ -175,6 +185,7 @@ int main(int argc, char* argv[]) { bool skip_prefit = false; bool skip_proc_errs = false; bool total_shapes = false; + bool save_yields = false; bool verbose = false; std::vector reverse_bins_; std::vector bins_; @@ -218,6 +229,9 @@ int main(int argc, char* argv[]) { ("no-sampling", po::value(&no_sampling)->default_value(no_sampling)->implicit_value(true), "Do not use the cov. matrix sampling method for the post-fit uncertainty") + ("save-yields,y", + po::value(&save_yields)->default_value(save_yields)->implicit_value(true), + "Save process yields in addition to the shapes. This will generate RooRealVar objects with the yield as value and the (sampled) uncertainties as error") ("samples", po::value(&samples)->default_value(samples), "Number of samples to make in each evaluate call") @@ -393,26 +407,23 @@ int main(int argc, char* argv[]) { // We can always do the prefit version, // Loop through the bins writing the shapes to the output file map> pre_shapes; - map> pre_yields; TString helper; if (!skip_prefit) { if(total_shapes){ pre_shapes["total"] = map(); - pre_yields["total"] = vector(); - do_complete_set(cmb, outfile, "total_prefit", pre_shapes["total"], pre_yields["total"], merged_procs, skip_procs, cmb_card, datacard, skip_proc_errs, factors); + do_complete_set(cmb, outfile, "total_prefit", pre_shapes["total"], merged_procs, skip_procs, cmb_card, save_yields, datacard, skip_proc_errs, factors); } for (auto bin : bins) { pre_shapes[bin] = map(); - pre_yields[bin] = vector(); do_complete_set( cmb.cp().bin({bin}), outfile, bin+"_prefit", pre_shapes[bin], - pre_yields[bin], merged_procs, skip_procs, - cmb_card.cp().bin({bin}), + cmb_card.cp().bin({bin}), + save_yields, datacard, skip_proc_errs, factors, @@ -430,7 +441,6 @@ int main(int argc, char* argv[]) { // Calculate the post-fit fractional background uncertainty in each bin map> post_shapes; - map> post_yields; map post_yield_cov; map post_yield_cor; @@ -439,17 +449,16 @@ int main(int argc, char* argv[]) { if(total_shapes){ - post_yields["total"] = vector(); post_shapes["total"] = map(); do_complete_set( cmb, outfile, "total_postfit", post_shapes["total"], - post_yields["total"], merged_procs, skip_procs, cmb_card, + save_yields, datacard, skip_proc_errs, factors, @@ -462,7 +471,6 @@ int main(int argc, char* argv[]) { for (auto bin : bins) { - post_yields[bin] = vector(); post_shapes[bin] = map(); auto cmb_bin = cmb.cp().bin({bin}); do_complete_set( @@ -470,10 +478,10 @@ int main(int argc, char* argv[]) { outfile, bin+"_postfit", post_shapes[bin], - post_yields[bin], merged_procs, skip_procs, cmb_card.cp().bin({bin}), + save_yields, datacard, skip_proc_errs, factors, From f79ffa9af809b1319f2138fbe654880d1889ec79 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Wed, 8 Mar 2023 14:54:11 +0100 Subject: [PATCH 08/13] optimized verbose mode for GetShapeWithUncertainties, introduced new threshold parameter to report in verbose mode (not excessible right now) --- .../bin/PostFitShapesFromWorkspace.cpp | 2 +- CombineTools/interface/CombineHarvester.h | 4 +- CombineTools/src/CombineHarvester_Evaluate.cc | 50 +++++++++++-------- 3 files changed, 31 insertions(+), 25 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index e0802db8325..bc433985b72 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -49,7 +49,7 @@ void load_shapes_and_yields( } else { shape_container[proc] = no_sampling ? cmb.GetShapeWithUncertainty() - : cmb.GetShapeWithUncertainty(*res, samples); + : cmb.GetShapeWithUncertainty(*res, samples, verbose); } // fill the yield container if there is one if(yield_container != nullptr){ diff --git a/CombineTools/interface/CombineHarvester.h b/CombineTools/interface/CombineHarvester.h index c8118c562e0..97906efb44b 100644 --- a/CombineTools/interface/CombineHarvester.h +++ b/CombineTools/interface/CombineHarvester.h @@ -363,8 +363,8 @@ class CombineHarvester { * version of this method instead - this method will be removed in an * upcoming release */ - TH1F GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples, const bool& verbose = false); - TH1F GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples, const bool& verbose = false); + TH1F GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples, const bool& verbose = false, const float& fraction_threshold=0.5); + TH1F GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples, const bool& verbose = false, const float& fraction_threshold=0.5); TH1F GetObservedShape(); TH2F GetRateCovariance(RooFitResult const& fit, unsigned n_samples); diff --git a/CombineTools/src/CombineHarvester_Evaluate.cc b/CombineTools/src/CombineHarvester_Evaluate.cc index f260a58c4df..c3195d5294e 100644 --- a/CombineTools/src/CombineHarvester_Evaluate.cc +++ b/CombineTools/src/CombineHarvester_Evaluate.cc @@ -128,20 +128,23 @@ TH1F CombineHarvester::GetShapeWithUncertainty() { TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples, - const bool& verbose) { - return GetShapeWithUncertainty(*fit, n_samples, verbose); + const bool& verbose, + const float& fraction_threshold) { + return GetShapeWithUncertainty(*fit, n_samples, verbose, fraction_threshold); } TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples, - const bool& verbose) + const bool& verbose, + const float& fraction_threshold + ) { auto lookup = GenerateProcSystMap(); TH1F shape = GetShapeInternal(lookup); std::vector full_rand_shape_summary; for (int i = 1; i <= shape.GetNbinsX(); ++i) { shape.SetBinError(i, 0.0); - full_rand_shape_summary.push_back(""); + if (verbose) full_rand_shape_summary.push_back(""); } // Create a backup copy of the current parameter values auto backup = GetParameters(); @@ -162,25 +165,26 @@ TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, } TString tmp; - + std::string rand_shape_summary = ""; // Main loop through n_samples for (unsigned i = 0; i < n_samples; ++i) { // Randomise and update values fit.randomizePars(); - std::string rand_shape_summary = ""; for (int n = 0; n < n_pars; ++n) { if (p_vec[n] && !p_vec[n]->frozen()) { p_vec[n]->set_val(r_vec[n]->getVal()); - tmp.Form("setting parameter '%s' to %f\n", p_vec[n]->name().c_str(), r_vec[n]->getVal()); - // std::cout << tmp.Data(); - rand_shape_summary = tmp.Data(); - // std::cout << rand_shape_summary.c_str() << std::endl; - TH1F rand_shape = this->GetShapeInternal(lookup); - for (int j = 1; j <= rand_shape.GetNbinsX(); ++j) { - tmp.Form("\tBin Content %i: %f\n", j, rand_shape.GetBinContent(j)); + if(verbose){ + tmp.Form("setting parameter '%s' to %f\n", p_vec[n]->name().c_str(), r_vec[n]->getVal()); // std::cout << tmp.Data(); - full_rand_shape_summary.at(j-1) = rand_shape_summary; - full_rand_shape_summary.at(j-1) += tmp.Data(); + rand_shape_summary = tmp.Data(); + // std::cout << rand_shape_summary.c_str() << std::endl; + TH1F rand_shape = this->GetShapeInternal(lookup); + for (int j = 1; j <= rand_shape.GetNbinsX(); ++j) { + tmp.Form("\tBin Content %i: %f\n", j, rand_shape.GetBinContent(j)); + // std::cout << tmp.Data(); + full_rand_shape_summary.at(j-1) = rand_shape_summary; + full_rand_shape_summary.at(j-1) += tmp.Data(); + } } } // else if (!p_vec[n]) { @@ -192,13 +196,15 @@ TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, double err = std::fabs(rand_shape.GetBinContent(j) - shape.GetBinContent(j)); - if (std::fabs(err/shape.GetBinContent(j)) >= 0.5 and verbose){ - std::cout << "rel. error > 0.5 detected in bin " << j << " for toy " << i << std::endl; - std::cout << "nominal bin content (" << j << "):\t" << shape.GetBinContent(j) <= fraction_threshold){ + std::cout << "rel. error > " << fraction_threshold << " detected in bin " << j << " for toy " << i << std::endl; + std::cout << "nominal bin content (" << j << "):\t" << shape.GetBinContent(j) < Date: Tue, 28 Jul 2020 11:47:16 +0200 Subject: [PATCH 09/13] Add option for pattern-based merged bins in PostFitShapesFromWorkspace Conflicts: CombineTools/bin/PostFitShapesFromWorkspace.cpp --- .../bin/PostFitShapesFromWorkspace.cpp | 35 ++++++++++++++----- .../python/combine/EnhancedCombine.py | 1 - CombineTools/python/plotting.py | 6 ++++ CombineTools/src/CombineHarvester_Evaluate.cc | 4 ++- 4 files changed, 36 insertions(+), 10 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index bc433985b72..0b0f88ce2c2 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -194,6 +194,7 @@ int main(int argc, char* argv[]) { std::vector input_merge_procs_; std::map merged_procs; + std::vector merged_bins_; po::options_description help_config("Help"); help_config.add_options() @@ -261,7 +262,8 @@ int main(int argc, char* argv[]) { ("merge-procs,p", po::value>(&input_merge_procs_)->multitoken(), "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'") ("skip-procs,s", po::value>(&skip_procs_)->multitoken(), - "Skip these processes. Regex expression allowed. Format: 'expression'. Can be called multiple times"); + "Skip these processes. Regex expression allowed. Format: 'expression'. Can be called multiple times") + ("merged-bins", po::value>(&merged_bins_)->multitoken(), "List of [label]=[regex] for merged bins"); // if (!no_sampling && !postfit) { // throw logic_error( @@ -387,8 +389,7 @@ int main(int argc, char* argv[]) { vector bins; if (bins_.size() == 0) { - auto bin_set = cmb.cp().bin_set(); - std::copy(bin_set.begin(), bin_set.end(), std::back_inserter(bins)); + bins = ch::Set2Vec(cmb.cp().bin_set()); } else{ bins = bins_; @@ -401,6 +402,16 @@ int main(int argc, char* argv[]) { skip_procs.insert(process_set.begin(), process_set.end()); } + auto bin_patterns = bins; + for (unsigned i = 0; i < merged_bins_.size(); ++i) { + vector parts; + boost::split(parts, merged_bins_[i], boost::is_any_of("=")); + if (parts.size() == 2) { + bins.push_back(parts[0]); + bin_patterns.push_back(parts[1]); + } + } + TFile outfile(output.c_str(), "RECREATE"); TH1::AddDirectory(false); @@ -413,16 +424,21 @@ int main(int argc, char* argv[]) { pre_shapes["total"] = map(); do_complete_set(cmb, outfile, "total_prefit", pre_shapes["total"], merged_procs, skip_procs, cmb_card, save_yields, datacard, skip_proc_errs, factors); } - for (auto bin : bins) { + + + for (unsigned ib = 0; ib < bins.size(); ++ib) { + std::string bin = bins[ib]; + std::string bin_pattern = bin_patterns[ib]; pre_shapes[bin] = map(); + do_complete_set( - cmb.cp().bin({bin}), + cmb.cp().bin({bin_pattern}), outfile, bin+"_prefit", pre_shapes[bin], merged_procs, skip_procs, - cmb_card.cp().bin({bin}), + cmb_card.cp().bin({bin_pattern}), save_yields, datacard, skip_proc_errs, @@ -471,8 +487,11 @@ int main(int argc, char* argv[]) { for (auto bin : bins) { + for (unsigned ib = 0; ib < bins.size(); ++ib) { + std::string bin = bins[ib]; + std::string bin_pattern = bin_patterns[ib]; post_shapes[bin] = map(); - auto cmb_bin = cmb.cp().bin({bin}); + auto cmb_bin = cmb.cp().bin({bin_pattern}); do_complete_set( cmb_bin, outfile, @@ -480,7 +499,7 @@ int main(int argc, char* argv[]) { post_shapes[bin], merged_procs, skip_procs, - cmb_card.cp().bin({bin}), + cmb_card.cp().bin({bin_pattern}), save_yields, datacard, skip_proc_errs, diff --git a/CombineTools/python/combine/EnhancedCombine.py b/CombineTools/python/combine/EnhancedCombine.py index a69e701337b..1b17459653f 100755 --- a/CombineTools/python/combine/EnhancedCombine.py +++ b/CombineTools/python/combine/EnhancedCombine.py @@ -43,7 +43,6 @@ def attach_intercept_args(self, group): group.add_argument( '--setParameterRanges', help='Some other options will modify or add to the list of parameter ranges') - def attach_args(self, group): CombineToolBase.attach_args(self, group) group.add_argument( diff --git a/CombineTools/python/plotting.py b/CombineTools/python/plotting.py index ee734c2496a..d090c9d0334 100644 --- a/CombineTools/python/plotting.py +++ b/CombineTools/python/plotting.py @@ -383,6 +383,7 @@ def MultiRatioSplit(split_points, gaps_low, gaps_high): pad.Draw() pads.append(pad) pads.reverse() + pads[0].cd() return pads @@ -1191,6 +1192,9 @@ def FixTopRange(pad, fix_y, fraction): if ymin == 0.: print('Cannot adjust log-scale y-axis range if the minimum is zero!') return + if fix_y <= 0: + print 'Cannot adjust log-scale y-axis range if the maximum is zero!' + return maxval = (math.log10(fix_y) - fraction * math.log10(ymin)) / \ (1 - fraction) maxval = math.pow(10, maxval) @@ -1355,6 +1359,8 @@ def FixBoxPadding(pad, box, frac): # Convert this to a normalised frame value f_max_h = (a_max_h - ymin) / (ymax - ymin); if R.gPad.GetLogy() and f_max_h > 0.: + if a_max_h <= 0. or ymin <= 0. or ymax <= 0.: + return f_max_h = (math.log10(a_max_h) - math.log10(ymin)) / (math.log10(ymax) - math.log10(ymin)) if f_y1 - f_max_h < frac: diff --git a/CombineTools/src/CombineHarvester_Evaluate.cc b/CombineTools/src/CombineHarvester_Evaluate.cc index c3195d5294e..9f698ca853a 100644 --- a/CombineTools/src/CombineHarvester_Evaluate.cc +++ b/CombineTools/src/CombineHarvester_Evaluate.cc @@ -487,7 +487,9 @@ TH1F CombineHarvester::GetObservedShape() { tmp->SetBinErrorOption(TH1::kPoisson); proc_shape = *tmp; delete tmp; - proc_shape.Scale(1. / proc_shape.Integral()); + if (proc_shape.Integral() > 0.) { + proc_shape.Scale(1. / proc_shape.Integral()); + } } proc_shape.Scale(p_rate); if (!shape_init) { From 521207cbca5dfa5ce525cfa5249ccebcb3cea800 Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Mon, 6 Mar 2023 14:54:36 +0100 Subject: [PATCH 10/13] finalized cherry-pick --- CombineTools/bin/PostFitShapesFromWorkspace.cpp | 7 +++---- CombineTools/python/plotting.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/CombineTools/bin/PostFitShapesFromWorkspace.cpp b/CombineTools/bin/PostFitShapesFromWorkspace.cpp index 0b0f88ce2c2..75553f968c4 100644 --- a/CombineTools/bin/PostFitShapesFromWorkspace.cpp +++ b/CombineTools/bin/PostFitShapesFromWorkspace.cpp @@ -263,7 +263,7 @@ int main(int argc, char* argv[]) { "Merge these processes. Regex expression allowed. Format: NEWPROCESSNAME='expression'") ("skip-procs,s", po::value>(&skip_procs_)->multitoken(), "Skip these processes. Regex expression allowed. Format: 'expression'. Can be called multiple times") - ("merged-bins", po::value>(&merged_bins_)->multitoken(), "List of [label]=[regex] for merged bins"); + ("merged-channels", po::value>(&merged_bins_)->multitoken(), "List of [label]=[regex] for merged channels"); // if (!no_sampling && !postfit) { // throw logic_error( @@ -426,7 +426,7 @@ int main(int argc, char* argv[]) { } - for (unsigned ib = 0; ib < bins.size(); ++ib) { + for (unsigned ib = 0; ib < bins.size(); ++ib) { std::string bin = bins[ib]; std::string bin_pattern = bin_patterns[ib]; pre_shapes[bin] = map(); @@ -486,8 +486,7 @@ int main(int argc, char* argv[]) { } - for (auto bin : bins) { - for (unsigned ib = 0; ib < bins.size(); ++ib) { + for (unsigned ib = 0; ib < bins.size(); ++ib) { std::string bin = bins[ib]; std::string bin_pattern = bin_patterns[ib]; post_shapes[bin] = map(); diff --git a/CombineTools/python/plotting.py b/CombineTools/python/plotting.py index d090c9d0334..2afa00da6f2 100644 --- a/CombineTools/python/plotting.py +++ b/CombineTools/python/plotting.py @@ -1193,7 +1193,7 @@ def FixTopRange(pad, fix_y, fraction): print('Cannot adjust log-scale y-axis range if the minimum is zero!') return if fix_y <= 0: - print 'Cannot adjust log-scale y-axis range if the maximum is zero!' + print('Cannot adjust log-scale y-axis range if the maximum is zero!') return maxval = (math.log10(fix_y) - fraction * math.log10(ymin)) / \ (1 - fraction) From a742eb5e030b6ceecd23c480fc298ad0ddd613dd Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Wed, 28 Aug 2019 10:15:30 +0200 Subject: [PATCH 11/13] removed LimitCompare --- CombineTools/bin/LimitCompare.cpp | 149 ------------------------------ 1 file changed, 149 deletions(-) delete mode 100644 CombineTools/bin/LimitCompare.cpp diff --git a/CombineTools/bin/LimitCompare.cpp b/CombineTools/bin/LimitCompare.cpp deleted file mode 100644 index 9c6567dc202..00000000000 --- a/CombineTools/bin/LimitCompare.cpp +++ /dev/null @@ -1,149 +0,0 @@ -#include -#include -#include -#include -#include "TGraph.h" -#include "TCanvas.h" -#include "TStyle.h" -#include "TLegend.h" -#include "TTree.h" -#include "TFile.h" -#include "TAxis.h" -#include "TLatex.h" -#include "boost/lexical_cast.hpp" -#include "boost/algorithm/string.hpp" -#include "boost/format.hpp" -#include "boost/program_options.hpp" -#include "CombineHarvester/CombineTools/interface/Plotting.h" -#include "CombineHarvester/CombineTools/interface/Plotting_Style.h" -// #include "Core/interface/TextElement.h" -// #include "Utilities/interface/SimpleParamParser.h" -// #include "Utilities/interface/FnRootTools.h" - -namespace po = boost::program_options; - -using namespace std; - -void SetTGraphStyle(TGraph * gr, int color) { - gr->SetLineStyle(11.); - gr->SetLineWidth( 3.); - gr->SetLineColor(color); - gr->SetMarkerStyle(20); - gr->SetMarkerSize(1.3); - gr->SetMarkerColor(color); -} - -// "mh:limit","quantileExpected==0.5" - -TGraph ExtractExpected(TTree *limit) { - int n = limit->Draw("mh:limit", "quantileExpected==0.5", "goff"); - TGraph gr(n, limit->GetV1(), limit->GetV2()); - gr.Sort(); - return gr; -} - -struct Infos { - std::string file; - std::string label; - int color; -}; - -int main(int /*argc*/, char* argv[]) { - - std::string ch = argv[1]; - vector input; - input.push_back({"parametric/higgsCombineExpected"+ch+".Asymptotic.All.root", "Parametric", 2}); - input.push_back({"parametric-binned/higgsCombineExpected"+ch+".Asymptotic.All.root", "Binned", 4}); - // input.push_back({"parametric-asimov/higgsCombineExpected"+ch+".Asymptotic.All.root", "Parametric (asimov)", 5}); - // input.push_back({"parametric-binned-asimov/higgsCombineExpected"+ch+".Asymptotic.All.root", "Binned (asimov)", 6}); - - - vector graphs; - - for (auto const& in : input) { - TFile f(in.file.c_str()); - TTree *tree = (TTree*)f.Get("limit"); - graphs.push_back(ExtractExpected(tree)); - } - - for (unsigned i = 0; i < input.size(); ++i) { - SetTGraphStyle(&graphs[i], input[i].color); - } - - - ModTDRStyle(); - TCanvas canv("canv","canv",600,600); - canv.SetLogy(true); - gStyle->SetNdivisions(505); - canv.SetGridx(1); - canv.SetGridy(1); - TGraph & exp1 = graphs[0]; - exp1.GetXaxis()->SetLimits(exp1.GetX()[0]-.1, exp1.GetX()[exp1.GetN()-1]+.1); - - double min = 9999.; - double max = 0.; - for (int i = 0; i < exp1.GetN(); ++i) { - double x, y; - exp1.GetPoint(i, x, y); - if (y < min) min = y; - if (y > max) max = y; - } - exp1.SetMaximum(max*2.); - exp1.SetMinimum(min*0.5); - exp1.GetYaxis()->SetTitleSize(0.055); - exp1.GetYaxis()->SetTitleOffset(1.200); - exp1.GetYaxis()->SetLabelOffset(0.014); - exp1.GetYaxis()->SetLabelSize(0.040); - exp1.GetYaxis()->SetLabelFont(42); - exp1.GetYaxis()->SetTitleFont(42); - exp1.GetXaxis()->SetTitleSize(0.055); - exp1.GetXaxis()->SetTitleOffset(1.100); - exp1.GetXaxis()->SetLabelOffset(0.014); - exp1.GetXaxis()->SetLabelSize(0.040); - exp1.GetXaxis()->SetLabelFont(42); - exp1.GetXaxis()->SetTitleFont(42); - exp1.Draw("APL"); - - exp1.GetXaxis()->SetTitle("m_{#Phi} [GeV]"); - exp1.GetYaxis()->SetTitle("95% CL Limit on #sigma #times B [pb]"); - for (unsigned i = 1; i < input.size(); ++i) { - graphs[i].Draw("PLsame"); - } - // SetTGraphStyle(&obs1, 6); - // SetTGraphStyle(&exp2, 7); - // SetTGraphStyle(&obs2, 8); - - // obs1.Draw("PLsame"); - // // exp2.Draw("PLsame"); - // // obs2.Draw("PLsame"); - - TLegend *legend = new TLegend(0.62, 0.62, 0.93, 0.78, "", "brNDC"); - for (unsigned i = 0; i < input.size(); ++i) { - legend->AddEntry(&graphs[i], input[i].label.c_str(), "lP"); - } - legend->SetBorderSize(1); - legend->SetTextFont(42); - legend->SetLineColor(0); - legend->SetLineStyle(1); - legend->SetLineWidth(1); - legend->SetFillColor(0); - legend->SetFillStyle(1001); - legend->Draw(); - - TLatex *title_latex = new TLatex(); - title_latex->SetNDC(); - title_latex->SetTextSize(0.035); - title_latex->SetTextFont(62); - title_latex->SetTextAlign(31); - double height = 0.94; - title_latex->DrawLatex(0.95,height,"19.7 fb^{-1} (8 TeV)"); - title_latex->SetTextAlign(11); - title_latex->DrawLatex(0.17,height,"H#rightarrow#tau#tau"); - title_latex->SetTextSize(0.06); - title_latex->DrawLatex(0.78, 0.84, "#mu#tau_{h}"); - - canv.SaveAs((ch+"_comparison.pdf").c_str()); - - return 0; -} - From 69ba4f3e8070af704d66a88475b3049e850edb07 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Thu, 27 Apr 2023 15:11:04 +0200 Subject: [PATCH 12/13] removed LimitCompare from Build xml --- CombineTools/bin/BuildFile.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/CombineTools/bin/BuildFile.xml b/CombineTools/bin/BuildFile.xml index c967a15b5e0..95b972e0e58 100644 --- a/CombineTools/bin/BuildFile.xml +++ b/CombineTools/bin/BuildFile.xml @@ -3,7 +3,6 @@ - From 9397ef7ce9ceda93107d1576798de5bf469dc098 Mon Sep 17 00:00:00 2001 From: Philip Daniel Keicher Date: Thu, 1 Jun 2023 15:35:18 +0200 Subject: [PATCH 13/13] resolve minor merge conflict issues with main --- CombineTools/src/CombineHarvester_Evaluate.cc | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/CombineTools/src/CombineHarvester_Evaluate.cc b/CombineTools/src/CombineHarvester_Evaluate.cc index 769aea237b5..08e217203cf 100644 --- a/CombineTools/src/CombineHarvester_Evaluate.cc +++ b/CombineTools/src/CombineHarvester_Evaluate.cc @@ -209,13 +209,14 @@ TH1F CombineHarvester::GetShapeWithUncertainty(RooFitResult const& fit, for (auto x :sample_yields_perbin[i-1]){ double err = std::fabs(x - yield_mean); if(verbose){ - if (std::fabs(err/shape.GetBinContent(j)) >= fraction_threshold){ - std::cout << "rel. error > " << fraction_threshold << " detected in bin " << j << " for toy " << i << std::endl; - std::cout << "nominal bin content (" << j << "):\t" << shape.GetBinContent(j) <= fraction_threshold){ + std::cout << "rel. error > " << fraction_threshold << " detected in bin " << i << std::endl; + std::cout << "nominal bin content (" << i << "):\t" << orig_content <