diff --git a/CombineHarvester/CombineTools/interface/BinByBin.h b/CombineHarvester/CombineTools/interface/BinByBin.h new file mode 100644 index 00000000..efda0e14 --- /dev/null +++ b/CombineHarvester/CombineTools/interface/BinByBin.h @@ -0,0 +1,51 @@ +#ifndef CombineTools_BinByBin_h +#define CombineTools_BinByBin_h +#include +#include "CombineTools/interface/CombineHarvester.h" + +namespace ch { + +class BinByBinFactory { + public: + BinByBinFactory(); + + void MergeBinErrors(CombineHarvester &cb); + + void AddBinByBin(CombineHarvester &src, CombineHarvester &dest); + + void MergeAndAdd(CombineHarvester &src, CombineHarvester &dest); + + inline BinByBinFactory& SetVerbosity(unsigned verbosity) { + v_ = verbosity; + return *this; + } + + inline BinByBinFactory& SetAddThreshold(double val) { + bbb_threshold_ = val; + return *this; + } + inline BinByBinFactory& SetMergeThreshold(double val) { + merge_threshold_ = val; + return *this; + } + + inline BinByBinFactory& SetPattern(std::string const& pattern) { + pattern_ = pattern; + return *this; + } + + inline BinByBinFactory& SetFixNorm(bool fix) { + fix_norm_ = fix; + return *this; + } + + private: + std::string pattern_; + unsigned v_; + double bbb_threshold_; + double merge_threshold_; + bool fix_norm_; +}; +} + +#endif diff --git a/CombineHarvester/CombineTools/interface/CombineHarvester.h b/CombineHarvester/CombineTools/interface/CombineHarvester.h index 9c952fbf..8df00d88 100644 --- a/CombineHarvester/CombineTools/interface/CombineHarvester.h +++ b/CombineHarvester/CombineTools/interface/CombineHarvester.h @@ -363,16 +363,32 @@ class CombineHarvester { void InsertObservation(ch::Observation const& obs); void InsertProcess(ch::Process const& proc); void InsertSystematic(ch::Systematic const& sys); + void CreateParameterIfEmpty(std::string const& name); /** * Create bin-by-bin uncertainties * - * \deprecated Please use the AddBinByBin(double, bool, CombineHarvester &) - * interface instead + * \deprecated Please use the standalone ch::BinByBinFactory tool, defined + * in CombineTools/interface/BinByBin.h */ void AddBinByBin(double threshold, bool fixed_norm, CombineHarvester* other); + + + /** + * Create bin-by-bin uncertainties + * + * \deprecated Please use the standalone ch::BinByBinFactory tool, defined + * in CombineTools/interface/BinByBin.h + */ void AddBinByBin(double threshold, bool fixed_norm, CombineHarvester & other); + + /** + * Merge bin errors within a bin + * + * \deprecated Please use the standalone ch::BinByBinFactory tool, defined + * in CombineTools/interface/BinByBin.h + */ void MergeBinErrors(double bbb_threshold, double merge_threshold); /**@}*/ @@ -472,8 +488,6 @@ void FillHistMappings(std::vector & mappings); TH1 const* high, bool linear); void ShapeDiff(double x, TH1F* target, RooDataHist const* nom, RooDataHist const* low, RooDataHist const* high); - - void CreateParameterIfEmpty(CombineHarvester *cmb, std::string const& name); }; diff --git a/CombineHarvester/CombineTools/src/BinByBin.cc b/CombineHarvester/CombineTools/src/BinByBin.cc new file mode 100644 index 00000000..b50f4bb3 --- /dev/null +++ b/CombineHarvester/CombineTools/src/BinByBin.cc @@ -0,0 +1,157 @@ +#include "CombineTools/interface/BinByBin.h" +#include +#include +#include +#include "boost/format.hpp" +#include "boost/lexical_cast.hpp" + +namespace ch { + +BinByBinFactory::BinByBinFactory() + : pattern_("CMS_$ANALYSIS_$CHANNEL_$BIN_$ERA_$PROCESS_bin_$#"), + v_(0), + bbb_threshold_(0.), + merge_threshold_(0.), + fix_norm_(true) {} + + +void BinByBinFactory::MergeBinErrors(CombineHarvester &cb) { + // Reduce merge_threshold very slightly to avoid numerical issues + // E.g. two backgrounds each with bin error 1.0. merge_threshold of + // 0.5 should not result in merging - but can do depending on + // machine and compiler + double merge_threshold = merge_threshold_ - 1E-9 * merge_threshold_; + auto bins = cb.bin_set(); + for (auto const& bin : bins) { + unsigned bbb_added = 0; + unsigned bbb_removed = 0; + CombineHarvester tmp = std::move(cb.cp().bin({bin}).histograms()); + std::vector procs; + tmp.ForEachProc([&](Process *p) { + procs.push_back(p); + }); + if (procs.size() == 0) continue; + + std::vector> h_copies(procs.size()); + for (unsigned i = 0; i < h_copies.size(); ++i) { + h_copies[i] = procs[i]->ClonedScaledShape(); + } + + for (int i = 1; i <= h_copies[0]->GetNbinsX(); ++i) { + double tot_bbb_added = 0.0; + std::vector> result; + for (unsigned j = 0; j < h_copies.size(); ++j) { + double val = h_copies[j]->GetBinContent(i); + double err = h_copies[j]->GetBinError(i); + if (val == 0.0 && err == 0.0) continue; + if (val == 0 || (err/val) > bbb_threshold_) { + bbb_added += 1; + tot_bbb_added += (err * err); + result.push_back(std::make_pair(err*err, h_copies[j].get())); + } + } + if (tot_bbb_added == 0.0) continue; + std::sort(result.begin(), result.end()); + double removed = 0.0; + for (unsigned r = 0; r < result.size(); ++r) { + if ((result[r].first + removed) < (merge_threshold * tot_bbb_added) && + r < (result.size() - 1)) { + bbb_removed += 1; + removed += result[r].first; + result[r].second->SetBinError(i, 0.0); + } + } + double expand = std::sqrt(1. / (1. - (removed / tot_bbb_added))); + for (unsigned r = 0; r < result.size(); ++r) { + result[r] + .second->SetBinError(i, result[r].second->GetBinError(i) * expand); + } + } + for (unsigned i = 0; i < h_copies.size(); ++i) { + procs[i]->set_shape(std::move(h_copies[i]), false); + } + if (v_ > 0) { + std::cout << "BIN: " << bin << std::endl; + std::cout << "Total bbb added: " << bbb_added << "\n"; + std::cout << "Total bbb removed: " << bbb_removed << "\n"; + std::cout << "Total bbb =======>: " << bbb_added-bbb_removed << "\n"; + } + } +} + +void BinByBinFactory::AddBinByBin(CombineHarvester &src, CombineHarvester &dest) { + unsigned bbb_added = 0; + std::vector procs; + src.ForEachProc([&](Process *p) { + procs.push_back(p); + }); + for (unsigned i = 0; i < procs.size(); ++i) { + if (!procs[i]->shape()) continue; + TH1 const* h = procs[i]->shape(); + unsigned n_pop_bins = 0; + for (int j = 1; j <= h->GetNbinsX(); ++j) { + if (h->GetBinContent(j) > 0.0) ++n_pop_bins; + } + if (n_pop_bins <= 1 && fix_norm_) { + if (v_ >= 1) { + std::cout << "Requested fixed_norm but template has <= 1 populated " + "bins, skipping\n"; + std::cout << Process::PrintHeader << *(procs[i]) << "\n"; + } + continue; + } + for (int j = 1; j <= h->GetNbinsX(); ++j) { + bool do_bbb = false; + double val = h->GetBinContent(j); + double err = h->GetBinError(j); + if (val == 0. && err > 0.) do_bbb = true; + if (val > 0. && (err / val) > bbb_threshold_) do_bbb = true; + // if (h->GetBinContent(j) <= 0.0) { + // if (h->GetBinError(j) > 0.0) { + // std::cout << *(procs_[i]) << "\n"; + // std::cout << "Bin with content <= 0 and error > 0 found, skipping\n"; + // } + // continue; + // } + if (do_bbb) { + ++bbb_added; + ch::Systematic sys; + ch::SetProperties(&sys, procs[i]); + sys.set_type("shape"); + std::string name = pattern_; + boost::replace_all(name, "$ANALYSIS", sys.analysis()); + boost::replace_all(name, "$CHANNEL", sys.channel()); + boost::replace_all(name, "$BIN", sys.bin()); + boost::replace_all(name, "$BINID", boost::lexical_cast(sys.bin_id())); + boost::replace_all(name, "$ERA", sys.era()); + boost::replace_all(name, "$PROCESS", sys.process()); + boost::replace_all(name, "$MASS", sys.mass()); + boost::replace_all(name, "$#", boost::lexical_cast(j)); + sys.set_name(name); + sys.set_asymm(true); + std::unique_ptr h_d(static_cast(h->Clone())); + std::unique_ptr h_u(static_cast(h->Clone())); + h_d->SetBinContent(j, val - err); + if (h_d->GetBinContent(j) < 0.) h_d->SetBinContent(j, 0.); + h_u->SetBinContent(j, val + err); + if (fix_norm_) { + sys.set_value_d(1.0); + sys.set_value_u(1.0); + } else { + sys.set_value_d(h_d->Integral()/h->Integral()); + sys.set_value_u(h_u->Integral()/h->Integral()); + } + sys.set_shapes(std::move(h_u), std::move(h_d), nullptr); + dest.CreateParameterIfEmpty(sys.name()); + dest.InsertSystematic(sys); + } + } + } + // std::cout << "bbb added: " << bbb_added << std::endl; +} + +void BinByBinFactory::MergeAndAdd(CombineHarvester &src, CombineHarvester &dest) { + MergeBinErrors(src); + AddBinByBin(src, dest); +} +} diff --git a/CombineHarvester/CombineTools/src/CombineHarvester_Creation.cc b/CombineHarvester/CombineTools/src/CombineHarvester_Creation.cc index 722280f7..15e6a47f 100644 --- a/CombineHarvester/CombineTools/src/CombineHarvester_Creation.cc +++ b/CombineHarvester/CombineTools/src/CombineHarvester_Creation.cc @@ -14,6 +14,7 @@ #include "CombineTools/interface/Parameter.h" #include "CombineTools/interface/Utilities.h" #include "CombineTools/interface/Logging.h" +#include "CombineTools/interface/BinByBin.h" namespace ch { void CombineHarvester::AddObservations( @@ -97,7 +98,7 @@ void CombineHarvester::AddSystFromProc(Process const& proc, sys->set_value_d(1.0); sys->set_scale(val_u); } - CreateParameterIfEmpty(this, sys->name()); + CreateParameterIfEmpty(sys->name()); if (sys->type() == "lnU") { params_.at(sys->name())->set_err_d(0.); params_.at(sys->name())->set_err_u(0.); @@ -178,138 +179,28 @@ void CombineHarvester::AddBinByBin(double threshold, bool fixed_norm, void CombineHarvester::AddBinByBin(double threshold, bool fixed_norm, CombineHarvester *other) { - unsigned bbb_added = 0; - for (unsigned i = 0; i < procs_.size(); ++i) { - if (!procs_[i]->shape()) continue; - TH1 const* h = procs_[i]->shape(); - unsigned n_pop_bins = 0; - for (int j = 1; j <= h->GetNbinsX(); ++j) { - if (h->GetBinContent(j) > 0.0) ++n_pop_bins; - } - if (n_pop_bins <= 1 && fixed_norm) { - if (verbosity_ >= 1) { - std::cout << "Requested fixed_norm but template has <= 1 populated " - "bins, skipping\n"; - std::cout << Process::PrintHeader << *(procs_[i]) << "\n"; - } - continue; - } - for (int j = 1; j <= h->GetNbinsX(); ++j) { - bool do_bbb = false; - double val = h->GetBinContent(j); - double err = h->GetBinError(j); - if (val == 0. && err > 0.) do_bbb = true; - if (val > 0. && (err / val) > threshold) do_bbb = true; - // if (h->GetBinContent(j) <= 0.0) { - // if (h->GetBinError(j) > 0.0) { - // std::cout << *(procs_[i]) << "\n"; - // std::cout << "Bin with content <= 0 and error > 0 found, skipping\n"; - // } - // continue; - // } - if (do_bbb) { - ++bbb_added; - auto sys = std::make_shared(); - ch::SetProperties(sys.get(), procs_[i].get()); - sys->set_type("shape"); - // sys->set_name("CMS_" + sys->bin() + "_" + sys->process() + "_bin_" + - // boost::lexical_cast(j)); - sys->set_name("CMS_" + sys->analysis() + "_" + sys->channel() + "_" + - sys->bin() + "_" + sys->era() + "_" + sys->process() + - "_bin_" + boost::lexical_cast(j)); - sys->set_asymm(true); - std::unique_ptr h_d(static_cast(h->Clone())); - std::unique_ptr h_u(static_cast(h->Clone())); - h_d->SetBinContent(j, val - err); - if (h_d->GetBinContent(j) < 0.) h_d->SetBinContent(j, 0.); - h_u->SetBinContent(j, val + err); - if (fixed_norm) { - sys->set_value_d(1.0); - sys->set_value_u(1.0); - } else { - sys->set_value_d(h_d->Integral()/h->Integral()); - sys->set_value_u(h_u->Integral()/h->Integral()); - } - sys->set_shapes(std::move(h_u), std::move(h_d), nullptr); - CombineHarvester::CreateParameterIfEmpty(other ? other : this, - sys->name()); - if (other) { - other->systs_.push_back(sys); - } else { - systs_.push_back(sys); - } - } - } - } - // std::cout << "bbb added: " << bbb_added << std::endl; + auto bbb_factory = ch::BinByBinFactory() + .SetAddThreshold(threshold) + .SetFixNorm(fixed_norm) + .SetVerbosity(verbosity_); + bbb_factory.AddBinByBin(*this, *other); } -void CombineHarvester::CreateParameterIfEmpty(CombineHarvester *cmb, - std::string const &name) { +void CombineHarvester::CreateParameterIfEmpty(std::string const &name) { if (!params_.count(name)) { auto param = std::make_shared(Parameter()); param->set_name(name); - (*cmb).params_.insert({name, param}); + params_.insert({name, param}); } } void CombineHarvester::MergeBinErrors(double bbb_threshold, double merge_threshold) { - // Reduce merge_threshold very slightly to avoid numerical issues - // E.g. two backgrounds each with bin error 1.0. merge_threshold of - // 0.5 should not result in merging - but can do depending on - // machine and compiler - merge_threshold -= 1E-9 * merge_threshold; - auto bins = this->bin_set(); - for (auto const& bin : bins) { - // unsigned bbb_added = 0; - // unsigned bbb_removed = 0; - CombineHarvester tmp = std::move(this->cp().bin({bin}).histograms()); - if (tmp.procs_.size() == 0) continue; - - std::vector> h_copies(tmp.procs_.size()); - for (unsigned i = 0; i < h_copies.size(); ++i) { - h_copies[i] = tmp.procs_[i]->ClonedScaledShape(); - } - - for (int i = 1; i <= h_copies[0]->GetNbinsX(); ++i) { - double tot_bbb_added = 0.0; - std::vector> result; - for (unsigned j = 0; j < h_copies.size(); ++j) { - double val = h_copies[j]->GetBinContent(i); - double err = h_copies[j]->GetBinError(i); - if (val == 0.0 && err == 0.0) continue; - if (val == 0 || (err/val) > bbb_threshold) { - // bbb_added += 1; - tot_bbb_added += (err * err); - result.push_back(std::make_pair(err*err, h_copies[j].get())); - } - } - if (tot_bbb_added == 0.0) continue; - std::sort(result.begin(), result.end()); - double removed = 0.0; - for (unsigned r = 0; r < result.size(); ++r) { - if ((result[r].first + removed) < (merge_threshold * tot_bbb_added) && - r < (result.size() - 1)) { - // bbb_removed += 1; - removed += result[r].first; - result[r].second->SetBinError(i, 0.0); - } - } - double expand = std::sqrt(1. / (1. - (removed / tot_bbb_added))); - for (unsigned r = 0; r < result.size(); ++r) { - result[r] - .second->SetBinError(i, result[r].second->GetBinError(i) * expand); - } - } - for (unsigned i = 0; i < h_copies.size(); ++i) { - tmp.procs_[i]->set_shape(std::move(h_copies[i]), false); - } - // std::cout << "BIN: " << bin << std::endl; - // std::cout << "Total bbb added: " << bbb_added << "\n"; - // std::cout << "Total bbb removed: " << bbb_removed << "\n"; - // std::cout << "Total bbb =======>: " << bbb_added-bbb_removed << "\n"; - } + auto bbb_factory = ch::BinByBinFactory() + .SetAddThreshold(bbb_threshold) + .SetMergeThreshold(merge_threshold) + .SetVerbosity(verbosity_); + bbb_factory.MergeBinErrors(*this); } void CombineHarvester::InsertObservation(ch::Observation const& obs) { diff --git a/CombineHarvester/CombineTools/src/CombineHarvester_Datacards.cc b/CombineHarvester/CombineTools/src/CombineHarvester_Datacards.cc index d4a554c1..27ba9303 100644 --- a/CombineHarvester/CombineTools/src/CombineHarvester_Datacards.cc +++ b/CombineHarvester/CombineTools/src/CombineHarvester_Datacards.cc @@ -318,7 +318,7 @@ int CombineHarvester::ParseDatacard(std::string const& filename, if (sys->type() == "shape" || sys->type() == "shapeN2") sys->set_asymm(true); - CombineHarvester::CreateParameterIfEmpty(this, sys->name()); + CombineHarvester::CreateParameterIfEmpty(sys->name()); if (sys->type() == "lnU") { params_.at(sys->name())->set_err_d(0.); params_.at(sys->name())->set_err_u(0.); diff --git a/CombineHarvester/CombineTools/src/CombineHarvester_Python.cc b/CombineHarvester/CombineTools/src/CombineHarvester_Python.cc index 77e76c86..cef2ca30 100644 --- a/CombineHarvester/CombineTools/src/CombineHarvester_Python.cc +++ b/CombineHarvester/CombineTools/src/CombineHarvester_Python.cc @@ -2,6 +2,7 @@ #include "CombineTools/interface/CombineHarvester.h" #include "CombineTools/interface/Observation.h" #include "CombineTools/interface/CardWriter.h" +#include "CombineTools/interface/BinByBin.h" #include "CombineTools/interface/CopyTools.h" #include "boost/python.hpp" #include "TFile.h" @@ -14,6 +15,7 @@ using ch::Observation; using ch::Process; using ch::Systematic; using ch::CardWriter; +using ch::BinByBinFactory; void FilterAllPy(ch::CombineHarvester & cb, boost::python::object func) { auto lambda = [func](ch::Object *obj) -> bool { @@ -338,4 +340,20 @@ BOOST_PYTHON_MODULE(libCHCombineTools) py::def("CloneProcs", CloneProcsPy); py::def("CloneSysts", CloneSystsPy); py::def("CloneProcsAndSysts", CloneProcsAndSystsPy); + + py::class_("BinByBinFactory") + .def("MergeBinErrors", &BinByBinFactory::MergeBinErrors) + .def("AddBinByBin", &BinByBinFactory::AddBinByBin) + .def("MergeAndAdd", &BinByBinFactory::MergeAndAdd) + .def("SetVerbosity", &BinByBinFactory::SetVerbosity, + py::return_internal_reference<>()) + .def("SetAddThreshold", &BinByBinFactory::SetAddThreshold, + py::return_internal_reference<>()) + .def("SetMergeThreshold", &BinByBinFactory::SetMergeThreshold, + py::return_internal_reference<>()) + .def("SetPattern", &BinByBinFactory::SetPattern, + py::return_internal_reference<>()) + .def("SetFixNorm", &BinByBinFactory::SetFixNorm, + py::return_internal_reference<>()) + ; } \ No newline at end of file diff --git a/CombineHarvester/CombineTools/test/Example2.cpp b/CombineHarvester/CombineTools/test/Example2.cpp index 999b243c..5bc68726 100644 --- a/CombineHarvester/CombineTools/test/Example2.cpp +++ b/CombineHarvester/CombineTools/test/Example2.cpp @@ -98,15 +98,15 @@ int main() { aux_shapes + "Imperial/htt_mt.inputs-sm-8TeV-hcg.root", "$BIN/$PROCESS$MASS", "$BIN/$PROCESS$MASS_$SYSTEMATIC"); + //! [part7] + + //! [part8] + cb.cp().backgrounds().AddBinByBin(0.1, true, &cb); // This function modifies every entry to have a standardised bin name of // the form: {analysis}_{channel}_{bin_id}_{era} // which is commonly used in the htt analyses ch::SetStandardBinNames(cb); - //! [part7] - - //! [part8] - cb.cp().backgrounds().AddBinByBin(0.1, true, &cb); //! [part8] //! [part9] diff --git a/CombineHarvester/CombineTools/test/SMLegacyExample.cpp b/CombineHarvester/CombineTools/test/SMLegacyExample.cpp index d9c396e3..c4f86a92 100644 --- a/CombineHarvester/CombineTools/test/SMLegacyExample.cpp +++ b/CombineHarvester/CombineTools/test/SMLegacyExample.cpp @@ -11,6 +11,7 @@ #include "CombineTools/interface/HttSystematics.h" #include "CombineTools/interface/CardWriter.h" #include "CombineTools/interface/CopyTools.h" +#include "CombineTools/interface/BinByBin.h" using namespace std; @@ -158,95 +159,44 @@ int main() { } } - cout << ">> Merging bin errors...\n"; - ch::CombineHarvester cb_et = move(cb.cp().channel({"et"})); - for (string era : {"7TeV", "8TeV"}) { - cb_et.cp().era({era}).bin_id({1, 2}).process({"ZL", "ZJ", "QCD", "W"}) - .MergeBinErrors(0.1, 0.5); - cb_et.cp().era({era}).bin_id({3, 5}).process({"W"}) - .MergeBinErrors(0.1, 0.5); - } - cb_et.cp().era({"7TeV"}).bin_id({6}).process({"ZL", "ZJ", "W", "ZTT"}) - .MergeBinErrors(0.1, 0.5); - cb_et.cp().era({"8TeV"}).bin_id({7}).process({"ZL", "ZJ", "W", "ZTT"}) - .MergeBinErrors(0.1, 0.5); - cb_et.cp().era({"8TeV"}).bin_id({6}).process({"ZL", "ZJ", "W"}) - .MergeBinErrors(0.1, 0.5); - - ch::CombineHarvester cb_mt = move(cb.cp().channel({"mt"})); - for (string era : {"7TeV", "8TeV"}) { - cb_mt.cp().era({era}).bin_id({1, 2, 3, 4}).process({"W", "QCD"}) - .MergeBinErrors(0.1, 0.5); - } - cb_mt.cp().era({"7TeV"}).bin_id({5}).process({"W"}) - .MergeBinErrors(0.1, 0.5); - cb_mt.cp().era({"7TeV"}).bin_id({6}).process({"W", "ZTT"}) - .MergeBinErrors(0.1, 0.5); - cb_mt.cp().era({"8TeV"}).bin_id({5, 6}).process({"W"}) - .MergeBinErrors(0.1, 0.5); - cb_mt.cp().era({"8TeV"}).bin_id({7}).process({"W", "ZTT"}) - .MergeBinErrors(0.1, 0.5); - - ch::CombineHarvester cb_em = move(cb.cp().channel({"em"})); - for (string era : {"7TeV", "8TeV"}) { - cb_em.cp().era({era}).bin_id({1, 3}).process({"Fakes"}) - .MergeBinErrors(0.1, 0.5); - } - cb_em.cp().era({"7TeV"}).bin_id({4}).process({"Fakes", "EWK", "Ztt"}) - .MergeBinErrors(0.1, 0.5); - cb_em.cp().era({"8TeV"}).bin_id({5}).process({"Fakes", "EWK", "Ztt"}) - .MergeBinErrors(0.1, 0.5); - cb_em.cp().era({"8TeV"}).bin_id({4}).process({"Fakes", "EWK"}) - .MergeBinErrors(0.1, 0.5); - - ch::CombineHarvester cb_ee_mm = move(cb.cp().channel({"ee", "mm"})); - for (string era : {"7TeV", "8TeV"}) { - cb_ee_mm.cp().era({era}).bin_id({1, 3, 4}) - .process({"ZTT", "ZEE", "ZMM", "TTJ"}) - .MergeBinErrors(0.0, 0.5); - } - - ch::CombineHarvester cb_tt = move(cb.cp().channel({"tt"})); - cb_tt.cp().bin_id({0, 1, 2}).era({"8TeV"}).process({"ZTT", "QCD"}) - .MergeBinErrors(0.1, 0.5); - - cout << ">> Generating bbb uncertainties...\n"; - cb_mt.cp().bin_id({0, 1, 2, 3, 4}).process({"W", "QCD"}) - .AddBinByBin(0.1, true, &cb); - cb_mt.cp().era({"7TeV"}).bin_id({5}).process({"W"}) - .AddBinByBin(0.1, true, &cb); - cb_mt.cp().era({"7TeV"}).bin_id({6}).process({"W", "ZTT"}) - .AddBinByBin(0.1, true, &cb); - cb_mt.cp().era({"8TeV"}).bin_id({5, 6}).process({"W"}) - .AddBinByBin(0.1, true, &cb); - cb_mt.cp().era({"8TeV"}).bin_id({7}).process({"W", "ZTT"}) - .AddBinByBin(0.1, true, &cb); - - cb_et.cp().bin_id({1, 2}).process({"ZL", "ZJ", "QCD", "W"}) - .AddBinByBin(0.1, true, &cb); - cb_et.cp().bin_id({3, 5}).process({"W"}) - .AddBinByBin(0.1, true, &cb); - cb_et.cp().era({"7TeV"}).bin_id({6}).process({"ZL", "ZJ", "W", "ZTT"}) - .AddBinByBin(0.1, true, &cb); - cb_et.cp().era({"8TeV"}).bin_id({7}).process({"ZL", "ZJ", "W", "ZTT"}) - .AddBinByBin(0.1, true, &cb); - cb_et.cp().era({"8TeV"}).bin_id({6}).process({"ZL", "ZJ", "W"}) - .AddBinByBin(0.1, true, &cb); - - cb_em.cp().bin_id({1, 3}).process({"Fakes"}) - .AddBinByBin(0.1, true, &cb); - cb_em.cp().era({"7TeV"}).bin_id({4}).process({"Fakes", "EWK", "Ztt"}) - .AddBinByBin(0.1, true, &cb); - cb_em.cp().era({"8TeV"}).bin_id({5}).process({"Fakes", "EWK", "Ztt"}) - .AddBinByBin(0.1, true, &cb); - cb_em.cp().era({"8TeV"}).bin_id({4}).process({"Fakes", "EWK"}) - .AddBinByBin(0.1, true, &cb); - - cb_ee_mm.cp().bin_id({1, 3, 4}).process({"ZTT", "ZEE", "ZMM", "TTJ"}) - .AddBinByBin(0.0, true, &cb); - - cb_tt.cp().bin_id({0, 1, 2}).era({"8TeV"}).process({"QCD", "ZTT"}) - .AddBinByBin(0.1, true, &cb); + cout << ">> Merging bin errors and generating bbb uncertainties...\n"; + + auto bbb = ch::BinByBinFactory() + .SetAddThreshold(0.1) + .SetMergeThreshold(0.5) + .SetFixNorm(true); + + ch::CombineHarvester cb_et = cb.cp().channel({"et"}); + bbb.MergeAndAdd(cb_et.cp().era({"7TeV"}).bin_id({1, 2}).process({"ZL", "ZJ", "QCD", "W"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"7TeV"}).bin_id({3, 5}).process({"W"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"8TeV"}).bin_id({1, 2}).process({"ZL", "ZJ", "QCD", "W"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"8TeV"}).bin_id({3, 5}).process({"W"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"7TeV"}).bin_id({6}).process({"ZL", "ZJ", "W", "ZTT"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"8TeV"}).bin_id({6}).process({"ZL", "ZJ", "W"}), cb); + bbb.MergeAndAdd(cb_et.cp().era({"8TeV"}).bin_id({7}).process({"ZL", "ZJ", "W", "ZTT"}), cb); + + ch::CombineHarvester cb_mt = cb.cp().channel({"mt"}); + bbb.MergeAndAdd(cb_mt.cp().era({"7TeV"}).bin_id({1, 2, 3, 4}).process({"W", "QCD"}), cb); + bbb.MergeAndAdd(cb_mt.cp().era({"8TeV"}).bin_id({1, 2, 3, 4}).process({"W", "QCD"}), cb); + bbb.MergeAndAdd(cb_mt.cp().era({"7TeV"}).bin_id({5}).process({"W"}), cb); + bbb.MergeAndAdd(cb_mt.cp().era({"7TeV"}).bin_id({6}).process({"W", "ZTT"}), cb); + bbb.MergeAndAdd(cb_mt.cp().era({"8TeV"}).bin_id({5, 6}).process({"W"}), cb); + bbb.MergeAndAdd(cb_mt.cp().era({"8TeV"}).bin_id({7}).process({"W", "ZTT"}), cb); + + ch::CombineHarvester cb_em = cb.cp().channel({"em"}); + bbb.MergeAndAdd(cb_em.cp().era({"7TeV"}).bin_id({1, 3}).process({"Fakes"}), cb); + bbb.MergeAndAdd(cb_em.cp().era({"8TeV"}).bin_id({1, 3}).process({"Fakes"}), cb); + bbb.MergeAndAdd(cb_em.cp().era({"7TeV"}).bin_id({4}).process({"Fakes", "EWK", "Ztt"}), cb); + bbb.MergeAndAdd(cb_em.cp().era({"8TeV"}).bin_id({5}).process({"Fakes", "EWK", "Ztt"}), cb); + bbb.MergeAndAdd(cb_em.cp().era({"8TeV"}).bin_id({4}).process({"Fakes", "EWK"}), cb); + + ch::CombineHarvester cb_tt = cb.cp().channel({"tt"}); + bbb.MergeAndAdd(cb_tt.cp().era({"8TeV"}).bin_id({0, 1, 2}).process({"ZTT", "QCD"}), cb); + + bbb.SetAddThreshold(0.); // ee and mm use a different threshold + ch::CombineHarvester cb_ll = cb.cp().channel({"ee", "mm"}); + bbb.MergeAndAdd(cb_ll.cp().era({"7TeV"}).bin_id({1, 3, 4}).process({"ZTT", "ZEE", "ZMM", "TTJ"}), cb); + bbb.MergeAndAdd(cb_ll.cp().era({"8TeV"}).bin_id({1, 3, 4}).process({"ZTT", "ZEE", "ZMM", "TTJ"}), cb); cout << ">> Setting standardised bin names...\n"; ch::SetStandardBinNames(cb);