@@ -133,132 +130,132 @@
25 flags_[
"workspace-uuid-recycle"] =
true;
26 flags_[
"import-parameter-err"] =
true;
27 flags_[
"filters-use-regex"] =
false;
-
-
-
-
-
-
-
-
-
-
-
-
39 swap(first.obs_, second.obs_);
-
40 swap(first.procs_, second.procs_);
-
41 swap(first.systs_, second.systs_);
-
42 swap(first.params_, second.params_);
-
43 swap(first.wspaces_, second.wspaces_);
-
44 swap(first.verbosity_, second.verbosity_);
-
45 swap(first.flags_, second.flags_);
-
46 swap(first.post_lines_, second.post_lines_);
-
47 swap(first.log_, second.log_);
-
48 swap(first.auto_stats_settings_, second.auto_stats_settings_);
-
-
-
-
-
-
-
55 params_(other.params_),
-
56 wspaces_(other.wspaces_),
-
-
58 auto_stats_settings_(other.auto_stats_settings_),
-
59 post_lines_(other.post_lines_),
-
60 verbosity_(other.verbosity_),
-
-
-
-
-
+
28 flags_[
"fix-rooconstvar"] =
false;
+
+
+
+
+
+
+
+
+
+
+
+
40 swap(first.obs_, second.obs_);
+
41 swap(first.procs_, second.procs_);
+
42 swap(first.systs_, second.systs_);
+
43 swap(first.params_, second.params_);
+
44 swap(first.wspaces_, second.wspaces_);
+
45 swap(first.verbosity_, second.verbosity_);
+
46 swap(first.flags_, second.flags_);
+
47 swap(first.post_lines_, second.post_lines_);
+
48 swap(first.log_, second.log_);
+
49 swap(first.auto_stats_settings_, second.auto_stats_settings_);
+
+
+
+
+
+
+
56 params_(other.params_),
+
57 wspaces_(other.wspaces_),
+
+
59 auto_stats_settings_(other.auto_stats_settings_),
+
60 post_lines_(other.post_lines_),
+
61 verbosity_(other.verbosity_),
+
+
+
+
-
-
68 auto it = flags_.find(flag);
-
69 if (it != flags_.end()) {
-
70 FNLOG(std::cout) <<
"Changing value of flag \"" << it->first <<
"\" from " << it->second <<
" to " << value <<
"\n";
-
-
-
73 FNLOG(std::cout) <<
"Created new flag \"" << flag <<
"\" with value " << value <<
"\n";
-
-
-
-
-
-
79 auto it = flags_.find(flag);
-
80 if (it != flags_.end()) {
-
-
-
83 throw std::runtime_error(
FNERROR(
"Flag " + flag +
" is not defined"));
-
-
-
-
-
-
-
-
91 cpy.obs_.resize(obs_.size());
-
92 cpy.procs_.resize(procs_.size());
-
93 cpy.systs_.resize(systs_.size());
-
-
95 cpy.verbosity_ = verbosity_;
-
96 cpy.post_lines_ = post_lines_;
-
-
-
-
100 std::map<RooAbsData const*, RooAbsData *> dat_map;
-
101 std::map<RooAbsReal const*, RooAbsReal *> pdf_map;
-
102 std::map<RooRealVar const*, RooRealVar *> var_map;
-
103 std::map<RooAbsReal const*, RooAbsReal *> fun_map;
-
-
-
106 for (
auto const& it : wspaces_) {
-
107 std::string ws_name = it.first;
-
108 RooWorkspace
const& o_wsp = *(it.second.get());
-
109 cpy.wspaces_.insert({ws_name, std::make_shared<RooWorkspace>(o_wsp)});
-
110 RooWorkspace & n_wsp = *(cpy.wspaces_.at(ws_name).get());
-
-
112 std::list<RooAbsData *> o_dat = o_wsp.allData();
-
113 RooArgSet
const& o_pdf = o_wsp.allPdfs();
-
114 RooArgSet
const& o_var = o_wsp.allVars();
-
115 RooArgSet
const& o_fun = o_wsp.allFunctions();
-
-
117 std::list<RooAbsData *> n_dat = n_wsp.allData();
-
118 RooArgSet
const& n_pdf = n_wsp.allPdfs();
-
119 RooArgSet
const& n_var = n_wsp.allVars();
-
120 RooArgSet
const& n_fun = n_wsp.allFunctions();
-
-
122 auto o_dat_it = o_dat.begin();
-
123 auto n_dat_it = n_dat.begin();
-
124 for (; o_dat_it != o_dat.end(); ++o_dat_it, ++n_dat_it) {
-
125 dat_map[*o_dat_it] = *n_dat_it;
-
-
-
128 auto o_pdf_it = o_pdf.createIterator();
-
129 auto n_pdf_it = n_pdf.createIterator();
-
-
131 RooAbsPdf *o_pdf_ptr = static_cast<RooAbsPdf*>(**o_pdf_it);
-
132 RooAbsPdf *n_pdf_ptr = static_cast<RooAbsPdf*>(**n_pdf_it);
-
133 if (o_pdf_ptr && n_pdf_ptr) pdf_map[o_pdf_ptr] = n_pdf_ptr;
-
-
135 }
while (o_pdf_it->Next());
-
-
137 auto o_var_it = o_var.createIterator();
-
138 auto n_var_it = n_var.createIterator();
-
-
140 RooRealVar *o_var_ptr = static_cast<RooRealVar*>(**o_var_it);
-
141 RooRealVar *n_var_ptr = static_cast<RooRealVar*>(**n_var_it);
-
142 if (o_var_ptr && n_var_ptr) var_map[o_var_ptr] = n_var_ptr;
-
-
144 }
while (o_var_it->Next());
+
+
+
69 auto it = flags_.find(flag);
+
70 if (it != flags_.end()) {
+
71 FNLOG(std::cout) <<
"Changing value of flag \"" << it->first <<
"\" from " << it->second <<
" to " << value <<
"\n";
+
+
+
74 FNLOG(std::cout) <<
"Created new flag \"" << flag <<
"\" with value " << value <<
"\n";
+
+
+
+
+
+
80 auto it = flags_.find(flag);
+
81 if (it != flags_.end()) {
+
+
+
84 throw std::runtime_error(
FNERROR(
"Flag " + flag +
" is not defined"));
+
+
+
+
+
+
+
+
92 cpy.obs_.resize(obs_.size());
+
93 cpy.procs_.resize(procs_.size());
+
94 cpy.systs_.resize(systs_.size());
+
+
96 cpy.verbosity_ = verbosity_;
+
97 cpy.post_lines_ = post_lines_;
+
+
+
+
101 std::map<RooAbsData const*, RooAbsData *> dat_map;
+
102 std::map<RooAbsReal const*, RooAbsReal *> pdf_map;
+
103 std::map<RooRealVar const*, RooRealVar *> var_map;
+
104 std::map<RooAbsReal const*, RooAbsReal *> fun_map;
+
+
+
107 for (
auto const& it : wspaces_) {
+
108 std::string ws_name = it.first;
+
109 RooWorkspace
const& o_wsp = *(it.second.get());
+
110 cpy.wspaces_.insert({ws_name, std::make_shared<RooWorkspace>(o_wsp)});
+
111 RooWorkspace & n_wsp = *(cpy.wspaces_.at(ws_name).get());
+
+
113 std::list<RooAbsData *> o_dat = o_wsp.allData();
+
114 RooArgSet
const& o_pdf = o_wsp.allPdfs();
+
115 RooArgSet
const& o_var = o_wsp.allVars();
+
116 RooArgSet
const& o_fun = o_wsp.allFunctions();
+
+
118 std::list<RooAbsData *> n_dat = n_wsp.allData();
+
119 RooArgSet
const& n_pdf = n_wsp.allPdfs();
+
120 RooArgSet
const& n_var = n_wsp.allVars();
+
121 RooArgSet
const& n_fun = n_wsp.allFunctions();
+
+
123 auto o_dat_it = o_dat.begin();
+
124 auto n_dat_it = n_dat.begin();
+
125 for (; o_dat_it != o_dat.end(); ++o_dat_it, ++n_dat_it) {
+
126 dat_map[*o_dat_it] = *n_dat_it;
+
+
+
129 Int_t nPdf = o_pdf.getSize();
+
+
131 for (Int_t i=0; i<nPdf; ++i) {
+
+
133 RooAbsReal* o_pdf_ptr =
dynamic_cast<RooAbsReal *
>(o_pdf[i]);
+
134 RooAbsReal* n_pdf_ptr =
dynamic_cast<RooAbsReal *
>(n_pdf[i]);
+
135 if (o_pdf_ptr && n_pdf_ptr) pdf_map[o_pdf_ptr] = n_pdf_ptr;
+
+
+
138 Int_t nVar = o_var.getSize();
+
+
140 for (Int_t i=0; i<nVar; ++i) {
+
141 RooRealVar *o_var_ptr =
dynamic_cast<RooRealVar*
>(o_var[i]);
+
142 RooRealVar *n_var_ptr =
dynamic_cast<RooRealVar*
>(n_var[i]);
+
143 if (o_var_ptr && n_var_ptr) var_map[o_var_ptr] = n_var_ptr;
+
-
146 auto o_fun_it = o_fun.createIterator();
-
147 auto n_fun_it = n_fun.createIterator();
-
-
149 RooAbsReal *o_fun_ptr = static_cast<RooAbsReal*>(**o_fun_it);
-
150 RooAbsReal *n_fun_ptr = static_cast<RooAbsReal*>(**n_fun_it);
+
146 Int_t nFun = o_fun.getSize();
+
+
148 for (Int_t i=0; i<nFun; ++i) {
+
149 RooAbsReal* o_fun_ptr =
dynamic_cast<RooAbsReal *
>(o_fun[i]);
+
150 RooAbsReal* n_fun_ptr =
dynamic_cast<RooAbsReal *
>(n_fun[i]);
151 if (o_fun_ptr && n_fun_ptr) fun_map[o_fun_ptr] = n_fun_ptr;
-
-
153 }
while (o_fun_it->Next());
+
+
@@ -289,8 +286,8 @@
181 cpy.systs_[i] = std::make_shared<Systematic>(*(systs_[i]));
182 if (systs_[i]->data_u() || systs_[i]->data_d()) {
183 cpy.systs_[i]->set_data(
-
184 static_cast<RooDataHist*>(dat_map.at(systs_[i]->data_u())),
-
185 static_cast<RooDataHist*>(dat_map.at(systs_[i]->data_d())),
+
184 static_cast<RooDataHist*
>(dat_map.at(systs_[i]->data_u())),
+
185 static_cast<RooDataHist*
>(dat_map.at(systs_[i]->data_d())),
@@ -522,7 +519,7 @@
-
469 RooRealVar* norm_var = dynamic_cast<RooRealVar*>(norm);
+
469 RooRealVar* norm_var =
dynamic_cast<RooRealVar*
>(norm);
471 RooArgSet tmp_set(*norm_var);
472 ImportParameters(&tmp_set);
@@ -588,7 +585,7 @@
532 boost::replace_all(mapping.pattern,
"$PROCESS", entry->process());
533 boost::replace_all(mapping.pattern,
"$MASS", entry->mass());
-
535 mapping.
IsPdf() ? mapping.SystWorkspaceObj() : mapping.syst_pattern;
+
535 mapping.IsPdf() ? mapping.SystWorkspaceObj() : mapping.syst_pattern;
536 boost::replace_all(p_s,
"$CHANNEL", entry->bin());
537 boost::replace_all(p_s,
"$BIN", entry->bin());
538 boost::replace_all(p_s,
"$PROCESS", entry->process());
@@ -634,11 +631,11 @@
-
581 (mapping.sys_ws) ? dynamic_cast<RooDataHist*>(mapping.sys_ws->data(mapping.WorkspaceObj().c_str())) :
nullptr;
+
581 (mapping.sys_ws) ?
dynamic_cast<RooDataHist*
>(mapping.sys_ws->data(mapping.WorkspaceObj().c_str())) :
nullptr;
-
583 (mapping.sys_ws) ? dynamic_cast<RooDataHist*>(mapping.sys_ws->data(p_s_hi.c_str())):
nullptr;
+
583 (mapping.sys_ws) ?
dynamic_cast<RooDataHist*
>(mapping.sys_ws->data(p_s_hi.c_str())):
nullptr;
-
585 (mapping.sys_ws) ? dynamic_cast<RooDataHist*>(mapping.sys_ws->data(p_s_lo.c_str())):
nullptr;
+
585 (mapping.sys_ws) ?
dynamic_cast<RooDataHist*
>(mapping.sys_ws->data(p_s_lo.c_str())):
nullptr;
586 RooAbsReal* pdf = (mapping.sys_ws) ? mapping.sys_ws->function(mapping.WorkspaceObj().c_str()):
nullptr;
587 RooAbsReal* pdf_u = (mapping.sys_ws) ? mapping.sys_ws->function(p_s_hi.c_str()):
nullptr;
588 RooAbsReal* pdf_d = (mapping.sys_ws) ? mapping.sys_ws->function(p_s_lo.c_str()):
nullptr;
@@ -718,312 +715,344 @@
682 if (
GetFlag(
"workspaces-use-clone")) {
683 wspaces_[std::string(ws.GetName())] = std::shared_ptr<RooWorkspace>(
-
684 reinterpret_cast<RooWorkspace*>(ws.Clone()));
-
-
686 wspaces_[std::string(ws.GetName())] =
-
687 std::make_shared<RooWorkspace>(RooWorkspace(ws));
-
-
689 return wspaces_.at(ws.GetName());
-
-
-
-
-
-
695 throw std::runtime_error(
FNERROR(
"A different workspace with name " +
-
696 std::string(ws.GetName()) +
-
-
-
-
-
701 std::set<int> used_ints = {0};
-
702 std::string src_name(ws.GetName());
-
703 for (
auto const& it : wspaces_) {
-
704 std::string test_name(it.second->GetName());
-
705 if (test_name.find(src_name) == 0) {
-
706 std::string postfix = test_name.substr(src_name.size());
-
-
708 int number = boost::lexical_cast<int>(postfix);
-
709 used_ints.insert(number);
-
710 }
catch (boost::bad_lexical_cast & e) {
-
-
-
-
714 std::string new_name =
-
715 src_name + boost::lexical_cast<std::string>(*(used_ints.rbegin()) + 1);
-
716 FNLOGC(log(), verbosity_ >= 1) <<
"Workspace with name " << src_name
-
717 <<
" already defined, renaming to " << new_name
-
-
-
720 std::shared_ptr<RooWorkspace> new_wsp;
-
721 if (
GetFlag(
"workspaces-use-clone")) {
-
722 new_wsp = std::shared_ptr<RooWorkspace>(
-
723 reinterpret_cast<RooWorkspace*>(ws.Clone(new_name.c_str())));
-
-
725 new_wsp = std::make_shared<RooWorkspace>(RooWorkspace(ws));
-
-
727 new_wsp->SetName(new_name.c_str());
-
728 wspaces_[new_name] = new_wsp;
-
729 return wspaces_.at(new_name);
-
-
-
732 void CombineHarvester::ImportParameters(RooArgSet *vars) {
-
733 auto x = vars->createIterator();
-
-
735 RooRealVar *y = dynamic_cast<RooRealVar*>(**x);
-
-
737 if (!params_.count(y->GetName())) {
-
738 if (verbosity_ >= 1) {
-
739 log() <<
"[ImportParameters] Creating parameter from RooRealVar:\n";
-
740 y->printStream(log(), y->defaultPrintContents(0),
-
741 y->defaultPrintStyle(0),
"[ImportParameters] ");
-
-
-
744 par.set_name(y->GetName());
-
745 par.set_val(y->getVal());
-
746 if ((y->hasError() || y->hasAsymError()) &&
-
747 flags_[
"import-parameter-err"]) {
-
748 par.set_err_d(y->getErrorLo());
-
749 par.set_err_u(y->getErrorHi());
-
-
-
-
-
754 params_[y->GetName()] = std::make_shared<Parameter>(par);
-
-
-
757 log() <<
"[ImportParameters] Parameter \"" << y->GetName()
-
758 <<
"\" already exists\n";
-
-
760 Parameter *param = params_[y->GetName()].get();
-
761 std::vector<RooRealVar *> & p_vars = param->vars();
-
762 if (std::find(p_vars.begin(), p_vars.end(), y) == p_vars.end()) {
-
-
-
765 log() <<
"[ImportParameters] Parameter now stores " << p_vars.size()
-
766 <<
" link(s) to RooRealVar objects\n";
-
-
-
769 log() <<
"[ImportParameters] Parameter already stores link to "
-
770 "RooRealVar object\n";
-
-
-
-
-
-
776 RooAbsData
const* CombineHarvester::FindMatchingData(Process
const* proc) {
-
777 RooAbsData
const* data_obj =
nullptr;
-
778 for (
unsigned i = 0; i < obs_.size(); ++i) {
-
779 if (proc->bin() == obs_[i]->bin() &&
-
780 proc->bin_id() == obs_[i]->bin_id()) {
-
781 data_obj = obs_[i]->data();
-
-
-
-
-
-
787 ch::Parameter* CombineHarvester::SetupRateParamVar(std::string
const& name,
-
788 double val,
bool is_ext_arg) {
-
789 RooWorkspace *ws =
nullptr;
-
790 if (!wspaces_.count(
"_rateParams")) {
-
791 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
-
-
793 ws = wspaces_.at(
"_rateParams").get();
-
-
-
796 RooRealVar *var = ws->var(name.c_str());
-
-
-
799 RooRealVar tmp_var(name.c_str(), name.c_str(), 0);
-
-
801 var = ws->var(name.c_str());
-
802 FNLOGC(log(), verbosity_ > 1)
-
803 <<
"Created new RooRealVar for rateParam: " << name <<
"\n";
-
804 if (verbosity_ > 1) var->Print();
-
-
806 FNLOGC(log(), verbosity_ > 1)
-
807 <<
"Reusing existing RooRealVar for rateParam: " << name <<
"\n";
-
-
809 if (is_ext_arg) var->setAttribute(
"extArg");
-
810 Parameter * param =
nullptr;
-
811 if (!params_.count(name)) {
-
812 params_[name] = std::make_shared<Parameter>(Parameter());
-
813 param = params_.at(name).get();
-
814 param->set_name(name);
-
815 param->set_err_u(0.);
-
816 param->set_err_d(0.);
-
-
818 param = params_.at(name).get();
-
-
-
821 bool var_in_par =
false;
-
822 for (
auto const& ptr : param->vars()) {
-
-
-
-
-
-
-
829 param->vars().push_back(var);
-
-
-
832 FNLOGC(log(), verbosity_ > 1) <<
"Updating parameter value from "
-
-
-
835 if (verbosity_ > 1) log() <<
" to " << param->val() <<
"\n";
-
-
-
-
839 void CombineHarvester::SetupRateParamFunc(std::string
const& name,
-
840 std::string
const& formula,
-
841 std::string
const& pars) {
-
842 RooWorkspace *ws =
nullptr;
-
843 if (!wspaces_.count(
"_rateParams")) {
-
844 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
-
-
846 ws = wspaces_.at(
"_rateParams").get();
-
-
848 RooAbsReal *form = ws->function(name.c_str());
-
-
-
851 RooFormulaVar formularvar(name.c_str(), name.c_str(),
-
-
853 RooArgList(ws->argSet(pars.c_str())));
-
854 ws->import(formularvar);
-
855 form = ws->function(name.c_str());
-
856 FNLOGC(log(), verbosity_ > 1)
-
857 <<
"Created new RooFormulaVar for rateParam: " << name <<
"\n";
-
858 if (verbosity_ > 1) form->Print();
-
-
-
-
862 void CombineHarvester::SetupRateParamWspObj(std::string
const& name,
-
863 std::string
const& obj,
bool is_ext_arg) {
-
864 RooWorkspace *ws =
nullptr;
-
865 if (!wspaces_.count(
"_rateParams")) {
-
866 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
-
-
868 ws = wspaces_.at(
"_rateParams").get();
-
-
870 ws->import((obj+
":"+name).c_str(), RooFit::RecycleConflictNodes());
-
871 ws->arg(name.c_str())->setStringAttribute(
"wspSource", obj.c_str());
-
872 if (is_ext_arg) ws->arg(name.c_str())->setAttribute(
"extArg");
-
-
-
-
-
877 target.auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, sig, mode);
-
-
-
-
-
882 auto it = auto_stats_settings_.find(oldname);
-
883 if (it != auto_stats_settings_.end()) {
-
884 auto_stats_settings_[newname] = it->second;
-
885 auto_stats_settings_.erase(it);
-
-
-
-
-
890 std::set<std::string> result;
-
891 for (
auto const& it : auto_stats_settings_) {
-
892 result.insert(it.first);
-
-
-
-
-
-
-
-
-
-
-
-
904 double r=0.;
bool found=
false;
-
905 for (
auto & item : wspaces_) {
-
906 if (item.second.get()->var(name.c_str())) {
-
907 if (found) std::cout<<
"WARNING-DUPLICATE: ALREADY FOUND "<<name<<
"in an other ws"<<r<<std::endl;
-
908 r=item.second.get()->var(name.c_str())->getVal();
-
-
-
-
-
-
-
-
916 for (
auto & item : wspaces_) {
-
917 if (item.second.get()->var(name.c_str())){
-
918 if ( item.second.get()->var(name.c_str())->getMin() >value) item.second.get()->var(name.c_str())->setMin(value-0.001);
-
919 if ( item.second.get()->var(name.c_str())->getMax() <value) item.second.get()->var(name.c_str())->setMax(value+0.001);
-
920 item.second.get()->var(name.c_str())->setVal(value);
-
-
-
-
-
-
-
927 for (
auto & item : wspaces_) {
-
928 if ( (wsName==
"" or item.first == wsName ) and item.second.get()->var(name.c_str())){
-
929 item.second.get()->var(name.c_str())->SetName(newName.c_str());
-
-
-
-
+
684 reinterpret_cast<RooWorkspace*
>(ws.Clone(ws.GetName()))
+
+
+
687 wspaces_[std::string(ws.GetName())] =
+
688 std::make_shared<RooWorkspace>(RooWorkspace(ws));
+
+
690 return wspaces_.at(ws.GetName());
+
+
+
+
+
+
696 throw std::runtime_error(
FNERROR(
"A different workspace with name " +
+
697 std::string(ws.GetName()) +
+
+
+
+
+
702 std::set<int> used_ints = {0};
+
703 std::string src_name(ws.GetName());
+
704 for (
auto const& it : wspaces_) {
+
705 std::string test_name(it.second->GetName());
+
706 if (test_name.find(src_name) == 0) {
+
707 std::string postfix = test_name.substr(src_name.size());
+
+
709 int number = boost::lexical_cast<int>(postfix);
+
710 used_ints.insert(number);
+
711 }
catch (boost::bad_lexical_cast & e) {
+
+
+
+
715 std::string new_name =
+
716 src_name + boost::lexical_cast<std::string>(*(used_ints.rbegin()) + 1);
+
717 FNLOGC(log(), verbosity_ >= 1) <<
"Workspace with name " << src_name
+
718 <<
" already defined, renaming to " << new_name
+
+
+
721 std::shared_ptr<RooWorkspace> new_wsp;
+
722 if (
GetFlag(
"workspaces-use-clone")) {
+
723 new_wsp = std::shared_ptr<RooWorkspace>(
+
724 reinterpret_cast<RooWorkspace*
>(ws.Clone(new_name.c_str()))
+
+
+
727 new_wsp = std::make_shared<RooWorkspace>(RooWorkspace(ws));
+
+
729 new_wsp->SetName(new_name.c_str());
+
730 wspaces_[new_name] = new_wsp;
+
731 return wspaces_.at(new_name);
+
+
+
734 void CombineHarvester::ImportParameters(RooArgSet *vars) {
+
735 for (RooAbsArg *x : *vars) {
+
736 RooRealVar *y =
dynamic_cast<RooRealVar*
>(x);
+
+
738 if (!params_.count(y->GetName())) {
+
739 if (verbosity_ >= 1) {
+
740 log() <<
"[ImportParameters] Creating parameter from RooRealVar:\n";
+
741 y->printStream(log(), y->defaultPrintContents(0),
+
742 y->defaultPrintStyle(0),
"[ImportParameters] ");
+
+
+
745 par.set_name(y->GetName());
+
746 par.set_val(y->getVal());
+
747 if ((y->hasError() || y->hasAsymError()) &&
+
748 flags_[
"import-parameter-err"]) {
+
749 par.set_err_d(y->getErrorLo());
+
750 par.set_err_u(y->getErrorHi());
+
+
+
+
+
755 params_[y->GetName()] = std::make_shared<Parameter>(par);
+
+
+
758 log() <<
"[ImportParameters] Parameter \"" << y->GetName()
+
759 <<
"\" already exists\n";
+
+
761 Parameter *param = params_[y->GetName()].get();
+
762 std::vector<RooRealVar *> & p_vars = param->vars();
+
763 if (std::find(p_vars.begin(), p_vars.end(), y) == p_vars.end()) {
+
+
+
766 log() <<
"[ImportParameters] Parameter now stores " << p_vars.size()
+
767 <<
" link(s) to RooRealVar objects\n";
+
+
+
770 log() <<
"[ImportParameters] Parameter already stores link to "
+
771 "RooRealVar object\n";
+
+
+
+
+
+
777 RooAbsData
const* CombineHarvester::FindMatchingData(Process
const* proc) {
+
778 RooAbsData
const* data_obj =
nullptr;
+
779 for (
unsigned i = 0; i < obs_.size(); ++i) {
+
780 if (proc->bin() == obs_[i]->bin() &&
+
781 proc->bin_id() == obs_[i]->bin_id()) {
+
782 data_obj = obs_[i]->data();
+
+
+
+
+
+
788 ch::Parameter* CombineHarvester::SetupRateParamVar(std::string
const& name,
+
789 double val,
bool is_ext_arg) {
+
790 RooWorkspace *ws =
nullptr;
+
791 if (!wspaces_.count(
"_rateParams")) {
+
792 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
+
+
794 ws = wspaces_.at(
"_rateParams").get();
+
+
+
797 RooRealVar *var = ws->var(name.c_str());
+
+
+
800 RooRealVar tmp_var(name.c_str(), name.c_str(), 0);
+
+
802 var = ws->var(name.c_str());
+
803 FNLOGC(log(), verbosity_ > 1)
+
804 <<
"Created new RooRealVar for rateParam: " << name <<
"\n";
+
805 if (verbosity_ > 1) var->Print();
+
+
807 FNLOGC(log(), verbosity_ > 1)
+
808 <<
"Reusing existing RooRealVar for rateParam: " << name <<
"\n";
+
+
810 if (is_ext_arg) var->setAttribute(
"extArg");
+
811 Parameter * param =
nullptr;
+
812 if (!params_.count(name)) {
+
813 params_[name] = std::make_shared<Parameter>(Parameter());
+
814 param = params_.at(name).get();
+
815 param->set_name(name);
+
816 param->set_err_u(0.);
+
817 param->set_err_d(0.);
+
+
819 param = params_.at(name).get();
+
+
+
822 bool var_in_par =
false;
+
823 for (
auto const& ptr : param->vars()) {
+
+
+
+
+
+
+
830 param->vars().push_back(var);
+
+
+
833 FNLOGC(log(), verbosity_ > 1) <<
"Updating parameter value from "
+
+
+
836 if (verbosity_ > 1) log() <<
" to " << param->val() <<
"\n";
+
+
+
+
840 void CombineHarvester::SetupRateParamFunc(std::string
const& name,
+
841 std::string
const& formula,
+
842 std::string
const& pars) {
+
843 RooWorkspace *ws =
nullptr;
+
844 if (!wspaces_.count(
"_rateParams")) {
+
845 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
+
+
847 ws = wspaces_.at(
"_rateParams").get();
+
+
849 RooAbsReal *form = ws->function(name.c_str());
+
+
+
852 RooFormulaVar formularvar(name.c_str(), name.c_str(),
+
+
854 RooArgList(ws->argSet(pars.c_str())));
+
855 ws->import(formularvar);
+
856 form = ws->function(name.c_str());
+
857 FNLOGC(log(), verbosity_ > 1)
+
858 <<
"Created new RooFormulaVar for rateParam: " << name <<
"\n";
+
859 if (verbosity_ > 1) form->Print();
+
+
+
+
863 void CombineHarvester::SetupRateParamWspObj(std::string
const& name,
+
864 std::string
const& obj,
bool is_ext_arg) {
+
865 RooWorkspace *ws =
nullptr;
+
866 if (!wspaces_.count(
"_rateParams")) {
+
867 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
+
+
869 ws = wspaces_.at(
"_rateParams").get();
+
+
871 ws->import((obj+
":"+name).c_str(), RooFit::RecycleConflictNodes());
+
872 ws->arg(name.c_str())->setStringAttribute(
"wspSource", obj.c_str());
+
873 if (is_ext_arg) ws->arg(name.c_str())->setAttribute(
"extArg");
+
+
+
876 bool CombineHarvester::SetupRateParamWspObjFromWsStore(std::string
const& name, std::string
const& obj,
+
877 std::map<std::string, std::shared_ptr<RooWorkspace>>
const& ws_store) {
+
878 bool found_arg_in_ws =
false;
+
879 RooWorkspace *ws =
nullptr;
+
880 if (!wspaces_.count(
"_rateParams")) {
+
881 ws = this->SetupWorkspace(RooWorkspace(
"_rateParams",
"_rateParams")).get();
+
+
883 ws = wspaces_.at(
"_rateParams").get();
+
+
885 std::vector<std::string> tokens;
+
886 boost::split(tokens, obj, boost::is_any_of(
":"));
+
+
888 std::string file_name_extArg = tokens[0];
+
889 std::string wsp_name_extArg = tokens[1];
+
+
891 RooAbsArg* warg =
nullptr;
+
892 for (
auto const& it : ws_store) {
+
893 if (boost::ends_with(it.first, file_name_extArg+wsp_name_extArg)){
+
894 warg = it.second->arg(name.c_str()) ;
+
+
+
+
+
899 if(warg && !ws->arg(name.c_str())) {
+
900 ws->import(*warg, RooFit::RecycleConflictNodes());
+
901 ws->arg(name.c_str())->setStringAttribute(
"wspSource", obj.c_str());
+
902 ws->arg(name.c_str())->setAttribute(
"extArg");
+
903 found_arg_in_ws =
true;
+
+
905 return found_arg_in_ws;
+
+
+
+
+
910 target.auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, sig, mode);
+
+
+
+
+
915 auto it = auto_stats_settings_.find(oldname);
+
916 if (it != auto_stats_settings_.end()) {
+
917 auto_stats_settings_[newname] = it->second;
+
918 auto_stats_settings_.erase(it);
+
+
+
+
+
923 std::set<std::string> result;
+
924 for (
auto const& it : auto_stats_settings_) {
+
925 result.insert(it.first);
+
+
+
+
+
+
+
+
-
@@ -127,1411 +124,1503 @@
19 #include "RooRealVar.h"
20 #include "RooFormulaVar.h"
21 #include "RooCategory.h"
-
-
-
-
-
-
-
-
-
-
31 #include "CombineHarvester/CombineTools/interface/zstr.hpp"
-
-
-
-
-
-
37 std::string parse_rules) {
-
38 boost::replace_all(parse_rules,
"$ANALYSIS",
"(?<ANALYSIS>[\\w\\.]+)");
-
39 boost::replace_all(parse_rules,
"$ERA",
"(?<ERA>[\\w\\.]+)");
-
40 boost::replace_all(parse_rules,
"$CHANNEL",
"(?<CHANNEL>[\\w\\.]+)");
-
41 boost::replace_all(parse_rules,
"$BINID",
"(?<BINID>[\\w\\.]+)");
-
42 boost::replace_all(parse_rules,
"$MASS",
"(?<MASS>[\\w\\.]+)");
-
43 boost::regex rgx(parse_rules);
-
44 boost::smatch matches;
-
45 boost::regex_search(filename, matches, rgx);
-
-
47 matches.str(
"ANALYSIS"),
-
-
49 matches.str(
"CHANNEL"),
-
50 matches.str(
"BINID").length() ?
-
51 boost::lexical_cast<int>(matches.str(
"BINID")) : 0,
-
-
-
-
-
-
57 std::string
const& analysis,
-
58 std::string
const&
era,
-
-
-
61 std::string
const& mass) {
-
62 TH1::AddDirectory(kFALSE);
-
-
-
-
-
-
-
69 std::vector<std::vector<std::string>> words;
-
70 for (
unsigned i = 0; i < lines.size(); ++i) {
-
71 boost::trim(lines[i]);
-
72 if (lines[i].size() == 0)
continue;
-
73 if (lines[i].at(0) ==
'#' || lines[i].at(0) ==
'-')
continue;
-
74 words.push_back(std::vector<std::string>());
-
75 boost::split(words.back(), lines[i], boost::is_any_of(
"\t "),
-
76 boost::token_compress_on);
-
-
-
79 std::vector<HistMapping> hist_mapping;
-
-
81 std::map<std::string, std::shared_ptr<TFile>> file_store;
-
82 std::map<std::string, std::shared_ptr<RooWorkspace>> ws_store;
-
-
84 bool start_nuisance_scan =
false;
-
-
-
-
-
-
90 std::shared_ptr<ch::Observation> single_obs =
nullptr;
-
91 std::set<std::string> bin_names;
-
-
-
94 std::map<std::string, std::set<std::string>> groups;
+
22 #include "RooConstVar.h"
+
+
+
+
+
+
+
+
+
+
32 #include "CombineHarvester/CombineTools/interface/zstr.hpp"
+
+
+
+
36 RooWorkspace* CombineHarvester::fixRooConstVar(RooWorkspace *win,
bool useRooRealVar,
bool clean)
+
+
38 if ( !
GetFlag(
"fix-rooconstvar"))
return win;
+
39 std::cout<<
"[INFO]"<<
"Running fixRooConstVar on "<<win->GetName()<<std::endl;
+
+
+
42 std::string name = win->GetName();
+
43 RooWorkspace *wout =
new RooWorkspace(name.c_str(), (std::string( win->GetTitle() ) +
"_RCVfix").c_str() ) ;
+
44 for (
auto cobj : win->components()) {
+
45 if (cobj->IsA() == RooConstVar::Class() )
+
+
47 RooConstVar* cvar =
dynamic_cast<RooConstVar*
> (cobj);
+
48 std::cout<<
"[INFO] Found RooConstVar '"<<cvar->GetName()<<
"' with value "<< cvar->getVal()<<
" in workspace "<<name<<std::endl;
+
+
+
51 RooRealVar rrv = RooRealVar(cvar->GetName(),cvar->GetTitle(),cvar->getVal());
+
+
+
54 wout ->
import ( rrv, RooFit::RecycleConflictNodes());
+
+
+
57 RooConstVar rcv = RooConstVar(cvar->GetName(),cvar->GetTitle(),cvar->getVal());
+
58 wout ->
import ( rcv, RooFit::RecycleConflictNodes());
+
+
+
61 std::cout<<
" -> imported in workspace"<<std::endl;
+
62 wout -> obj( cvar->GetName() ) -> Print(
"T");
+
+
+
+
+
+
68 wout ->
import(win->components(), RooFit::RecycleConflictNodes());
+
69 for (
auto d : win->allData() ) wout ->
import(*d);
+
70 if (clean)
delete win;
+
+
+
+
+
+
+
77 std::string parse_rules) {
+
78 boost::replace_all(parse_rules,
"$ANALYSIS",
"(?<ANALYSIS>[\\w\\.]+)");
+
79 boost::replace_all(parse_rules,
"$ERA",
"(?<ERA>[\\w\\.]+)");
+
80 boost::replace_all(parse_rules,
"$CHANNEL",
"(?<CHANNEL>[\\w\\.]+)");
+
81 boost::replace_all(parse_rules,
"$BINID",
"(?<BINID>[\\w\\.]+)");
+
82 boost::replace_all(parse_rules,
"$MASS",
"(?<MASS>[\\w\\.]+)");
+
83 boost::regex rgx(parse_rules);
+
84 boost::smatch matches;
+
85 boost::regex_search(filename, matches, rgx);
+
+
87 matches.str(
"ANALYSIS"),
+
+
89 matches.str(
"CHANNEL"),
+
90 matches.str(
"BINID").length() ?
+
91 boost::lexical_cast<int>(matches.str(
"BINID")) : 0,
+
+
+
-
-
-
98 for (
unsigned i = 0; i < words.size(); ++i) {
-
-
100 if (words[i].size() <= 1)
continue;
-
-
-
-
104 if (boost::iequals(words[i][0],
"shapes") && words[i].size() >= 5) {
-
-
-
-
-
-
-
-
112 std::size_t slash = filename.find_last_of(
'/');
-
113 if (slash != filename.npos) {
-
114 dc_path = filename.substr(0, slash) +
"/" + words[i][3];
-
-
116 dc_path = words[i][3];
-
-
118 if (!file_store.count(dc_path))
-
119 file_store[dc_path] = std::make_shared<TFile>(dc_path.c_str());
-
120 mapping.
file = file_store.at(dc_path);
-
-
122 if (words[i].size() > 5) mapping.
syst_pattern = words[i][5];
+
+
97 std::string
const& analysis,
+
98 std::string
const&
era,
+
+
+
101 std::string
const& mass) {
+
102 TH1::AddDirectory(kFALSE);
+
+
+
+
+
+
+
109 std::vector<std::vector<std::string>> words;
+
110 for (
unsigned i = 0; i < lines.size(); ++i) {
+
111 boost::trim(lines[i]);
+
112 if (lines[i].size() == 0)
continue;
+
113 if (lines[i].at(0) ==
'#' || lines[i].at(0) ==
'-')
continue;
+
114 words.push_back(std::vector<std::string>());
+
115 boost::split(words.back(), lines[i], boost::is_any_of(
"\t "),
+
116 boost::token_compress_on);
+
+
+
119 std::vector<HistMapping> hist_mapping;
+
+
121 std::map<std::string, std::shared_ptr<TFile>> file_store;
+
122 std::map<std::string, std::shared_ptr<RooWorkspace>> ws_store;
-
124 if (mapping.
IsPdf()) {
-
125 std::string store_key =
-
-
127 if (!ws_store.count(store_key)) {
-
-
129 std::shared_ptr<RooWorkspace> ptr(dynamic_cast<RooWorkspace*>(
-
-
-
132 throw std::runtime_error(
FNERROR(
"Workspace not found in file"));
-
-
134 ws_store[store_key] = ptr;
-
-
136 mapping.
ws = SetupWorkspace(*(ws_store[store_key]),
true);
-
-
-
139 std::string store_key =
-
-
141 if (!ws_store.count(store_key)) {
-
-
143 std::shared_ptr<RooWorkspace> ptr(dynamic_cast<RooWorkspace*>(
-
-
-
146 throw std::runtime_error(
FNERROR(
"Workspace not found in file"));
-
-
148 ws_store[store_key] = ptr;
-
-
150 mapping.
sys_ws = SetupWorkspace(*(ws_store[store_key]),
true);
-
-
-
-
-
-
156 if (boost::iequals(words[i][0],
"shapes") && words[i].size() == 4 &&
-
157 boost::iequals(words[i][3],
"FAKE")) {
-
-
-
-
-
-
-
-
-
166 for (
unsigned i = 0; i < words.size(); ++i) {
-
167 if (words[i].size() >= 3 &&
-
168 boost::iequals(words[i][1],
"extArg")) {
-
169 if (verbosity_ > 1) {
-
170 FNLOG(log()) <<
"Processing extArg line:\n";
-
171 for (
auto const& str : words[i]) {
-
172 log() << str <<
"\t";
-
-
-
-
-
177 bool has_range = words[i].size() == 4 && words[i][3][0] ==
'[';
-
178 std::string param_name = words[i][0];
-
179 bool is_wsp_rateparam =
false;
-
-
181 boost::lexical_cast<double>(words[i][2]);
-
182 }
catch (boost::bad_lexical_cast &) {
-
183 is_wsp_rateparam =
true;
-
-
185 if ((!is_wsp_rateparam) && (words[i].size() == 3 || has_range)) {
-
-
187 param_name, boost::lexical_cast<double>(words[i][2]),
true);
-
-
-
-
191 std::vector<std::string> tokens;
-
192 boost::split(tokens, words[i][3], boost::is_any_of(
"[],"));
-
193 if (tokens.size() == 4) {
-
194 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
-
195 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
-
196 FNLOGC(log(), verbosity_ > 1) <<
"Setting parameter range to " << words[i][2];
-
-
-
199 }
else if (words[i].size() == 3 && is_wsp_rateparam) {
-
200 SetupRateParamWspObj(param_name, words[i][2],
true);
-
-
-
-
-
-
206 for (
unsigned i = 0; i < words.size(); ++i) {
-
-
208 if (words[i].size() <= 1)
continue;
-
-
-
-
-
-
-
215 if ( boost::iequals(words[i][0],
"observation") &&
-
216 boost::iequals(words[i-1][0],
"bin") &&
-
217 words[i].size() == words[i-1].size()) {
-
218 for (
unsigned p = 1; p < words[i].size(); ++p) {
-
219 auto obs = std::make_shared<Observation>();
-
220 obs->set_bin(words[i-1][p]);
-
221 obs->set_rate(boost::lexical_cast<double>(words[i][p]));
-
-
-
-
-
-
-
228 LoadShapes(obs.get(), hist_mapping);
-
-
-
-
-
-
-
235 if (boost::iequals(words[i][0],
"observation") &&
-
236 !boost::iequals(words[i-1][0],
"bin") &&
-
237 words[i].size() == 2 &&
-
238 single_obs.get() ==
nullptr) {
-
239 for (
unsigned p = 1; p < words[i].size(); ++p) {
-
240 single_obs = std::make_shared<Observation>();
-
241 single_obs->set_bin(
"");
-
242 single_obs->set_rate(boost::lexical_cast<double>(words[i][p]));
-
-
244 single_obs->set_era(
era);
-
245 single_obs->set_channel(
channel);
-
246 single_obs->set_bin_id(
bin_id);
-
247 single_obs->set_mass(
mass);
-
-
-
-
-
-
-
-
-
-
-
258 if ( boost::iequals(words[i][0],
"rate") &&
-
259 boost::iequals(words[i-1][0],
"process") &&
-
260 boost::iequals(words[i-2][0],
"process") &&
-
261 boost::iequals(words[i-3][0],
"bin") &&
-
262 words[i].size() == words[i-1].size() &&
-
263 words[i].size() == words[i-2].size() &&
-
264 words[i].size() == words[i-3].size()) {
-
265 for (
unsigned p = 1; p < words[i].size(); ++p) {
-
266 auto proc = std::make_shared<Process>();
-
267 proc->set_bin(words[i-3][p]);
-
268 bin_names.insert(words[i-3][p]);
-
-
270 int process_id = boost::lexical_cast<int>(words[i-2][p]);
-
271 proc->set_signal(process_id <= 0);
-
272 proc->set_process(words[i-1][p]);
-
273 }
catch(boost::bad_lexical_cast &) {
-
274 int process_id = boost::lexical_cast<int>(words[i-1][p]);
-
275 proc->set_signal(process_id <= 0);
-
276 proc->set_process(words[i-2][p]);
-
-
278 proc->set_rate(boost::lexical_cast<double>(words[i][p]));
-
-
-
-
-
283 proc->set_mass(
mass);
-
-
285 LoadShapes(proc.get(), hist_mapping);
-
-
287 procs_.push_back(proc);
-
-
-
290 start_nuisance_scan =
true;
-
-
-
-
294 if (start_nuisance_scan && words[i].size() >= 4) {
-
295 if (boost::iequals(words[i][1],
"param")) {
-
296 std::string param_name = words[i][0];
-
297 if (!params_.count(param_name))
-
298 params_[param_name] = std::make_shared<Parameter>(
Parameter());
-
299 Parameter * param = params_.at(param_name).get();
-
-
301 param->
set_val(boost::lexical_cast<double>(words[i][2]));
-
302 std::size_t slash_pos = words[i][3].find(
"/");
-
303 if (slash_pos != words[i][3].npos) {
-
-
305 boost::lexical_cast<double>(words[i][3].substr(0, slash_pos)));
-
-
307 boost::lexical_cast<double>(words[i][3].substr(slash_pos+1)));
-
-
309 param->
set_err_u(+1.0 * boost::lexical_cast<double>(words[i][3]));
-
310 param->
set_err_d(-1.0 * boost::lexical_cast<double>(words[i][3]));
-
-
312 if (words[i].size() >= 5) {
-
-
314 std::vector<std::string> tokens;
-
315 boost::split(tokens, words[i][4], boost::is_any_of(
"[],"));
-
316 if (tokens.size() == 4) {
-
317 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
-
318 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
-
-
-
-
-
-
-
325 if (start_nuisance_scan && words[i].size() >= 5 &&
-
326 boost::iequals(words[i][1],
"rateParam")) {
-
327 if (verbosity_ > 1) {
-
328 FNLOG(log()) <<
"Processing rateParam line:\n";
-
329 for (
auto const& str : words[i]) {
-
330 log() << str <<
"\t";
-
-
-
+
124 bool start_nuisance_scan =
false;
+
+
+
+
+
+
130 std::shared_ptr<ch::Observation> single_obs =
nullptr;
+
131 std::set<std::string> bin_names;
+
+
+
134 std::map<std::string, std::set<std::string>> groups;
+
+
+
+
138 for (
unsigned i = 0; i < words.size(); ++i) {
+
+
140 if (words[i].size() <= 1)
continue;
+
+
+
+
144 if (boost::iequals(words[i][0],
"shapes") && words[i].size() >= 5) {
+
+
+
+
+
+
+
+
152 std::size_t slash = filename.find_last_of(
'/');
+
153 if (slash != filename.npos) {
+
154 dc_path = filename.substr(0, slash) +
"/" + words[i][3];
+
+
156 dc_path = words[i][3];
+
+
158 if (!file_store.count(dc_path))
+
159 file_store[dc_path] = std::make_shared<TFile>(dc_path.c_str());
+
160 mapping.
file = file_store.at(dc_path);
+
+
162 if (words[i].size() > 5) mapping.
syst_pattern = words[i][5];
+
+
164 if (mapping.
IsPdf()) {
+
165 std::string store_key =
+
+
167 if (!ws_store.count(store_key)) {
+
+
169 std::shared_ptr<RooWorkspace> ptr(
+
+
171 dynamic_cast<RooWorkspace*
>(gDirectory->Get(mapping.
WorkspaceName().c_str()))
+
+
+
+
175 throw std::runtime_error(
FNERROR(
"Workspace not found in file"));
+
+
177 ws_store[store_key] = ptr;
+
+
179 mapping.
ws = SetupWorkspace(*(ws_store[store_key]),
true);
+
+
+
182 std::string store_key =
+
+
184 if (!ws_store.count(store_key)) {
+
+
186 std::shared_ptr<RooWorkspace> ptr(
+
+
188 dynamic_cast<RooWorkspace*
>(gDirectory->Get(mapping.
SystWorkspaceName().c_str()))
+
+
+
+
192 throw std::runtime_error(
FNERROR(
"Workspace not found in file"));
+
+
194 ws_store[store_key] = ptr;
+
+
196 mapping.
sys_ws = SetupWorkspace(*(ws_store[store_key]),
true);
+
+
+
+
+
+
202 if (boost::iequals(words[i][0],
"shapes") && words[i].size() == 4 &&
+
203 boost::iequals(words[i][3],
"FAKE")) {
+
+
+
+
+
+
+
+
+
212 for (
unsigned i = 0; i < words.size(); ++i) {
+
213 if (words[i].size() >= 3 &&
+
214 boost::iequals(words[i][1],
"extArg")) {
+
215 if (verbosity_ > 1) {
+
216 FNLOG(log()) <<
"Processing extArg line:\n";
+
217 for (
auto const& str : words[i]) {
+
218 log() << str <<
"\t";
+
+
+
+
+
223 bool has_range = words[i].size() == 4 && words[i][3][0] ==
'[';
+
224 std::string param_name = words[i][0];
+
225 bool is_wsp_rateparam =
false;
+
+
227 boost::lexical_cast<double>(words[i][2]);
+
228 }
catch (boost::bad_lexical_cast &) {
+
229 is_wsp_rateparam =
true;
+
+
231 if ((!is_wsp_rateparam) && (words[i].size() == 3 || has_range)) {
+
+
233 param_name, boost::lexical_cast<double>(words[i][2]),
true);
+
+
+
+
237 std::vector<std::string> tokens;
+
238 boost::split(tokens, words[i][3], boost::is_any_of(
"[],"));
+
239 if (tokens.size() == 4) {
+
240 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
+
241 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
+
242 FNLOGC(log(), verbosity_ > 1) <<
"Setting parameter range to " << words[i][2];
+
+
+
245 }
else if (words[i].size() == 3 && is_wsp_rateparam) {
+
246 if (!SetupRateParamWspObjFromWsStore(param_name, words[i][2], ws_store)) {
+
247 SetupRateParamWspObj(param_name, words[i][2],
true);
+
+
+
+
+
+
+
254 for (
unsigned i = 0; i < words.size(); ++i) {
+
+
256 if (words[i].size() <= 1)
continue;
+
+
+
+
+
+
+
263 if ( boost::iequals(words[i][0],
"observation") &&
+
264 boost::iequals(words[i-1][0],
"bin") &&
+
265 words[i].size() == words[i-1].size()) {
+
266 for (
unsigned p = 1; p < words[i].size(); ++p) {
+
267 auto obs = std::make_shared<Observation>();
+
268 obs->set_bin(words[i-1][p]);
+
269 obs->set_rate(boost::lexical_cast<double>(words[i][p]));
+
+
+
+
+
+
+
276 LoadShapes(obs.get(), hist_mapping);
+
+
+
+
+
+
+
283 if (boost::iequals(words[i][0],
"observation") &&
+
284 !boost::iequals(words[i-1][0],
"bin") &&
+
285 words[i].size() == 2 &&
+
286 single_obs.get() ==
nullptr) {
+
287 for (
unsigned p = 1; p < words[i].size(); ++p) {
+
288 single_obs = std::make_shared<Observation>();
+
289 single_obs->set_bin(
"");
+
290 single_obs->set_rate(boost::lexical_cast<double>(words[i][p]));
+
+
292 single_obs->set_era(
era);
+
293 single_obs->set_channel(
channel);
+
294 single_obs->set_bin_id(
bin_id);
+
295 single_obs->set_mass(
mass);
+
+
+
+
+
+
+
+
+
+
+
306 if ( boost::iequals(words[i][0],
"rate") &&
+
307 boost::iequals(words[i-1][0],
"process") &&
+
308 boost::iequals(words[i-2][0],
"process") &&
+
309 boost::iequals(words[i-3][0],
"bin") &&
+
310 words[i].size() == words[i-1].size() &&
+
311 words[i].size() == words[i-2].size() &&
+
312 words[i].size() == words[i-3].size()) {
+
313 for (
unsigned p = 1; p < words[i].size(); ++p) {
+
314 auto proc = std::make_shared<Process>();
+
315 proc->set_bin(words[i-3][p]);
+
316 bin_names.insert(words[i-3][p]);
+
+
318 int process_id = boost::lexical_cast<int>(words[i-2][p]);
+
319 proc->set_signal(process_id <= 0);
+
320 proc->set_process(words[i-1][p]);
+
321 }
catch(boost::bad_lexical_cast &) {
+
322 int process_id = boost::lexical_cast<int>(words[i-1][p]);
+
323 proc->set_signal(process_id <= 0);
+
324 proc->set_process(words[i-2][p]);
+
+
326 proc->set_rate(boost::lexical_cast<double>(words[i][p]));
+
+
+
+
+
331 proc->set_mass(
mass);
+
+
333 LoadShapes(proc.get(), hist_mapping);
-
335 bool has_range = words[i].size() == 6 && words[i][5][0] ==
'[';
-
336 std::string param_name = words[i][0];
-
-
-
-
-
341 bool is_wsp_rateparam =
false;
-
-
343 boost::lexical_cast<double>(words[i][4]);
-
344 }
catch (boost::bad_lexical_cast &) {
-
345 is_wsp_rateparam =
true;
-
-
347 if ((!is_wsp_rateparam) && (words[i].size() == 5 || has_range)) {
-
-
349 param_name, boost::lexical_cast<double>(words[i][4]));
-
-
-
-
353 std::vector<std::string> tokens;
-
354 boost::split(tokens, words[i][5], boost::is_any_of(
"[],"));
-
355 if (tokens.size() == 4) {
-
356 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
-
357 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
-
358 FNLOGC(log(), verbosity_ > 1) <<
"Setting parameter range to " << words[i][5];
-
-
-
361 }
else if (words[i].size() == 6 && !has_range) {
-
362 SetupRateParamFunc(param_name, words[i][4], words[i][5]);
-
363 }
else if (words[i].size() == 5 && is_wsp_rateparam) {
-
364 SetupRateParamWspObj(param_name, words[i][4]);
-
-
366 for (
unsigned p = 1; p < words[r].size(); ++p) {
-
367 bool matches_bin =
false;
-
368 bool matches_proc =
false;
-
-
-
371 std::string
bin = words[r-3][p];
-
-
-
374 process_id = boost::lexical_cast<int>(words[r-2][p]);
-
-
376 }
catch(boost::bad_lexical_cast &) {
-
377 process_id = boost::lexical_cast<int>(words[r-1][p]);
-
-
-
-
381 if (words[i][2] ==
"*" || fnmatch(words[i][2].c_str(),
bin.c_str(), 0) == 0) {
-
-
-
-
385 if (words[i][3] ==
"*" || fnmatch(words[i][3].c_str(),
process.c_str(), 0) == 0) {
-
-
-
388 if (!matches_bin || !matches_proc)
continue;
-
389 auto sys = std::make_shared<Systematic>();
-
-
391 sys->set_signal(process_id <= 0);
-
-
393 sys->set_name(param_name);
-
394 sys->set_type(
"rateParam");
-
-
-
-
-
-
400 systs_.push_back(sys);
-
-
-
-
-
405 if (start_nuisance_scan && words[i].size() >= 4 &&
-
406 boost::iequals(words[i][1],
"group")) {
-
407 if (!groups.count(words[i][0])) {
-
408 groups[words[i][0]] = std::set<std::string>();
-
-
410 for (
unsigned ig = 3; ig < words[i].size(); ++ig) {
-
411 groups[words[i][0]].insert(words[i][ig]);
-
-
-
-
-
416 if (start_nuisance_scan && words[i].size() >= 3 &&
-
417 boost::iequals(words[i][1],
"autoMCStats")) {
-
418 std::vector<std::string> for_bins;
-
419 if (words[i][0] ==
"*") {
-
-
-
422 for_bins.push_back(words[i][0]);
-
-
424 for (
auto const&
bin : for_bins) {
-
425 double thresh = boost::lexical_cast<double>(words[i][2]);
-
426 if (words[i].size() == 3) {
-
427 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh);
-
428 }
else if (words[i].size() == 4) {
-
429 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, boost::lexical_cast<int>(words[i][3]));
-
-
431 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, boost::lexical_cast<int>(words[i][3]), boost::lexical_cast<int>(words[i][4]));
-
-
-
-
-
436 if (start_nuisance_scan && words[i].size()-1 == words[r].size() && !boost::iequals(words[i][1],
"autoMCStats")) {
-
437 for (
unsigned p = 2; p < words[i].size(); ++p) {
-
438 if (words[i][p] ==
"-")
continue;
-
439 auto sys = std::make_shared<Systematic>();
-
440 sys->set_bin(words[r-3][p-1]);
-
-
442 int process_id = boost::lexical_cast<int>(words[r-2][p-1]);
-
443 sys->set_signal(process_id <= 0);
-
444 sys->set_process(words[r-1][p-1]);
-
445 }
catch(boost::bad_lexical_cast &) {
-
446 int process_id = boost::lexical_cast<int>(words[r-1][p-1]);
-
447 sys->set_signal(process_id <= 0);
-
448 sys->set_process(words[r-2][p-1]);
-
-
450 sys->set_name(words[i][0]);
-
451 std::string type = words[i][1];
-
452 if (!
contains(std::vector<std::string>{
"shape",
"shape?",
"shapeN2",
"shapeU",
"lnN",
"lnU"},
-
-
454 throw std::runtime_error(
-
455 FNERROR(
"Systematic type " + type +
" not supported"));
-
-
457 sys->set_type(words[i][1]);
-
-
-
-
-
-
-
464 std::size_t slash_pos = words[i][p].find(
"/");
-
465 if (slash_pos != words[i][p].npos) {
-
-
-
468 boost::lexical_cast<double>(words[i][p].substr(0, slash_pos)));
-
-
470 boost::lexical_cast<double>(words[i][p].substr(slash_pos+1)));
-
471 sys->set_asymm(
true);
-
-
473 sys->set_value_u(boost::lexical_cast<double>(words[i][p]));
-
474 sys->set_asymm(
false);
-
-
476 if (sys->type() ==
"shape" || sys->type() ==
"shapeN2" ||
-
477 sys->type() ==
"shapeU") {
-
478 sys->set_scale(boost::lexical_cast<double>(words[i][p]));
-
479 LoadShapes(sys.get(), hist_mapping);
-
480 }
else if (sys->type() ==
"shape?") {
-
-
-
483 LoadShapes(sys.get(), hist_mapping);
-
484 }
catch (std::exception & e) {
-
485 if (verbosity_ > 0) {
-
486 LOGLINE(log(),
"Systematic with shape? did not resolve to a shape");
-
-
-
-
490 if (!sys->shape_u() || !sys->shape_d()) {
-
491 sys->set_type(
"lnN");
-
-
493 sys->set_type(
"shape");
-
494 sys->set_scale(boost::lexical_cast<double>(words[i][p]));
-
-
-
497 if (sys->type() ==
"shape" || sys->type() ==
"shapeN2" ||
-
498 sys->type() ==
"shapeU")
-
499 sys->set_asymm(
true);
-
-
-
502 if (sys->type() ==
"lnU" || sys->type() ==
"shapeU") {
-
503 params_.at(sys->name())->set_err_d(0.);
-
504 params_.at(sys->name())->set_err_u(0.);
-
-
506 systs_.push_back(sys);
-
-
-
-
-
511 if (bin_names.size() == 1) {
-
512 single_obs->set_bin(*(bin_names.begin()));
-
513 LoadShapes(single_obs.get(), hist_mapping);
-
514 obs_.push_back(single_obs);
-
-
516 throw std::runtime_error(
FNERROR(
-
517 "Input card specifies a single observation entry without a bin, but "
-
518 "multiple bins defined elsewhere"));
-
-
-
-
-
523 for (
auto const& grp : groups) {
-
-
-
-
-
-
-
530 std::string
const& root_file) {
-
531 TFile file(root_file.c_str(),
"RECREATE");
-
-
-
-
-
-
-
-
-
-
-
542 void CombineHarvester::FillHistMappings(std::vector<HistMapping> & mappings) {
-
-
-
-
546 std::map<RooAbsData const*, RooWorkspace*> data_ws_map;
-
547 std::map<RooAbsReal const*, RooWorkspace*> pdf_ws_map;
-
548 for (
auto const& iter : wspaces_) {
-
549 auto dat = iter.second->allData();
-
-
551 data_ws_map[d] = iter.second.get();
-
-
553 RooArgSet vars = iter.second->allPdfs();
-
554 auto v = vars.createIterator();
-
-
556 RooAbsReal *y = dynamic_cast<RooAbsReal*>(**v);
-
557 if (y) pdf_ws_map[iter.second->pdf(y->GetName())] = iter.second.get();
-
-
559 RooArgSet fvars = iter.second->allFunctions();
-
560 auto fv = fvars.createIterator();
-
-
562 RooAbsReal *y = dynamic_cast<RooAbsReal*>(**fv);
-
563 if (y) pdf_ws_map[iter.second->function(y->GetName())] = iter.second.get();
-
564 }
while (fv->Next());
-
-
-
-
-
-
570 std::set<std::string> hist_bins;
-
-
572 for (
auto bin : bins) {
-
573 unsigned shape_count = std::count_if(procs_.begin(), procs_.end(),
-
574 [&](std::shared_ptr<ch::Process> p) {
-
575 return (p->bin() ==
bin && p->shape() && (!p->signal()));
-
-
577 shape_count += std::count_if(obs_.begin(), obs_.end(),
-
578 [&](std::shared_ptr<ch::Observation> p) {
-
579 return (p->bin() ==
bin && p->shape());
-
-
581 unsigned counting = std::count_if(
-
582 procs_.begin(), procs_.end(), [&](std::shared_ptr<ch::Process> p) {
-
583 return (p->bin() ==
bin && p->shape() ==
nullptr &&
-
584 p->pdf() ==
nullptr && p->
data() ==
nullptr);
-
-
586 counting += std::count_if(
-
587 obs_.begin(), obs_.end(), [&](std::shared_ptr<ch::Observation> p) {
-
588 return (p->bin() ==
bin && p->shape() ==
nullptr &&
-
589 p->
data() ==
nullptr);
-
-
-
592 if (shape_count > 0) {
-
593 mappings.emplace_back(
"*",
bin,
bin +
"/$PROCESS",
-
594 bin +
"/$PROCESS_$SYSTEMATIC");
-
595 hist_bins.insert(
bin);
-
596 }
else if (counting > 0) {
-
597 mappings.emplace_back(
"*",
bin,
"",
"");
-
598 mappings.back().is_fake =
true;
-
-
-
-
-
-
-
605 for (
auto sig_proc : sig_proc_set) {
-
-
-
608 auto masses =
Set2Vec(ch_signals.cp()
-
-
-
611 if (masses.size() != 1) {
-
612 throw std::runtime_error(
FNERROR(
"Process " + sig_proc +
" in bin " +
-
-
614 " has multiple entries with multiple "
-
615 "mass values, this is not supported"));
+
335 procs_.push_back(proc);
+
+
+
338 start_nuisance_scan =
true;
+
+
+
341 if (start_nuisance_scan && words[i].size() >= 4) {
+
342 if (boost::iequals(words[i][1],
"param")) {
+
343 std::string param_name = words[i][0];
+
344 Parameter * param = SetupRateParamVar(param_name, boost::lexical_cast<double>(words[i][2]));
+
345 param->
set_val(boost::lexical_cast<double>(words[i][2]));
+
346 std::size_t slash_pos = words[i][3].find(
"/");
+
347 if (slash_pos != words[i][3].npos) {
+
+
349 boost::lexical_cast<double>(words[i][3].substr(0, slash_pos)));
+
+
351 boost::lexical_cast<double>(words[i][3].substr(slash_pos+1)));
+
+
353 param->
set_err_u(+1.0 * boost::lexical_cast<double>(words[i][3]));
+
354 param->
set_err_d(-1.0 * boost::lexical_cast<double>(words[i][3]));
+
+
356 if (words[i].size() >= 5) {
+
+
358 std::vector<std::string> tokens;
+
359 boost::split(tokens, words[i][4], boost::is_any_of(
"[],"));
+
360 if (tokens.size() == 4) {
+
361 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
+
362 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
+
+
+
365 auto sys = std::make_shared<Systematic>();
+
366 sys->set_name(param_name);
+
367 sys->set_type(
"param");
+
+
+
+
+
+
+
374 std::string syst_str_ext = (boost::format(
"%g %g") % param->
val() % ((param->
err_u() - param->
err_d()) / 2.0)).str();
+
375 if (param->
range_d() != std::numeric_limits<double>::lowest() &&
+
376 param->
range_u() != std::numeric_limits<double>::max()) {
+
377 syst_str_ext += (boost::format(
" [%.4g,%.4g]") % param->
range_d() % param->
range_u()).str();
+
+
379 sys->set_param_str_ext(syst_str_ext);
+
380 systs_.push_back(sys);
+
+
+
+
+
385 if (start_nuisance_scan && words[i].size() >= 5 &&
+
386 boost::iequals(words[i][1],
"rateParam")) {
+
387 if (verbosity_ > 1) {
+
388 FNLOG(log()) <<
"Processing rateParam line:\n";
+
389 for (
auto const& str : words[i]) {
+
390 log() << str <<
"\t";
+
+
+
+
+
395 bool has_range = words[i].size() == 6 && words[i][5][0] ==
'[';
+
396 std::string param_name = words[i][0];
+
+
+
+
+
401 bool is_wsp_rateparam =
false;
+
+
403 boost::lexical_cast<double>(words[i][4]);
+
404 }
catch (boost::bad_lexical_cast &) {
+
405 is_wsp_rateparam =
true;
+
+
407 if ((!is_wsp_rateparam) && (words[i].size() == 5 || has_range)) {
+
+
409 param_name, boost::lexical_cast<double>(words[i][4]));
+
+
+
+
413 std::vector<std::string> tokens;
+
414 boost::split(tokens, words[i][5], boost::is_any_of(
"[],"));
+
415 if (tokens.size() == 4) {
+
416 param->
set_range_d(boost::lexical_cast<double>(tokens[1]));
+
417 param->
set_range_u(boost::lexical_cast<double>(tokens[2]));
+
418 FNLOGC(log(), verbosity_ > 1) <<
"Setting parameter range to " << words[i][5];
+
+
+
421 }
else if (words[i].size() == 6 && !has_range) {
+
422 SetupRateParamFunc(param_name, words[i][4], words[i][5]);
+
423 }
else if (words[i].size() == 5 && is_wsp_rateparam) {
+
424 SetupRateParamWspObj(param_name, words[i][4]);
+
+
426 for (
unsigned p = 1; p < words[r].size(); ++p) {
+
427 bool matches_bin =
false;
+
428 bool matches_proc =
false;
+
+
+
431 std::string
bin = words[r-3][p];
+
+
+
434 process_id = boost::lexical_cast<int>(words[r-2][p]);
+
+
436 }
catch(boost::bad_lexical_cast &) {
+
437 process_id = boost::lexical_cast<int>(words[r-1][p]);
+
+
+
+
441 if (words[i][2] ==
"*" || fnmatch(words[i][2].c_str(),
bin.c_str(), 0) == 0) {
+
+
+
+
445 if (words[i][3] ==
"*" || fnmatch(words[i][3].c_str(),
process.c_str(), 0) == 0) {
+
+
+
448 if (!matches_bin || !matches_proc)
continue;
+
449 auto sys = std::make_shared<Systematic>();
+
+
451 sys->set_signal(process_id <= 0);
+
+
453 sys->set_name(param_name);
+
454 sys->set_type(
"rateParam");
+
+
+
+
+
+
460 systs_.push_back(sys);
+
+
+
+
+
465 if (start_nuisance_scan && words[i].size() >= 4 &&
+
466 boost::iequals(words[i][1],
"group")) {
+
467 if (!groups.count(words[i][0])) {
+
468 groups[words[i][0]] = std::set<std::string>();
+
+
470 for (
unsigned ig = 3; ig < words[i].size(); ++ig) {
+
471 groups[words[i][0]].insert(words[i][ig]);
+
+
+
+
+
476 if (start_nuisance_scan && words[i].size() >= 3 &&
+
477 boost::iequals(words[i][1],
"autoMCStats")) {
+
478 std::vector<std::string> for_bins;
+
479 if (words[i][0] ==
"*") {
+
+
+
482 for_bins.push_back(words[i][0]);
+
+
484 for (
auto const&
bin : for_bins) {
+
485 double thresh = boost::lexical_cast<double>(words[i][2]);
+
486 if (words[i].size() == 3) {
+
487 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh);
+
488 }
else if (words[i].size() == 4) {
+
489 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, boost::lexical_cast<int>(words[i][3]));
+
+
491 auto_stats_settings_[
bin] = AutoMCStatsSettings(thresh, boost::lexical_cast<int>(words[i][3]), boost::lexical_cast<int>(words[i][4]));
+
+
+
+
+
496 if (start_nuisance_scan && words[i].size()-1 == words[r].size() && !boost::iequals(words[i][1],
"autoMCStats")) {
+
497 for (
unsigned p = 2; p < words[i].size(); ++p) {
+
498 if (words[i][p] ==
"-")
continue;
+
499 auto sys = std::make_shared<Systematic>();
+
500 sys->set_bin(words[r-3][p-1]);
+
+
502 int process_id = boost::lexical_cast<int>(words[r-2][p-1]);
+
503 sys->set_signal(process_id <= 0);
+
504 sys->set_process(words[r-1][p-1]);
+
505 }
catch(boost::bad_lexical_cast &) {
+
506 int process_id = boost::lexical_cast<int>(words[r-1][p-1]);
+
507 sys->set_signal(process_id <= 0);
+
508 sys->set_process(words[r-2][p-1]);
+
+
510 sys->set_name(words[i][0]);
+
511 std::string type = words[i][1];
+
512 if (!
contains(std::vector<std::string>{
"shape",
"shape?",
"shapeN2",
"shapeU",
"lnN",
"lnU"},
+
+
514 throw std::runtime_error(
+
515 FNERROR(
"Systematic type " + type +
" not supported"));
+
+
517 sys->set_type(words[i][1]);
+
+
+
+
+
+
+
524 std::size_t slash_pos = words[i][p].find(
"/");
+
525 if (slash_pos != words[i][p].npos) {
+
+
+
528 boost::lexical_cast<double>(words[i][p].substr(0, slash_pos)));
+
+
530 boost::lexical_cast<double>(words[i][p].substr(slash_pos+1)));
+
531 sys->set_asymm(
true);
+
+
533 sys->set_value_u(boost::lexical_cast<double>(words[i][p]));
+
534 sys->set_asymm(
false);
+
+
536 if (sys->type() ==
"shape" || sys->type() ==
"shapeN2" ||
+
537 sys->type() ==
"shapeU") {
+
538 sys->set_scale(boost::lexical_cast<double>(words[i][p]));
+
539 LoadShapes(sys.get(), hist_mapping);
+
540 }
else if (sys->type() ==
"shape?") {
+
+
+
543 LoadShapes(sys.get(), hist_mapping);
+
544 }
catch (std::exception & e) {
+
545 if (verbosity_ > 0) {
+
546 LOGLINE(log(),
"Systematic with shape? did not resolve to a shape");
+
+
+
+
550 if (!sys->shape_u() || !sys->shape_d()) {
+
551 sys->set_type(
"lnN");
+
+
553 sys->set_type(
"shape");
+
554 sys->set_scale(boost::lexical_cast<double>(words[i][p]));
+
+
+
557 if (sys->type() ==
"shape" || sys->type() ==
"shapeN2" ||
+
558 sys->type() ==
"shapeU")
+
559 sys->set_asymm(
true);
+
+
+
562 if (sys->type() ==
"lnU" || sys->type() ==
"shapeU") {
+
563 params_.at(sys->name())->set_err_d(0.);
+
564 params_.at(sys->name())->set_err_u(0.);
+
+
566 systs_.push_back(sys);
+
+
+
+
+
571 if (bin_names.size() == 1) {
+
572 single_obs->set_bin(*(bin_names.begin()));
+
573 LoadShapes(single_obs.get(), hist_mapping);
+
574 obs_.push_back(single_obs);
+
+
576 throw std::runtime_error(
FNERROR(
+
577 "Input card specifies a single observation entry without a bin, but "
+
578 "multiple bins defined elsewhere"));
+
+
+
+
+
583 for (
auto const& grp : groups) {
+
+
+
+
+
+
+
590 std::string
const& root_file) {
+
591 TFile file(root_file.c_str(),
"RECREATE");
+
+
+
+
+
+
+
+
+
+
+
602 void CombineHarvester::FillHistMappings(std::vector<HistMapping> & mappings) {
+
+
+
+
606 std::map<RooAbsData const*, RooWorkspace*> data_ws_map;
+
607 std::map<RooAbsReal const*, RooWorkspace*> pdf_ws_map;
+
608 for (
auto const& iter : wspaces_) {
+
609 auto dat = iter.second->allData();
+
+
611 data_ws_map[d] = iter.second.get();
+
+
613 RooArgSet vars = iter.second->allPdfs();
+
614 for (RooAbsArg *v : vars) {
+
615 RooAbsPdf *y =
dynamic_cast<RooAbsPdf*
>(v);
-
-
-
619 mappings.emplace_back(sig_proc,
bin,
bin +
"/" + sig_proc +
"$MASS",
-
620 bin +
"/" + sig_proc +
"$MASS_$SYSTEMATIC");
-
-
-
-
-
-
626 for (
auto bin : bins) {
-
-
628 for (
auto obs : ch_bin.obs_) {
-
629 if (!obs->data())
continue;
-
630 std::string obj_name = std::string(data_ws_map[obs->data()]->GetName()) +
-
631 ":" + std::string(obs->data()->GetName());
-
632 mappings.emplace_back(
"data_obs", obs->bin(), obj_name,
"");
-
-
-
635 bool prototype_ok =
false;
-
636 HistMapping prototype;
-
637 std::vector<HistMapping> full_list;
-
638 auto pmap = ch_bin.GenerateProcSystMap();
-
639 for (
unsigned i = 0; i < ch_bin.procs_.size(); ++i) {
-
-
641 if (!proc->
data() && !proc->
pdf())
continue;
-
642 std::string obj_name;
-
643 std::string obj_sys_name;
-
-
645 obj_name = std::string(proc->
data()->GetName());
-
646 boost::replace_all(obj_name, proc->
process(),
"$PROCESS");
-
647 obj_name = std::string(data_ws_map[proc->
data()]->GetName()) +
":" +
-
-
-
-
651 obj_name = std::string(proc->
pdf()->GetName());
-
652 boost::replace_all(obj_name, proc->
process(),
"$PROCESS");
-
653 obj_name = std::string(pdf_ws_map[proc->
pdf()]->GetName()) +
-
-
-
656 for (
unsigned j = 0; j < pmap[i].size(); ++j) {
-
-
-
659 obj_sys_name = std::string(sys->
data_u()->GetName());
-
660 boost::replace_all(obj_sys_name, sys->
name() +
"Up",
"$SYSTEMATIC");
-
661 boost::replace_all(obj_sys_name, sys->
process(),
"$PROCESS");
-
662 obj_sys_name = std::string(data_ws_map[sys->
data_u()]->GetName()) +
-
-
-
-
-
667 obj_sys_name = std::string(sys->
pdf_u()->GetName());
-
668 boost::replace_all(obj_sys_name, sys->
name() +
"Up",
"$SYSTEMATIC");
-
669 boost::replace_all(obj_sys_name, sys->
process(),
"$PROCESS");
-
670 obj_sys_name = std::string(pdf_ws_map[sys->
pdf_u()]->GetName()) +
-
-
-
-
+
617 if (y) pdf_ws_map[iter.second->pdf(y->GetName())] = iter.second.get();
+
+
619 RooArgSet fvars = iter.second->allFunctions();
+
620 for (RooAbsArg *fv : fvars) {
+
621 RooAbsReal *y =
dynamic_cast<RooAbsReal*
>(fv);
+
622 if (y) pdf_ws_map[iter.second->function(y->GetName())] = iter.second.get();
+
+
+
+
+
+
+
629 std::set<std::string> hist_bins;
+
+
631 for (
auto bin : bins) {
+
632 unsigned shape_count = std::count_if(procs_.begin(), procs_.end(),
+
633 [&](std::shared_ptr<ch::Process> p) {
+
634 return (p->bin() == bin && p->shape() && (!p->signal()));
+
+
636 shape_count += std::count_if(obs_.begin(), obs_.end(),
+
637 [&](std::shared_ptr<ch::Observation> p) {
+
638 return (p->bin() == bin && p->shape());
+
+
640 unsigned counting = std::count_if(
+
641 procs_.begin(), procs_.end(), [&](std::shared_ptr<ch::Process> p) {
+
642 return (p->bin() == bin && p->shape() == nullptr &&
+
643 p->pdf() == nullptr && p->data() == nullptr);
+
+
645 counting += std::count_if(
+
646 obs_.begin(), obs_.end(), [&](std::shared_ptr<ch::Observation> p) {
+
647 return (p->bin() == bin && p->shape() == nullptr &&
+
648 p->data() == nullptr);
+
+
+
651 if (shape_count > 0) {
+
652 mappings.emplace_back(
"*",
bin,
bin +
"/$PROCESS",
+
653 bin +
"/$PROCESS_$SYSTEMATIC");
+
654 hist_bins.insert(
bin);
+
655 }
else if (counting > 0) {
+
656 mappings.emplace_back(
"*",
bin,
"",
"");
+
657 mappings.back().is_fake =
true;
+
+
+
+
+
+
+
664 for (
auto sig_proc : sig_proc_set) {
+
+
+
667 auto masses =
Set2Vec(ch_signals.cp()
+
+
+
670 if (masses.size() != 1) {
+
671 throw std::runtime_error(
FNERROR(
"Process " + sig_proc +
" in bin " +
+
+
673 " has multiple entries with multiple "
+
674 "mass values, this is not supported"));
-
-
-
678 if (prototype.pattern.size() && prototype.pattern != obj_name) {
-
679 prototype_ok =
false;
+
+
+
678 mappings.emplace_back(sig_proc,
bin,
bin +
"/" + sig_proc +
"$MASS",
+
679 bin +
"/" + sig_proc +
"$MASS_$SYSTEMATIC");
-
-
682 if (prototype.syst_pattern.size() && obj_sys_name.size() &&
-
683 prototype.syst_pattern != obj_sys_name) {
-
684 prototype_ok =
false;
-
-
-
687 if (!prototype.pattern.size()) {
-
-
-
690 prototype.category =
bin;
-
691 prototype.pattern = obj_name;
-
-
693 if (!prototype.syst_pattern.size()) {
-
694 prototype.syst_pattern = obj_sys_name;
-
-
-
697 full_list.emplace_back(proc->
process(), proc->
bin(), obj_name,
-
-
-
-
-
-
-
704 if (!prototype_ok || hist_bins.count(
bin)) {
-
705 for (
auto m : full_list) {
-
706 mappings.push_back(m);
-
-
-
709 mappings.push_back(prototype);
-
-
-
-
-
-
-
-
-
-
-
720 bool is_counting =
true;
-
-
722 if (obs->
shape() !=
nullptr || obs->
data() !=
nullptr) {
-
-
-
-
-
-
728 if (proc->
shape() !=
nullptr || proc->
data() !=
nullptr ||
-
729 proc->
pdf() !=
nullptr) {
-
-
-
-
+
+
+
+
+
685 for (
auto bin : bins) {
+
+
687 for (
auto obs : ch_bin.obs_) {
+
688 if (!obs->data())
continue;
+
689 std::string obj_name = std::string(data_ws_map[obs->data()]->GetName()) +
+
690 ":" + std::string(obs->data()->GetName());
+
691 mappings.emplace_back(
"data_obs", obs->bin(), obj_name,
"");
+
+
+
694 bool prototype_ok =
false;
+
695 HistMapping prototype;
+
696 std::vector<HistMapping> full_list;
+
697 auto pmap = ch_bin.GenerateProcSystMap();
+
698 for (
unsigned i = 0; i < ch_bin.procs_.size(); ++i) {
+
+
700 if (!proc->
data() && !proc->
pdf())
continue;
+
701 std::string obj_name;
+
702 std::string obj_sys_name;
+
+
704 obj_name = std::string(proc->
data()->GetName());
+
705 boost::replace_all(obj_name, proc->
process(),
"$PROCESS");
+
706 obj_name = std::string(data_ws_map[proc->
data()]->GetName()) +
":" +
+
+
+
+
710 obj_name = std::string(proc->
pdf()->GetName());
+
711 boost::replace_all(obj_name, proc->
process(),
"$PROCESS");
+
712 obj_name = std::string(pdf_ws_map[proc->
pdf()]->GetName()) +
+
+
+
715 for (
unsigned j = 0; j < pmap[i].size(); ++j) {
+
+
+
718 obj_sys_name = std::string(sys->
data_u()->GetName());
+
719 boost::replace_all(obj_sys_name, sys->
name() +
"Up",
"$SYSTEMATIC");
+
720 boost::replace_all(obj_sys_name, sys->
process(),
"$PROCESS");
+
721 obj_sys_name = std::string(data_ws_map[sys->
data_u()]->GetName()) +
+
+
+
+
+
726 obj_sys_name = std::string(sys->
pdf_u()->GetName());
+
727 boost::replace_all(obj_sys_name, sys->
name() +
"Up",
"$SYSTEMATIC");
+
728 boost::replace_all(obj_sys_name, sys->
process(),
"$PROCESS");
+
729 obj_sys_name = std::string(pdf_ws_map[sys->
pdf_u()]->GetName()) +
+
+
+
+
-
-
736 if (!root_file.IsOpen() && !is_counting) {
-
737 throw std::runtime_error(
FNERROR(
-
738 std::string(
"Output ROOT file is not open: ") + root_file.GetName()));
-
+
+
+
737 if (prototype.pattern.size() && prototype.pattern != obj_name) {
+
738 prototype_ok =
false;
+
-
741 std::unique_ptr<std::ostream> txt_file_ptr =
nullptr;
-
-
-
744 std::string zip_ext =
".gz";
-
745 bool has_zip_ext = (name.length() >= zip_ext.length() && name.compare(name.length() - zip_ext.length(), zip_ext.length(), zip_ext) == 0);
-
-
-
748 txt_file_ptr = std::make_unique<zstr::ofstream>(name);
-
-
750 txt_file_ptr = std::make_unique<std::ofstream>(name);
-
-
752 if (txt_file_ptr->fail()) {
-
753 throw std::runtime_error(
FNERROR(
"Unable to create file: " + name));
-
+
741 if (prototype.syst_pattern.size() && obj_sys_name.size() &&
+
742 prototype.syst_pattern != obj_sys_name) {
+
743 prototype_ok =
false;
+
+
+
746 if (!prototype.pattern.size()) {
+
+
748 prototype.process =
"*";
+
749 prototype.category =
bin;
+
750 prototype.pattern = obj_name;
+
+
752 if (!prototype.syst_pattern.size()) {
+
753 prototype.syst_pattern = obj_sys_name;
+
-
756 std::ostream & txt_file = *txt_file_ptr;
-
-
-
-
-
761 std::string dashes(80,
'-');
-
-
-
-
765 std::set<std::string> sys_set;
-
766 std::set<std::string> rateparam_set;
-
-
768 if (sys->
type() ==
"rateParam") {
-
769 rateparam_set.insert(sys->
name());
-
-
771 sys_set.insert(sys->
name());
-
-
-
774 txt_file <<
"imax " <<
bin_set.size()
-
775 <<
" number of bins\n";
-
776 txt_file <<
"jmax " << proc_set.size() - 1
-
777 <<
" number of processes minus 1\n";
-
778 txt_file <<
"kmax " <<
"*" <<
" number of nuisance parameters\n";
-
779 txt_file << dashes <<
"\n";
-
-
-
782 std::vector<HistMapping> mappings;
-
783 FillHistMappings(mappings);
-
-
-
-
787 auto proc_sys_map = this->GenerateProcSystMap();
-
-
789 std::set<std::string> all_dependents_pars;
-
790 std::set<std::string> multipdf_cats;
-
791 for (
auto proc : procs_) {
-
792 if (!proc->
pdf())
continue;
-
-
797 RooAbsData
const* data_obj = FindMatchingData(proc.get());
-
-
799 RooArgSet norm_par_list;
-
-
-
-
-
-
-
-
807 RooRealVar mx(
"CMS_th1x" ,
"CMS_th1x", 0, 1);
-
808 RooArgSet tmp_set(mx);
-
-
-
-
-
-
814 RooFIter par_list_it = par_list.fwdIterator();
-
815 RooAbsArg *par_list_var =
nullptr;
-
816 while ((par_list_var = par_list_it.next())) {
-
817 if (dynamic_cast<RooRealVar*>(par_list_var)) {
-
818 all_dependents_pars.insert(par_list_var->GetName());
-
-
-
-
-
-
-
825 if (dynamic_cast<RooCategory*>(par_list_var) &&
-
826 proc->
pdf()->InheritsFrom(
"RooMultiPdf")) {
-
827 multipdf_cats.insert(par_list_var->GetName());
-
-
-
-
831 RooFIter nm_list_it = norm_par_list.fwdIterator();
-
832 RooAbsArg *nm_list_var =
nullptr;
-
833 while ((nm_list_var = nm_list_it.next())) {
-
834 if (dynamic_cast<RooRealVar*>(nm_list_var)) {
-
835 all_dependents_pars.insert(nm_list_var->GetName());
-
-
-
-
-
-
-
842 std::string file_name = root_file.GetName();
-
-
-
-
846 boost::filesystem::path root_file_path =
-
847 boost::filesystem::absolute(file_name);
-
-
849 boost::filesystem::path txt_file_path =
-
850 boost::filesystem::absolute(name).parent_path();
-
-
852 file_name =
make_relative(txt_file_path, root_file_path).string();
-
-
-
-
856 std::set<std::string> used_wsps;
-
-
858 for (
auto const& mapping : mappings) {
-
859 if (!mapping.is_fake) {
-
860 txt_file << format(
"shapes %s %s %s %s %s\n")
-
-
-
-
-
865 % mapping.syst_pattern;
-
866 if (mapping.IsPdf() || mapping.IsData()) {
-
867 used_wsps.insert(mapping.WorkspaceName());
-
868 if (mapping.syst_pattern !=
"") {
-
869 used_wsps.insert(mapping.SystWorkspaceName());
-
-
-
-
873 txt_file << format(
"shapes %s %s %s\n")
-
-
-
-
-
-
879 txt_file << dashes <<
"\n";
-
-
-
882 for (
auto ws_it : wspaces_) {
-
883 if (ws_it.first ==
"_rateParams")
continue;
-
-
885 if (!used_wsps.count(ws_it.second->GetName()))
continue;
-
-
-
-
-
-
891 if (obs_.size() > 0) {
-
-
893 for (
auto const& obs : obs_) {
-
894 txt_file << format(
"%-15s ") % obs->
bin();
-
-
896 bool add_dir = TH1::AddDirectoryStatus();
-
897 TH1::AddDirectory(
false);
-
898 std::unique_ptr<TH1> h((TH1*)(obs->
shape()->Clone()));
-
899 h->Scale(obs->
rate());
-
900 WriteHistToFile(h.get(), &root_file, mappings, obs->
bin(),
"data_obs",
-
-
902 TH1::AddDirectory(add_dir);
-
-
-
-
906 txt_file <<
"observation ";
-
-
-
-
-
911 std::string obs_fmt_int =
"%-15.1f ";
-
912 std::string obs_fmt_flt =
"%-15.4f ";
-
913 for (
auto const& obs : obs_) {
-
-
915 std::fabs(obs->
rate() - std::round(obs->
rate())) > 1E-4;
-
916 txt_file << format(
is_float ? obs_fmt_flt : obs_fmt_int)
-
-
-
-
920 txt_file << dashes <<
"\n";
-
-
-
923 unsigned sys_str_len = 14;
-
924 for (
auto const& sys : sys_set) {
-
925 if (sys.length() > sys_str_len) sys_str_len = sys.length();
-
-
927 for (
auto const& sys : all_dependents_pars) {
-
928 if (sys.length() > sys_str_len) sys_str_len = sys.length();
-
-
930 std::string sys_str_short = boost::lexical_cast<std::string>(sys_str_len);
-
931 std::string sys_str_long = boost::lexical_cast<std::string>(sys_str_len+9);
-
-
933 auto getProcLen = [](
auto const& proc) -> std::string {
-
934 std::size_t proc_len = 15;
-
935 proc_len = std::max(proc_len, proc->
process().size());
-
936 std::string proc_len_str = boost::lexical_cast<std::string>(proc_len);
-
-
+
756 full_list.emplace_back(proc->
process(), proc->
bin(), obj_name,
+
+
+
+
+
+
+
763 if (!prototype_ok || hist_bins.count(
bin)) {
+
764 for (
auto m : full_list) {
+
765 mappings.push_back(m);
+
+
+
768 mappings.push_back(prototype);
+
+
+
+
+
+
+
+
+
+
+
779 bool is_counting =
true;
+
+
781 if (obs->
shape() !=
nullptr || obs->
data() !=
nullptr) {
+
+
+
+
+
+
787 if (proc->
shape() !=
nullptr || proc->
data() !=
nullptr ||
+
788 proc->
pdf() !=
nullptr) {
+
+
+
+
+
+
+
795 if (!root_file.IsOpen() && !is_counting) {
+
796 throw std::runtime_error(
FNERROR(
+
797 std::string(
"Output ROOT file is not open: ") + root_file.GetName()));
+
+
+
800 std::unique_ptr<std::ostream> txt_file_ptr =
nullptr;
+
+
+
803 std::string zip_ext =
".gz";
+
804 bool has_zip_ext = (name.length() >= zip_ext.length() && name.compare(name.length() - zip_ext.length(), zip_ext.length(), zip_ext) == 0);
+
+
+
807 txt_file_ptr = std::make_unique<zstr::ofstream>(name);
+
+
809 txt_file_ptr = std::make_unique<std::ofstream>(name);
+
+
811 if (txt_file_ptr->fail()) {
+
812 throw std::runtime_error(
FNERROR(
"Unable to create file: " + name));
+
+
+
815 std::ostream & txt_file = *txt_file_ptr;
+
+
+
+
+
820 std::string dashes(80,
'-');
+
+
+
+
824 std::set<std::string> sys_set;
+
825 std::set<std::string> param_set;
+
826 std::set<std::string> rateparam_set;
+
+
828 if (sys->
type() ==
"rateParam") {
+
829 rateparam_set.insert(sys->name());
+
+
831 else if (sys->
type() ==
"param"){
+
832 param_set.insert(sys->name());
+
+
834 else sys_set.insert(sys->
name());
+
+
836 txt_file <<
"imax " <<
bin_set.size()
+
837 <<
" number of bins\n";
+
838 txt_file <<
"jmax " << proc_set.size() - 1
+
839 <<
" number of processes minus 1\n";
+
840 txt_file <<
"kmax " <<
"*" <<
" number of nuisance parameters\n";
+
841 txt_file << dashes <<
"\n";
+
+
+
844 std::vector<HistMapping> mappings;
+
845 FillHistMappings(mappings);
+
+
+
+
849 auto proc_sys_map = this->GenerateProcSystMap();
+
+
851 std::set<std::string> all_dependents_pars;
+
852 std::set<std::string> multipdf_cats;
+
853 for (
auto proc : procs_) {
+
854 if (!proc->
pdf())
continue;
+
+
859 RooAbsData
const* data_obj = FindMatchingData(proc.get());
+
+
861 RooArgSet norm_par_list;
+
+
+
+
+
+
+
+
869 RooRealVar mx(
"CMS_th1x" ,
"CMS_th1x", 0, 1);
+
870 RooArgSet tmp_set(mx);
+
+
+
+
+
+
876 for (RooAbsArg *par_list_var : par_list) {
+
877 if (
dynamic_cast<RooRealVar*
>(par_list_var)) {
+
878 all_dependents_pars.insert(par_list_var->GetName());
+
+
+
+
+
+
+
885 if (
dynamic_cast<RooCategory*
>(par_list_var) &&
+
886 proc->
pdf()->InheritsFrom(
"RooMultiPdf")) {
+
887 multipdf_cats.insert(par_list_var->GetName());
+
+
+
+
+
892 for (RooAbsArg *nm_list_var: norm_par_list) {
+
893 if (
dynamic_cast<RooRealVar*
>(nm_list_var)) {
+
894 all_dependents_pars.insert(nm_list_var->GetName());
+
+
+
+
+
+
+
901 std::string file_name = root_file.GetName();
+
+
+
+
905 boost::filesystem::path root_file_path =
+
906 boost::filesystem::absolute(file_name);
+
+
908 boost::filesystem::path txt_file_path =
+
909 boost::filesystem::absolute(name).parent_path();
+
+
911 file_name =
make_relative(txt_file_path, root_file_path).string();
+
+
+
+
915 std::set<std::string> used_wsps;
+
+
917 for (
auto const& mapping : mappings) {
+
918 if (!mapping.is_fake) {
+
919 txt_file << format(
"shapes %s %s %s %s %s\n")
+
+
+
+
+
924 % mapping.syst_pattern;
+
925 if (mapping.IsPdf() || mapping.IsData()) {
+
926 used_wsps.insert(mapping.WorkspaceName());
+
927 if (mapping.syst_pattern !=
"") {
+
928 used_wsps.insert(mapping.SystWorkspaceName());
+
+
+
+
932 txt_file << format(
"shapes %s %s %s\n")
+
+
+
+
+
+
938 txt_file << dashes <<
"\n";
-
940 txt_file << format(
"%-"+sys_str_long+
"s") %
"bin";
-
941 for (
auto const& proc : procs_) {
-
-
943 bool add_dir = TH1::AddDirectoryStatus();
-
944 TH1::AddDirectory(
false);
-
-
946 WriteHistToFile(h.get(), &root_file, mappings, proc->
bin(),
-
-
948 TH1::AddDirectory(add_dir);
-
-
950 txt_file << format(
"%-"+getProcLen(proc)+
"s ") % proc->
bin();
-
-
-
-
954 txt_file << format(
"%-"+sys_str_long+
"s") %
"process";
-
-
956 for (
auto const& proc : procs_) {
-
957 txt_file << format(
"%-"+getProcLen(proc)+
"s ") % proc->
process();
-
-
-
-
961 txt_file << format(
"%-"+sys_str_long+
"s") %
"process";
-
-
-
964 std::map<std::string, int> p_ids;
-
-
-
967 for (
unsigned p = 0; p < procs_.size(); ++p) {
-
968 if (!procs_[p]->signal()) {
-
969 if (p_ids.count(procs_[p]->process()) == 0) {
-
970 p_ids[procs_[p]->process()] = bkg_id;
-
-
-
-
-
975 if (p_ids[procs_[p]->
process()] <= 0)
throw std::runtime_error(
FNERROR(
"Ambiguous definition of process (" + procs_[p]->
process() +
") as both signal and background"));
-
+
+
941 for (
auto ws_it : wspaces_) {
+
942 if (ws_it.first ==
"_rateParams")
continue;
+
+
944 if (!used_wsps.count(ws_it.second->GetName()))
continue;
+
+
+
+
+
+
950 if (obs_.size() > 0) {
+
+
952 for (
auto const& obs : obs_) {
+
953 txt_file << format(
"%-15s ") % obs->
bin();
+
+
955 bool add_dir = TH1::AddDirectoryStatus();
+
956 TH1::AddDirectory(
false);
+
957 std::unique_ptr<TH1> h((TH1*)(obs->
shape()->Clone()));
+
958 h->Scale(obs->
rate());
+
959 WriteHistToFile(h.get(), &root_file, mappings, obs->
bin(),
"data_obs",
+
+
961 TH1::AddDirectory(add_dir);
+
+
+
+
965 txt_file <<
"observation ";
+
+
+
+
+
970 std::string obs_fmt_int =
"%-15.1f ";
+
971 std::string obs_fmt_flt =
"%-15.4f ";
+
972 for (
auto const& obs : obs_) {
+
+
974 std::fabs(obs->
rate() - std::round(obs->
rate())) > 1E-4;
+
975 txt_file << format(
is_float ? obs_fmt_flt : obs_fmt_int)
+
-
978 unsigned q = procs_.size() - (p + 1);
-
979 if (procs_[q]->signal()) {
-
980 if (p_ids.count(procs_[q]->process()) == 0) {
-
981 p_ids[procs_[q]->process()] = sig_id;
-
-
-
-
-
986 if (p_ids[procs_[q]->
process()] > 0)
throw std::runtime_error(
FNERROR(
"Ambiguous definition of process (" + procs_[q]->
process() +
") as both signal and background"));
-
-
-
-
990 for (
auto const& proc : procs_) {
-
991 txt_file << format(
"%-"+getProcLen(proc)+
"s ") % p_ids[proc->
process()];
-
-
-
-
-
996 txt_file << format(
"%-"+sys_str_long+
"s") %
"rate";
-
997 for (
auto const& proc : procs_) {
-
998 txt_file << format(
"%-"+getProcLen(proc)+
".6g ") % proc->
no_norm_rate();
-
-
-
1001 txt_file << dashes <<
"\n";
-
-
-
-
1005 for (
auto par : params_) {
-
-
1007 if (p->
err_d() != 0.0 && p->
err_u() != 0.0 &&
-
1008 all_dependents_pars.count(p->
name()) && sys_set.count(p->
name())) {
-
1009 txt_file << format((format(
"%%-%is param %%g %%g") % sys_str_len).str()) %
-
-
1011 if (p->
range_d() != std::numeric_limits<double>::lowest() &&
-
1012 p->
range_u() != std::numeric_limits<double>::max()) {
-
1013 txt_file << format(
" [%.4g,%.4g]") % p->
range_d() % p->
range_u();
-
-
-
+
+
979 txt_file << dashes <<
"\n";
+
+
+
982 unsigned sys_str_len = 14;
+
983 for (
auto const& sys : sys_set) {
+
984 if (sys.length() > sys_str_len) sys_str_len = sys.length();
+
+
986 for (
auto const& sys : all_dependents_pars) {
+
987 if (sys.length() > sys_str_len) sys_str_len = sys.length();
+
+
989 std::string sys_str_short = boost::lexical_cast<std::string>(sys_str_len);
+
990 std::string sys_str_long = boost::lexical_cast<std::string>(sys_str_len+9);
+
+
992 auto getProcLen = [](
auto const& proc) -> std::string {
+
993 std::size_t proc_len = 15;
+
994 proc_len = std::max(proc_len, proc->
process().size());
+
995 std::string proc_len_str = boost::lexical_cast<std::string>(proc_len);
+
+
+
+
999 txt_file << format(
"%-"+sys_str_long+
"s") %
"bin";
+
1000 for (
auto const& proc : procs_) {
+
1001 if (proc->
shape()) {
+
1002 bool add_dir = TH1::AddDirectoryStatus();
+
1003 TH1::AddDirectory(
false);
+
+
1005 WriteHistToFile(h.get(), &root_file, mappings, proc->
bin(),
+
+
1007 TH1::AddDirectory(add_dir);
+
+
1009 txt_file << format(
"%-"+getProcLen(proc)+
"s ") % proc->
bin();
+
+
+
+
1013 txt_file << format(
"%-"+sys_str_long+
"s") %
"process";
+
+
1015 for (
auto const& proc : procs_) {
+
1016 txt_file << format(
"%-"+getProcLen(proc)+
"s ") % proc->
process();
-
-
1019 for (
auto const& sys : sys_set) {
-
1020 std::vector<std::string> line(procs_.size() + 2);
-
-
1022 bool seen_lnN =
false;
-
1023 bool seen_lnU =
false;
-
1024 bool seen_shape =
false;
-
1025 bool seen_shapeN2 =
false;
-
1026 bool seen_shapeU =
false;
-
1027 for (
unsigned p = 0; p < procs_.size(); ++p) {
-
-
1029 for (
unsigned n = 0; n < proc_sys_map[p].size(); ++n) {
-
-
-
1032 if (ptr->
name() != sys)
continue;
-
1033 std::string tp = ptr->
type();
-
1034 if (tp ==
"lnN" || tp ==
"lnU") {
-
1035 if (tp ==
"lnN") seen_lnN =
true;
-
1036 if (tp ==
"lnU") seen_lnU =
true;
-
-
-
-
1040 : (format(
"%g") % ptr->
value_u()).str();
-
-
-
1043 if (tp ==
"shape" || tp ==
"shapeN2" || tp ==
"shapeU") {
-
1044 if (tp ==
"shape") seen_shape =
true;
-
1045 if (tp ==
"shapeN2") seen_shapeN2 =
true;
-
1046 if (tp ==
"shapeU") seen_shapeU =
true;
-
1047 line[p + 2] = (format(
"%g") % ptr->
scale()).str();
-
-
1049 bool add_dir = TH1::AddDirectoryStatus();
-
1050 TH1::AddDirectory(
false);
-
-
1052 h_d->Scale(procs_[p]->rate() * ptr->
value_d());
-
1053 WriteHistToFile(h_d.get(), &root_file, mappings, ptr->
bin(),
-
-
-
1056 h_u->Scale(procs_[p]->rate() * ptr->
value_u());
-
1057 WriteHistToFile(h_u.get(), &root_file, mappings, ptr->
bin(),
-
-
1059 TH1::AddDirectory(add_dir);
-
-
-
-
1063 if (!flags_.at(
"allow-missing-shapes")) {
-
1064 std::stringstream err;
-
1065 err <<
"Trying to write shape uncertainty with missing "
-
-
-
1068 throw std::runtime_error(
FNERROR(err.str()));
-
-
-
-
-
-
-
1075 line[1] =
"shapeN2";
-
1076 }
else if (seen_shapeU) {
-
-
1078 }
else if (seen_lnU) {
-
-
1080 }
else if (seen_lnN && !seen_shape) {
-
-
1082 }
else if (!seen_lnN && seen_shape) {
-
-
1084 }
else if (seen_lnN && seen_shape) {
-
-
-
1087 throw std::runtime_error(
FNERROR(
"Systematic type could not be deduced"));
-
-
1089 txt_file << format(
"%-" + sys_str_short +
"s %-7s ") % line[0] % line[1];
-
1090 for (
unsigned p = 0; p < procs_.size(); ++p) {
-
1091 txt_file << format(
"%-"+getProcLen(procs_[p])+
"s ") % line[p + 2];
-
-
-
-
-
-
-
-
-
1100 std::vector<std::vector<std::string> > floating_params;
-
1101 for (
auto const& sys : rateparam_set) {
-
1102 FNLOGC(log(), verbosity_ > 1) <<
"Doing rateParam: " << sys <<
"\n";
-
-
-
1105 FNLOGC(log(), verbosity_ > 1) <<
"Procs for this rateParam: \n";
-
1106 for (
auto const& proc : procs_rp) {
-
1107 FNLOGC(log(), verbosity_ > 1) <<
" - " << proc <<
"\n";
-
-
-
-
-
1112 if (bins_rp.size() > 0 && bins_rp.size() == bins_all.size()) {
-
1113 floating_params.push_back({sys,
"*", proc});
-
-
1115 for (
auto const&
bin : bins_rp) {
-
1116 floating_params.push_back({sys,
bin, proc});
-
-
-
1119 FNLOGC(log(), verbosity_ > 1)
-
1120 << bins_rp.size() <<
"/" << bins_all.size()
-
1121 <<
" bins with this process have a rateParam\n";
-
-
-
1124 std::set<std::string> frozen_params;
-
1125 for (
auto const& rp : floating_params) {
-
1126 if (params_.count(rp[0])) {
-
-
1128 txt_file << format(
"%-" + sys_str_short +
"s %-10s %-10s %-10s %g") %
-
1129 rp[0] %
"rateParam" % rp[1] % rp[2] % par->
val();
-
-
-
1132 if (par->
range_d() != std::numeric_limits<double>::lowest() &&
-
1133 par->
range_u() != std::numeric_limits<double>::max()) {
-
1134 txt_file << format(
" [%.4g,%.4g]") % par->
range_d() % par->
range_u();
-
-
-
-
1138 frozen_params.insert(rp[0]);
-
-
-
-
1142 for (
auto const& par : frozen_params) {
-
1143 txt_file <<
"nuisance edit freeze " << par <<
"\n";
-
-
1145 std::set<std::string> all_fn_param_args;
-
1146 for (
auto const& rp : floating_params) {
-
1147 if (!params_.count(rp[0])) {
-
-
-
1150 RooWorkspace *rp_ws = wspaces_.at(
"_rateParams").get();
-
1151 RooAbsArg *arg = rp_ws->arg(rp[0].c_str());
-
1152 if (arg->getStringAttribute(
"wspSource")) {
-
1153 txt_file << format(
"%-" + sys_str_short +
-
1154 "s %-10s %-10s %-10s %-20s\n") %
-
1155 rp[0] %
"rateParam" % rp[1] % rp[2] % std::string(arg->getStringAttribute(
"wspSource"));
-
-
-
1158 RooFormulaVar* form =
-
1159 dynamic_cast<RooFormulaVar*>(arg);
-
-
-
1162 std::stringstream ss;
-
1163 form->printMetaArgs(ss);
-
1164 std::string form_str = ss.str();
-
-
-
1167 std::size_t first = form_str.find_first_of(
'"');
-
1168 std::size_t last = form_str.find_last_of(
'"');
-
1169 form_str = form_str.substr(first + 1, last - first - 1);
-
-
-
1172 std::unique_ptr<RooArgSet>(form->getParameters(RooArgList()))->getSize();
-
1173 std::string args =
"";
-
1174 for (
unsigned i = 0; i < npars; ++i) {
-
1175 all_fn_param_args.insert(std::string(form->getParameter(i)->GetName()));
-
1176 args += std::string(form->getParameter(i)->GetName());
-
1177 if (i < (npars-1)) {
-
-
-
-
1181 txt_file << format(
"%-" + sys_str_short +
-
1182 "s %-10s %-10s %-10s %-20s %s\n") %
-
1183 rp[0] %
"rateParam" % rp[1] % rp[2] % form_str % args;
-
-
-
-
1187 if (wspaces_.count(
"_rateParams")) {
-
1188 RooWorkspace *rp_ws = wspaces_.at(
"_rateParams").get();
-
1189 RooArgSet vars = rp_ws->allVars();
-
1190 auto v = vars.createIterator();
-
-
1192 RooRealVar *y = dynamic_cast<RooRealVar*>(**v);
-
1193 if (y && y->getAttribute(
"extArg") && all_fn_param_args.count(std::string(y->GetName()))) {
-
1194 Parameter const* p = params_.at(y->GetName()).get();
-
1195 txt_file << format(
"%-" + sys_str_short +
"s %-10s %g") %
-
1196 y->GetName() %
"extArg" % p->
val();
-
1197 if (p->
range_d() != std::numeric_limits<double>::lowest() &&
-
1198 p->
range_u() != std::numeric_limits<double>::max()) {
-
1199 txt_file << format(
" [%.4g,%.4g]") % p->
range_d() % p->
range_u();
-
-
-
-
-
1204 }
while (v->Next());
-
1205 RooArgSet funcs = rp_ws->allFunctions();
-
1206 v = funcs.createIterator();
-
-
1208 RooAbsReal *y = dynamic_cast<RooAbsReal*>(**v);
-
1209 if (y && y->getAttribute(
"extArg") && y->getStringAttribute(
"wspSource") && all_fn_param_args.count(std::string(y->GetName()))) {
-
1210 txt_file << format(
"%-" + sys_str_short +
-
1211 "s %-10s %-20s\n") %
-
1212 y->GetName() %
"extArg" % std::string(y->getStringAttribute(
"wspSource"));
-
-
-
-
1216 }
while (v->Next());
-
-
-
1219 std::set<std::string> ws_vars;
-
1220 for (
auto iter : wspaces_) {
-
1221 RooArgSet vars = iter.second->allVars();
-
1222 auto v = vars.createIterator();
-
-
1224 RooRealVar *y = dynamic_cast<RooRealVar*>(**v);
-
-
1226 ws_vars.insert(y->GetName());
-
-
1228 }
while (v->Next());
-
-
-
-
-
-
-
-
-
-
1238 for (
auto par : params_) {
-
-
1240 if (p->
err_d() != 0.0 && p->
err_u() != 0.0 &&
-
1241 all_dependents_pars.count(p->
name()) && !sys_set.count(p->
name())) {
-
1242 txt_file << format((format(
"%%-%is param %%g %%g") % sys_str_len).str()) %
-
-
1244 if (p->
range_d() != std::numeric_limits<double>::lowest() &&
-
1245 p->
range_u() != std::numeric_limits<double>::max()) {
-
1246 txt_file << format(
" [%.4g,%.4g]") % p->
range_d() % p->
range_u();
-
-
-
-
-
-
1252 for (
auto par : multipdf_cats) {
-
1253 txt_file << format((format(
"%%-%is discrete\n") % sys_str_len).str()) % par;
-
-
-
1256 for (
auto stat_settings : auto_stats_settings_) {
-
1257 if (!
bin_set.count(stat_settings.first))
continue;
-
1258 txt_file << format(
"%s autoMCStats %g %i %i\n") % stat_settings.first %
-
1259 stat_settings.second.event_threshold %
-
1260 stat_settings.second.include_signal %
-
1261 stat_settings.second.hist_mode;
-
-
-
-
-
-
-
1268 std::map<std::string, std::set<std::string>> groups;
-
1269 for (
auto const& par : params_) {
-
-
1271 if (!all_dependents_pars.count(p->
name()) && !sys_set.count(p->
name())) {
-
-
-
1274 if (p->
err_d() == 0.0 && p->
err_u() == 0.0)
continue;
-
1275 for (
auto const& str : p->
groups()) {
-
1276 if (!groups.count(str)) {
-
1277 groups[str] = std::set<std::string>();
-
-
1279 groups[str].insert(p->
name());
-
-
-
1282 for (
auto const& gr : groups) {
-
1283 txt_file << format(
"%s group =") % gr.first;
-
1284 for (
auto const& np : gr.second) {
-
1285 txt_file <<
" " << np;
-
-
-
-
-
1290 for (
auto const& postl : post_lines_) {
-
1291 txt_file << postl <<
"\n";
-
-
-
-
1295 void CombineHarvester::WriteHistToFile(
-
-
-
1298 std::vector<HistMapping>
const& mappings,
-
1299 std::string
const& bin,
-
-
1301 std::string
const& mass,
-
1302 std::string
const& nuisance,
-
-
1304 StrPairVec attempts = this->GenerateShapeMapAttempts(
process,
bin);
-
1305 for (
unsigned a = 0; a < attempts.size(); ++a) {
-
1306 for (
unsigned m = 0; m < mappings.size(); ++m) {
-
1307 if ((attempts[a].first == mappings[m].
process) &&
-
1308 (attempts[a].second == mappings[m].category)) {
-
1309 std::string p = (type == 0 ?
-
1310 mappings[m].pattern : mappings[m].syst_pattern);
-
1311 boost::replace_all(p,
"$CHANNEL",
bin);
-
1312 boost::replace_all(p,
"$PROCESS",
process);
-
1313 boost::replace_all(p,
"$MASS",
mass);
-
1314 if (type == 1) boost::replace_all(p,
"$SYSTEMATIC", nuisance+
"Down");
-
1315 if (type == 2) boost::replace_all(p,
"$SYSTEMATIC", nuisance+
"Up");
-
-
-
-
-
-
-
-
@@ -189,7 +186,7 @@
81 std::vector<RooRealVar const*> r_vec(n_pars,
nullptr);
82 std::vector<ch::Parameter*> p_vec(n_pars,
nullptr);
83 for (
unsigned n = 0; n < p_vec.size(); ++n) {
-
84 r_vec[n] = dynamic_cast<RooRealVar const*>(rands.at(n));
+
84 r_vec[n] =
dynamic_cast<RooRealVar const*
>(rands.at(n));
@@ -242,740 +239,750 @@
134 unsigned n_samples) {
135 auto lookup = GenerateProcSystMap();
136 TH1F shape = GetShapeInternal(lookup);
-
137 for (
int i = 1; i <= shape.GetNbinsX(); ++i) {
-
138 shape.SetBinError(i, 0.0);
-
-
-
-
-
-
-
145 RooArgList
const& rands = fit.randomizePars();
-
-
-
-
149 int n_pars = rands.getSize();
-
150 std::vector<RooRealVar const*> r_vec(n_pars,
nullptr);
-
151 std::vector<ch::Parameter*> p_vec(n_pars,
nullptr);
-
152 for (
unsigned n = 0; n < p_vec.size(); ++n) {
-
153 r_vec[n] = dynamic_cast<RooRealVar const*>(rands.at(n));
-
-
-
-
-
158 for (
unsigned i = 0; i < n_samples; ++i) {
-
-
-
161 for (
int n = 0; n < n_pars; ++n) {
-
162 if (p_vec[n]) p_vec[n]->set_val(r_vec[n]->getVal());
-
-
-
165 TH1F rand_shape = this->GetShapeInternal(lookup);
-
166 for (
int i = 1; i <= shape.GetNbinsX(); ++i) {
-
-
168 std::fabs(rand_shape.GetBinContent(i) - shape.GetBinContent(i));
-
169 shape.SetBinError(i, err*err + shape.GetBinError(i));
-
-
-
172 for (
int i = 1; i <= shape.GetNbinsX(); ++i) {
-
173 shape.SetBinError(i, std::sqrt(shape.GetBinError(i)/double(n_samples)));
-
-
-
-
-
-
-
180 unsigned n_samples) {
-
181 auto lookup = GenerateProcSystMap();
-
-
183 unsigned n = procs_.size();
-
184 TH1F nom(
"nominal",
"nominal", n, 0, n);
-
185 TH2F res(
"covariance",
"covariance", n, 0, n, n, 0, n);
-
-
187 std::vector<CombineHarvester> ch_procs;
-
188 std::vector<std::string> labels;
-
189 unsigned nbins = this->
bin_set().size();
-
-
191 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
-
193 this->
cp().
bin({procs_[i]->bin()}).
process({procs_[i]->process()}));
-
-
195 labels.push_back(procs_[i]->
bin() +
"," + procs_[i]->
process());
-
-
197 labels.push_back(procs_[i]->
process());
-
-
-
200 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
201 nom.SetBinContent(i + 1, ch_procs[i].
GetRate());
-
-
-
-
-
-
207 RooArgList
const& rands = fit.randomizePars();
-
-
-
-
211 int n_pars = rands.getSize();
-
212 std::vector<RooRealVar const*> r_vec(n_pars,
nullptr);
-
213 std::vector<ch::Parameter*> p_vec(n_pars,
nullptr);
-
214 for (
unsigned n = 0; n < p_vec.size(); ++n) {
-
215 r_vec[n] = dynamic_cast<RooRealVar const*>(rands.at(n));
-
-
+
+
138 shape.Copy(yield_sum);
+
139 int n_bins = shape.GetNbinsX();
+
140 std::vector<std::vector<float>> sample_yields_perbin(n_bins, std::vector<float>(n_samples));
+
+
142 for (
int i = 1; i <= n_bins; ++i) {
+
143 shape.SetBinError(i, 0.0);
+
144 yield_sum.SetBinContent(i, 0.0);
+
+
+
+
+
+
+
151 RooArgList
const& rands = fit.randomizePars();
+
+
+
+
155 int n_pars = rands.getSize();
+
156 std::vector<RooRealVar const*> r_vec(n_pars,
nullptr);
+
157 std::vector<ch::Parameter*> p_vec(n_pars,
nullptr);
+
158 for (
unsigned n = 0; n < p_vec.size(); ++n) {
+
159 r_vec[n] =
dynamic_cast<RooRealVar const*
>(rands.at(n));
+
+
+
+
+
164 for (
unsigned i = 0; i < n_samples; ++i) {
+
+
+
167 for (
int n = 0; n < n_pars; ++n) {
+
168 if (p_vec[n]) p_vec[n]->set_val(r_vec[n]->getVal());
+
+
170 TH1F rand_shape = this->GetShapeInternal(lookup);
+
171 for (
int j = 1; j <= n_bins; ++j) {
+
172 double yield = rand_shape.GetBinContent(j);
+
173 sample_yields_perbin[j-1][i] = yield;
+
174 yield_sum.SetBinContent(j, yield_sum.GetBinContent(j) + yield);
+
+
+
177 for (
int i = 1; i <= n_bins; ++i) {
+
178 double yield_mean = yield_sum.GetBinContent(i)/double(n_samples);
+
179 for (
auto x :sample_yields_perbin[i-1]){
+
180 double err = std::fabs(x - yield_mean);
+
181 shape.SetBinError(i, err*err + shape.GetBinError(i));
+
+
183 shape.SetBinError(i, std::sqrt(shape.GetBinError(i)/
double(n_samples)));
+
+
+
+
+
+
+
190 unsigned n_samples) {
+
191 auto lookup = GenerateProcSystMap();
+
+
193 unsigned n = procs_.size();
+
194 TH1F nom(
"nominal",
"nominal", n, 0, n);
+
195 TH2F res(
"covariance",
"covariance", n, 0, n, n, 0, n);
+
+
197 std::vector<CombineHarvester> ch_procs;
+
198 std::vector<std::string> labels;
+
199 unsigned nbins = this->
bin_set().size();
+
+
201 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
+
203 this->
cp().
bin({procs_[i]->bin()}).
process({procs_[i]->process()}));
+
+
205 labels.push_back(procs_[i]->
bin() +
"," + procs_[i]->
process());
+
+
207 labels.push_back(procs_[i]->
process());
+
+
+
210 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
211 nom.SetBinContent(i + 1, ch_procs[i].
GetRate());
+
+
+
+
+
+
217 RooArgList
const& rands = fit.randomizePars();
-
-
220 for (
unsigned rnd = 0; rnd < n_samples; ++rnd) {
-
-
-
223 for (
int n = 0; n < n_pars; ++n) {
-
224 if (p_vec[n]) p_vec[n]->set_val(r_vec[n]->getVal());
-
-
-
227 for (
int i = 1; i <= nom.GetNbinsX(); ++i) {
-
228 for (
int j = 1; j <= nom.GetNbinsX(); ++j) {
-
-
230 int y = nom.GetNbinsX() - (i - 1);
-
231 res.GetXaxis()->SetBinLabel(x, labels[j - 1].c_str());
-
232 res.GetYaxis()->SetBinLabel(y, labels[i - 1].c_str());
-
-
234 x, y, res.GetBinContent(x, y) +
-
235 (ch_procs[i - 1].GetRate() - nom.GetBinContent(i)) *
-
236 (ch_procs[j - 1].
GetRate() - nom.GetBinContent(j)));
-
-
-
-
-
241 for (
int i = 1; i <= nom.GetNbinsX(); ++i) {
-
242 for (
int j = 1; j <= nom.GetNbinsX(); ++j) {
-
-
244 int y = nom.GetNbinsX() - (i - 1);
-
245 res.SetBinContent(x, y, res.GetBinContent(x, y) / double(n_samples));
-
-
-
-
-
-
-
-
253 unsigned n_samples) {
-
-
-
256 res.SetName(
"correlation");
-
257 res.SetTitle(
"correlation");
-
258 for (
int i = 1; i <= cov.GetNbinsX(); ++i) {
-
259 for (
int j = 1; j <= cov.GetNbinsX(); ++j) {
-
-
261 int y = cov.GetNbinsX() - (i - 1);
-
-
-
264 cov.GetBinContent(x, y) /
-
265 (std::sqrt(cov.GetBinContent(i, cov.GetNbinsX() - (i - 1))) *
-
266 std::sqrt(cov.GetBinContent(j, cov.GetNbinsX() - (j - 1)))));
-
-
-
-
-
-
-
273 auto lookup = GenerateProcSystMap();
-
274 return GetRateInternal(lookup);
-
-
-
-
278 auto lookup = GenerateProcSystMap();
-
279 return GetShapeInternal(lookup);
+
+
+
221 int n_pars = rands.getSize();
+
222 std::vector<RooRealVar const*> r_vec(n_pars,
nullptr);
+
223 std::vector<ch::Parameter*> p_vec(n_pars,
nullptr);
+
224 for (
unsigned n = 0; n < p_vec.size(); ++n) {
+
225 r_vec[n] =
dynamic_cast<RooRealVar const*
>(rands.at(n));
+
+
+
+
+
230 for (
unsigned rnd = 0; rnd < n_samples; ++rnd) {
+
+
+
233 for (
int n = 0; n < n_pars; ++n) {
+
234 if (p_vec[n]) p_vec[n]->set_val(r_vec[n]->getVal());
+
+
+
237 for (
int i = 1; i <= nom.GetNbinsX(); ++i) {
+
238 for (
int j = 1; j <= nom.GetNbinsX(); ++j) {
+
+
240 int y = nom.GetNbinsX() - (i - 1);
+
241 res.GetXaxis()->SetBinLabel(x, labels[j - 1].c_str());
+
242 res.GetYaxis()->SetBinLabel(y, labels[i - 1].c_str());
+
+
244 x, y, res.GetBinContent(x, y) +
+
245 (ch_procs[i - 1].GetRate() - nom.GetBinContent(i)) *
+
246 (ch_procs[j - 1].
GetRate() - nom.GetBinContent(j)));
+
+
+
+
+
251 for (
int i = 1; i <= nom.GetNbinsX(); ++i) {
+
252 for (
int j = 1; j <= nom.GetNbinsX(); ++j) {
+
+
254 int y = nom.GetNbinsX() - (i - 1);
+
255 res.SetBinContent(x, y, res.GetBinContent(x, y) / double(n_samples));
+
+
+
+
+
+
+
+
263 unsigned n_samples) {
+
+
+
266 res.SetName(
"correlation");
+
267 res.SetTitle(
"correlation");
+
268 for (
int i = 1; i <= cov.GetNbinsX(); ++i) {
+
269 for (
int j = 1; j <= cov.GetNbinsX(); ++j) {
+
+
271 int y = cov.GetNbinsX() - (i - 1);
+
+
+
274 cov.GetBinContent(x, y) /
+
275 (std::sqrt(cov.GetBinContent(i, cov.GetNbinsX() - (i - 1))) *
+
276 std::sqrt(cov.GetBinContent(j, cov.GetNbinsX() - (j - 1)))));
+
+
+
-
282 double CombineHarvester::GetRateInternal(ProcSystMap
const& lookup,
-
283 std::string
const& single_sys) {
-
-
-
286 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
287 double p_rate = procs_[i]->rate();
-
-
-
-
-
-
293 if (single_sys !=
"" && !procs_[i]->pdf()) {
-
-
295 return sys->
name() == single_sys;
-
-
-
298 for (
auto sys_it : lookup[i]) {
-
299 if (sys_it->type() ==
"rateParam") {
-
-
-
302 double x = params_[sys_it->name()]->val();
-
303 if (sys_it->asymm()) {
-
304 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
-
-
-
307 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
-
-
-
-
-
-
-
-
315 TH1F CombineHarvester::GetShapeInternal(ProcSystMap
const& lookup,
-
316 std::string
const& single_sys) {
-
-
318 bool shape_init =
false;
-
-
320 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
-
-
-
324 if (single_sys !=
"" && !procs_[i]->pdf()) {
-
325 if (!
ch::any_of(lookup[i], [&](Systematic
const* sys) {
-
326 return sys->name() == single_sys;
-
-
+
+
283 auto lookup = GenerateProcSystMap();
+
284 return GetRateInternal(lookup);
+
+
+
+
288 auto lookup = GenerateProcSystMap();
+
289 return GetShapeInternal(lookup);
+
+
+
292 double CombineHarvester::GetRateInternal(ProcSystMap
const& lookup,
+
293 std::string
const& single_sys) {
+
+
+
296 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
297 double p_rate = procs_[i]->rate();
+
+
+
+
+
+
303 if (single_sys !=
"" && !procs_[i]->pdf()) {
+
+
305 return sys->
name() == single_sys;
+
+
+
308 for (
auto sys_it : lookup[i]) {
+
309 if (sys_it->type() ==
"rateParam") {
+
+
+
312 double x = params_[sys_it->name()]->val();
+
313 if (sys_it->asymm()) {
+
314 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
+
+
+
317 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
+
+
+
+
+
+
+
+
325 TH1F CombineHarvester::GetShapeInternal(ProcSystMap
const& lookup,
+
326 std::string
const& single_sys) {
+
+
328 bool shape_init =
false;
-
330 double p_rate = procs_[i]->rate();
-
331 if (procs_[i]->shape() || procs_[i]->
data()) {
-
332 TH1F proc_shape = procs_[i]->ShapeAsTH1F();
-
333 for (
auto sys_it : lookup[i]) {
-
334 if (sys_it->type() ==
"rateParam") {
-
-
-
337 auto param_it = params_.find(sys_it->name());
-
338 if (param_it == params_.end()) {
-
339 throw std::runtime_error(
-
340 FNERROR(
"Parameter " + sys_it->name() +
-
341 " not found in CombineHarvester instance"));
-
-
343 double x = param_it->second->val();
-
344 if (sys_it->asymm()) {
-
345 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
-
-
347 if (sys_it->type() ==
"shape" || sys_it->type() ==
"shapeN2" ||
-
348 sys_it->type() ==
"shapeU") {
-
-
350 if (sys_it->type() ==
"shapeN2") linear =
false;
-
351 if (sys_it->shape_u() && sys_it->shape_d()) {
-
352 ShapeDiff(x * sys_it->scale(), &proc_shape, procs_[i]->shape(),
-
353 sys_it->shape_d(), sys_it->shape_u(), linear);
-
-
355 if (sys_it->data_u() && sys_it->data_d()) {
-
356 RooDataHist
const* nom =
-
357 dynamic_cast<RooDataHist const*>(procs_[i]->
data());
-
-
359 ShapeDiff(x * sys_it->scale(), &proc_shape, nom,
-
360 sys_it->data_d(), sys_it->data_u());
-
-
-
-
-
365 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
-
-
-
368 for (
int b = 1; b <= proc_shape.GetNbinsX(); ++b) {
-
369 if (proc_shape.GetBinContent(b) < 0.) proc_shape.SetBinContent(b, 0.);
-
-
371 proc_shape.Scale(p_rate);
-
-
373 proc_shape.Copy(shape);
-
-
-
-
377 shape.Add(&proc_shape);
-
378 }
else if (procs_[i]->pdf()) {
-
379 if (!procs_[i]->observable()) {
-
380 RooAbsData
const* data_obj = FindMatchingData(procs_[i].get());
-
381 std::string var_name =
"CMS_th1x";
-
382 if (data_obj) var_name = data_obj->get()->first()->GetName();
-
383 procs_[i]->set_observable((RooRealVar *)procs_[i]->pdf()->findServer(var_name.c_str()));
-
-
385 TH1::AddDirectory(
false);
-
386 TH1F* tmp = (TH1F*)procs_[i]->observable()->createHistogram(
"");
-
387 for (
int b = 1; b <= tmp->GetNbinsX(); ++b) {
-
388 procs_[i]->observable()->setVal(tmp->GetBinCenter(b));
-
389 tmp->SetBinContent(b, tmp->GetBinWidth(b) * procs_[i]->pdf()->getVal());
-
-
391 TH1F proc_shape = *tmp;
-
-
393 RooAbsPdf
const* aspdf = dynamic_cast<RooAbsPdf const*>(procs_[i]->pdf());
-
394 if ((aspdf && !aspdf->selfNormalized()) || (!aspdf)) {
-
-
-
397 if (proc_shape.Integral() > 0.) {
-
398 proc_shape.Scale(1. / proc_shape.Integral());
-
+
330 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
+
+
+
334 if (single_sys !=
"" && !procs_[i]->pdf()) {
+
335 if (!
ch::any_of(lookup[i], [&](Systematic
const* sys) {
+
336 return sys->name() == single_sys;
+
+
+
+
340 double p_rate = procs_[i]->rate();
+
341 if (procs_[i]->shape() || procs_[i]->
data()) {
+
342 TH1F proc_shape = procs_[i]->ShapeAsTH1F();
+
343 for (
auto sys_it : lookup[i]) {
+
344 if (sys_it->type() ==
"rateParam") {
+
+
+
347 auto param_it = params_.find(sys_it->name());
+
348 if (param_it == params_.end()) {
+
349 throw std::runtime_error(
+
350 FNERROR(
"Parameter " + sys_it->name() +
+
351 " not found in CombineHarvester instance"));
+
+
353 double x = param_it->second->val();
+
354 if (sys_it->asymm()) {
+
355 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
+
+
357 if (sys_it->type() ==
"shape" || sys_it->type() ==
"shapeN2" ||
+
358 sys_it->type() ==
"shapeU") {
+
+
360 if (sys_it->type() ==
"shapeN2") linear =
false;
+
361 if (sys_it->shape_u() && sys_it->shape_d()) {
+
362 ShapeDiff(x * sys_it->scale(), &proc_shape, procs_[i]->shape(),
+
363 sys_it->shape_d(), sys_it->shape_u(), linear);
+
+
365 if (sys_it->data_u() && sys_it->data_d()) {
+
366 RooDataHist
const* nom =
+
367 dynamic_cast<RooDataHist const*
>(procs_[i]->data());
+
+
369 ShapeDiff(x * sys_it->scale(), &proc_shape, nom,
+
370 sys_it->data_d(), sys_it->data_u());
+
+
+
+
+
375 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
+
+
+
378 for (
int b = 1; b <= proc_shape.GetNbinsX(); ++b) {
+
379 if (proc_shape.GetBinContent(b) < 0.) proc_shape.SetBinContent(b, 0.);
+
+
381 proc_shape.Scale(p_rate);
+
+
383 proc_shape.Copy(shape);
+
+
+
+
387 shape.Add(&proc_shape);
+
388 }
else if (procs_[i]->pdf()) {
+
389 if (!procs_[i]->observable()) {
+
390 RooAbsData
const* data_obj = FindMatchingData(procs_[i].get());
+
391 std::string var_name =
"CMS_th1x";
+
392 if (data_obj) var_name = data_obj->get()->first()->GetName();
+
393 procs_[i]->set_observable((RooRealVar *)procs_[i]->pdf()->findServer(var_name.c_str()));
+
+
395 TH1::AddDirectory(
false);
+
396 TH1F* tmp = (TH1F*)procs_[i]->observable()->createHistogram(
"");
+
397 for (
int b = 1; b <= tmp->GetNbinsX(); ++b) {
+
398 procs_[i]->observable()->setVal(tmp->GetBinCenter(b));
+
399 tmp->SetBinContent(b, tmp->GetBinWidth(b) * procs_[i]->pdf()->getVal());
-
401 for (
auto sys_it : lookup[i]) {
-
402 if (sys_it->type() ==
"rateParam") {
-
-
-
405 double x = params_[sys_it->name()]->val();
-
406 if (sys_it->asymm()) {
-
407 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
-
-
-
410 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
-
-
-
413 proc_shape.Scale(p_rate);
-
-
415 proc_shape.Copy(shape);
-
-
-
-
419 shape.Add(&proc_shape);
-
-
-
-
-
-
-
-
427 for (
unsigned i = 0; i < obs_.size(); ++i) {
-
428 rate += obs_[i]->rate();
-
-
-
-
-
-
-
435 bool shape_init =
false;
-
+
401 TH1F proc_shape = *tmp;
+
+
403 RooAbsPdf
const* aspdf =
dynamic_cast<RooAbsPdf const*
>(procs_[i]->pdf());
+
404 if ((aspdf && !aspdf->selfNormalized()) || (!aspdf)) {
+
+
+
407 if (proc_shape.Integral() > 0.) {
+
408 proc_shape.Scale(1. / proc_shape.Integral());
+
+
+
411 for (
auto sys_it : lookup[i]) {
+
412 if (sys_it->type() ==
"rateParam") {
+
+
+
415 double x = params_[sys_it->name()]->val();
+
416 if (sys_it->asymm()) {
+
417 p_rate *= logKappaForX(x * sys_it->scale(), sys_it->value_d(),
+
+
+
420 p_rate *= std::pow(sys_it->value_u(), x * sys_it->scale());
+
+
+
423 proc_shape.Scale(p_rate);
+
+
425 proc_shape.Copy(shape);
+
+
+
+
429 shape.Add(&proc_shape);
+
+
+
+
+
+
+
437 for (
unsigned i = 0; i < obs_.size(); ++i) {
-
-
439 double p_rate = obs_[i]->rate();
-
440 if (obs_[i]->shape()) {
-
441 proc_shape = obs_[i]->ShapeAsTH1F();
-
442 }
else if (obs_[i]->
data()) {
-
443 TH1F* tmp = dynamic_cast<TH1F*>(obs_[i]->
data()->createHistogram(
-
444 "", *(RooRealVar*)obs_[i]->
data()->get()->first()));
-
-
446 tmp->SetBinErrorOption(TH1::kPoisson);
-
-
-
449 proc_shape.Scale(1. / proc_shape.Integral());
-
-
451 proc_shape.Scale(p_rate);
-
-
453 proc_shape.Copy(shape);
-
-
-
-
457 shape.Add(&proc_shape);
-
-
-
-
-
462 void CombineHarvester::ShapeDiff(
double x,
-
-
-
-
-
-
468 double fx = smoothStepFunc(x);
-
469 for (
int i = 1; i <= target->GetNbinsX(); ++i) {
-
470 float h = high->GetBinContent(i);
-
471 float l = low->GetBinContent(i);
-
472 float n = nom->GetBinContent(i);
-
-
474 float t = target->GetBinContent(i);
-
475 target->SetBinContent(i, t > 0. ? std::log(t) : -999.);
-
476 h = (h > 0. && n > 0.) ? std::log(h/n) : 0.;
-
477 l = (l > 0. && n > 0.) ? std::log(l/n) : 0.;
-
478 target->SetBinContent(i, target->GetBinContent(i) +
-
479 0.5 * x * ((h - l) + (h + l) * fx));
-
480 target->SetBinContent(i, std::exp(target->GetBinContent(i)));
-
-
482 target->SetBinContent(i, target->GetBinContent(i) +
-
483 0.5 * x * ((h - l) + (h + l - 2. * n) * fx));
-
-
-
-
-
488 void CombineHarvester::ShapeDiff(
double x,
-
-
490 RooDataHist
const* nom,
-
491 RooDataHist
const* low,
-
492 RooDataHist
const* high) {
-
493 double fx = smoothStepFunc(x);
-
494 for (
int i = 1; i <= target->GetNbinsX(); ++i) {
-
-
-
-
-
-
500 float h = high->weight() / high->sumEntries();
-
501 float l = low->weight() / low->sumEntries();
-
502 float n = nom->weight() / nom->sumEntries();
-
503 target->SetBinContent(i, target->GetBinContent(i) +
-
504 0.5 * x * ((h - l) + (h + l - 2. * n) * fx));
-
-
-
-
-
-
-
-
-
-
-
-
516 std::string
const& newname) {
-
517 auto it = params_.find(oldname);
-
518 if (it != params_.end()) {
-
519 params_[newname] = it->second;
-
520 params_[newname]->set_name(newname);
-
-
-
+
438 rate += obs_[i]->rate();
+
+
+
+
+
+
+
445 bool shape_init =
false;
+
+
447 for (
unsigned i = 0; i < obs_.size(); ++i) {
+
+
449 double p_rate = obs_[i]->rate();
+
450 if (obs_[i]->shape()) {
+
451 proc_shape = obs_[i]->ShapeAsTH1F();
+
452 }
else if (obs_[i]->
data()) {
+
453 TH1F* tmp =
dynamic_cast<TH1F*
>(obs_[i]->data()->createHistogram(
+
454 "", *(RooRealVar*)obs_[i]->
data()->get()->first()));
+
+
456 tmp->SetBinErrorOption(TH1::kPoisson);
+
+
+
459 proc_shape.Scale(1. / proc_shape.Integral());
+
+
461 proc_shape.Scale(p_rate);
+
+
463 proc_shape.Copy(shape);
+
+
+
+
467 shape.Add(&proc_shape);
+
+
+
+
+
472 void CombineHarvester::ShapeDiff(
double x,
+
+
+
+
+
+
478 double fx = smoothStepFunc(x);
+
479 for (
int i = 1; i <= target->GetNbinsX(); ++i) {
+
480 float h = high->GetBinContent(i);
+
481 float l = low->GetBinContent(i);
+
482 float n = nom->GetBinContent(i);
+
+
484 float t = target->GetBinContent(i);
+
485 target->SetBinContent(i, t > 0. ? std::log(t) : -999.);
+
486 h = (h > 0. && n > 0.) ? std::log(h/n) : 0.;
+
487 l = (l > 0. && n > 0.) ? std::log(l/n) : 0.;
+
488 target->SetBinContent(i, target->GetBinContent(i) +
+
489 0.5 * x * ((h - l) + (h + l) * fx));
+
490 target->SetBinContent(i, std::exp(target->GetBinContent(i)));
+
+
492 target->SetBinContent(i, target->GetBinContent(i) +
+
493 0.5 * x * ((h - l) + (h + l - 2. * n) * fx));
+
+
+
+
+
498 void CombineHarvester::ShapeDiff(
double x,
+
+
500 RooDataHist
const* nom,
+
501 RooDataHist
const* low,
+
502 RooDataHist
const* high) {
+
503 double fx = smoothStepFunc(x);
+
504 for (
int i = 1; i <= target->GetNbinsX(); ++i) {
+
+
+
+
+
+
510 float h = high->weight() / high->sumEntries();
+
511 float l = low->weight() / low->sumEntries();
+
512 float n = nom->weight() / nom->sumEntries();
+
513 target->SetBinContent(i, target->GetBinContent(i) +
+
514 0.5 * x * ((h - l) + (h + l - 2. * n) * fx));
+
+
+
+
+
+
+
+
+
-
-
526 std::string
const& name)
const {
-
527 auto it = params_.find(name);
+
+
526 std::string
const& newname) {
+
527 auto it = params_.find(oldname);
528 if (it != params_.end()) {
-
529 return it->second.get();
-
-
+
529 params_[newname] = it->second;
+
530 params_[newname]->set_name(newname);
+
-
-
536 auto it = params_.find(name);
-
537 if (it != params_.end()) {
-
538 return it->second.get();
-
-
-
-
-
-
-
545 std::vector<ch::Parameter>
const& params) {
-
546 for (
unsigned i = 0; i < params.size(); ++i) {
-
547 auto it = params_.find(params[i].name());
-
548 if (it != params_.end()) {
-
549 it->second->set_val(params[i].val());
-
550 it->second->set_err_d(params[i].err_d());
-
551 it->second->set_err_u(params[i].err_u());
-
-
553 if (verbosity_ >= 1) {
-
554 LOGLINE(log(),
"Parameter " + params[i].name() +
" is not defined");
-
-
-
-
-
-
-
561 for (
int i = 0; i < fit.floatParsFinal().getSize(); ++i) {
-
562 RooRealVar
const* var =
-
563 dynamic_cast<RooRealVar const*>(fit.floatParsFinal().at(i));
-
-
565 auto it = params_.find(std::string(var->GetName()));
-
566 if (it != params_.end()) {
-
567 it->second->set_val(var->getVal());
-
568 it->second->set_err_d(var->getErrorLo());
-
569 it->second->set_err_u(var->getErrorHi());
-
-
571 if (verbosity_ >= 1) {
-
-
573 "Parameter " + std::string(var->GetName()) +
" is not defined");
-
-
-
-
-
-
-
-
-
-
-
584 std::vector<ch::Parameter> params;
-
585 for (
auto const& it : params_) {
-
586 params.push_back(*(it.second));
-
-
-
-
-
-
-
-
-
-
-
-
598 std::vector<double> prev_proc_rates(procs_.size());
-
-
-
-
602 std::vector<std::unique_ptr<TH1>> scaled_procs(procs_.size());
-
-
604 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
605 if (procs_[i]->shape()) {
-
-
607 std::unique_ptr<TH1> copy(procs_[i]->ClonedScaledShape());
-
-
609 prev_proc_rates[i] = procs_[i]->no_norm_rate();
-
610 std::unique_ptr<TH1> copy2(copy->Rebin(bins.size()-1,
"", &(bins[0])));
-
-
612 procs_[i]->set_shape(std::move(copy2),
true);
-
613 scaled_procs[i] = procs_[i]->ClonedScaledShape();
-
-
-
616 for (
unsigned i = 0; i < obs_.size(); ++i) {
-
617 if (obs_[i]->shape()) {
-
618 std::unique_ptr<TH1> copy(obs_[i]->ClonedScaledShape());
-
619 std::unique_ptr<TH1> copy2(copy->Rebin(bins.size()-1,
"", &(bins[0])));
-
620 obs_[i]->set_shape(std::move(copy2),
true);
-
-
-
623 for (
unsigned i = 0; i < systs_.size(); ++i) {
-
624 TH1
const* proc_hist =
nullptr;
-
625 double prev_rate = 0.;
-
626 for (
unsigned j = 0; j < procs_.size(); ++j) {
-
-
628 proc_hist = scaled_procs[j].get();
-
629 prev_rate = prev_proc_rates[j];
-
+
+
536 std::string
const& name)
const {
+
537 auto it = params_.find(name);
+
538 if (it != params_.end()) {
+
539 return it->second.get();
+
+
+
+
+
+
+
546 auto it = params_.find(name);
+
547 if (it != params_.end()) {
+
548 return it->second.get();
+
+
+
+
+
+
+
555 std::vector<ch::Parameter>
const& params) {
+
556 for (
unsigned i = 0; i < params.size(); ++i) {
+
557 auto it = params_.find(params[i].name());
+
558 if (it != params_.end()) {
+
559 it->second->set_val(params[i].val());
+
560 it->second->set_err_d(params[i].err_d());
+
561 it->second->set_err_u(params[i].err_u());
+
+
563 if (verbosity_ >= 1) {
+
564 LOGLINE(log(),
"Parameter " + params[i].name() +
" is not defined");
+
+
+
+
+
+
+
571 for (
int i = 0; i < fit.floatParsFinal().getSize(); ++i) {
+
572 RooRealVar
const* var =
+
573 dynamic_cast<RooRealVar const*
>(fit.floatParsFinal().at(i));
+
+
575 auto it = params_.find(std::string(var->GetName()));
+
576 if (it != params_.end()) {
+
577 it->second->set_val(var->getVal());
+
578 it->second->set_err_d(var->getErrorLo());
+
579 it->second->set_err_u(var->getErrorHi());
+
+
581 if (verbosity_ >= 1) {
+
+
583 "Parameter " + std::string(var->GetName()) +
" is not defined");
+
+
+
+
+
+
+
+
+
+
+
594 std::vector<ch::Parameter> params;
+
595 for (
auto const& it : params_) {
+
596 params.push_back(*(it.second));
+
+
+
+
+
+
+
+
+
+
+
+
608 std::vector<double> prev_proc_rates(procs_.size());
+
+
+
+
612 std::vector<std::unique_ptr<TH1>> scaled_procs(procs_.size());
+
+
614 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
615 if (procs_[i]->shape()) {
+
+
617 std::unique_ptr<TH1> copy(procs_[i]->ClonedScaledShape());
+
+
619 prev_proc_rates[i] = procs_[i]->no_norm_rate();
+
620 std::unique_ptr<TH1> copy2(copy->Rebin(bins.size()-1,
"", &(bins[0])));
+
+
622 procs_[i]->set_shape(std::move(copy2),
true);
+
623 scaled_procs[i] = procs_[i]->ClonedScaledShape();
+
+
+
626 for (
unsigned i = 0; i < obs_.size(); ++i) {
+
627 if (obs_[i]->shape()) {
+
628 std::unique_ptr<TH1> copy(obs_[i]->ClonedScaledShape());
+
629 std::unique_ptr<TH1> copy2(copy->Rebin(bins.size()-1,
"", &(bins[0])));
+
630 obs_[i]->set_shape(std::move(copy2),
true);
-
632 if (systs_[i]->shape_u() && systs_[i]->shape_d()) {
-
-
634 std::unique_ptr<TH1> copy_u(systs_[i]->ClonedShapeU());
-
635 std::unique_ptr<TH1> copy_d(systs_[i]->ClonedShapeD());
-
-
-
-
-
640 copy_u->Scale(systs_[i]->value_u() * prev_rate);
-
641 copy_d->Scale(systs_[i]->value_d() * prev_rate);
-
-
643 std::unique_ptr<TH1> copy2_u(
-
644 copy_u->Rebin(bins.size() - 1,
"", &(bins[0])));
-
645 std::unique_ptr<TH1> copy2_d(
-
646 copy_d->Rebin(bins.size() - 1,
"", &(bins[0])));
-
-
-
649 systs_[i]->set_shapes(std::move(copy2_u), std::move(copy2_d), proc_hist);
-
-
-
-
-
-
-
-
-
658 std::vector<double> prev_proc_rates(procs_.size());
-
-
-
-
-
663 std::vector<std::unique_ptr<TH1>> scaled_procs(procs_.size());
-
-
665 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
666 if (procs_[i]->shape()) {
-
-
668 std::unique_ptr<TH1> copy(procs_[i]->ClonedScaledShape());
-
-
670 prev_proc_rates[i] = procs_[i]->no_norm_rate();
-
671 for (
int j = 1; j <= copy->GetNbinsX();j ++){
-
672 if(copy->GetBinLowEdge(j) >= min && copy->GetBinLowEdge(j+1) <= max){
-
673 copy->SetBinContent(j,0.);
-
674 copy->SetBinError(j,0.);
-
-
-
-
678 procs_[i]->set_shape(std::move(copy),
true);
-
679 scaled_procs[i] = procs_[i]->ClonedScaledShape();
-
-
-
682 for (
unsigned i = 0; i < obs_.size(); ++i) {
-
683 if (obs_[i]->shape()) {
-
684 std::unique_ptr<TH1> copy(obs_[i]->ClonedScaledShape());
-
685 for (
int j = 1; j <= copy->GetNbinsX();j ++){
-
686 if(copy->GetBinLowEdge(j) >= min && copy->GetBinLowEdge(j+1) <= max){
-
687 copy->SetBinContent(j,0.);
-
688 copy->SetBinError(j,0.);
-
-
-
691 obs_[i]->set_shape(std::move(copy),
true);
-
-
-
694 for (
unsigned i = 0; i < systs_.size(); ++i) {
-
695 TH1
const* proc_hist =
nullptr;
-
696 double prev_rate = 0.;
-
697 for (
unsigned j = 0; j < procs_.size(); ++j) {
-
-
699 proc_hist = scaled_procs[j].get();
-
700 prev_rate = prev_proc_rates[j];
-
+
+
633 for (
unsigned i = 0; i < systs_.size(); ++i) {
+
634 TH1
const* proc_hist =
nullptr;
+
635 double prev_rate = 0.;
+
636 for (
unsigned j = 0; j < procs_.size(); ++j) {
+
+
638 proc_hist = scaled_procs[j].get();
+
639 prev_rate = prev_proc_rates[j];
+
+
+
642 if (systs_[i]->shape_u() && systs_[i]->shape_d()) {
+
+
644 std::unique_ptr<TH1> copy_u(systs_[i]->ClonedShapeU());
+
645 std::unique_ptr<TH1> copy_d(systs_[i]->ClonedShapeD());
+
+
+
+
+
650 copy_u->Scale(systs_[i]->value_u() * prev_rate);
+
651 copy_d->Scale(systs_[i]->value_d() * prev_rate);
+
+
653 std::unique_ptr<TH1> copy2_u(
+
654 copy_u->Rebin(bins.size() - 1,
"", &(bins[0])));
+
655 std::unique_ptr<TH1> copy2_d(
+
656 copy_d->Rebin(bins.size() - 1,
"", &(bins[0])));
+
+
+
659 systs_[i]->set_shapes(std::move(copy2_u), std::move(copy2_d), proc_hist);
+
+
+
+
+
+
+
+
+
668 std::vector<double> prev_proc_rates(procs_.size());
+
+
+
+
+
673 std::vector<std::unique_ptr<TH1>> scaled_procs(procs_.size());
+
+
675 for (
unsigned i = 0; i < procs_.size(); ++i) {
+
676 if (procs_[i]->shape()) {
+
+
678 std::unique_ptr<TH1> copy(procs_[i]->ClonedScaledShape());
+
+
680 prev_proc_rates[i] = procs_[i]->no_norm_rate();
+
681 for (
int j = 1; j <= copy->GetNbinsX();j ++){
+
682 if(copy->GetBinLowEdge(j) >= min && copy->GetBinLowEdge(j+1) <= max){
+
683 copy->SetBinContent(j,0.);
+
684 copy->SetBinError(j,0.);
+
+
+
+
688 procs_[i]->set_shape(std::move(copy),
true);
+
689 scaled_procs[i] = procs_[i]->ClonedScaledShape();
+
+
+
692 for (
unsigned i = 0; i < obs_.size(); ++i) {
+
693 if (obs_[i]->shape()) {
+
694 std::unique_ptr<TH1> copy(obs_[i]->ClonedScaledShape());
+
695 for (
int j = 1; j <= copy->GetNbinsX();j ++){
+
696 if(copy->GetBinLowEdge(j) >= min && copy->GetBinLowEdge(j+1) <= max){
+
697 copy->SetBinContent(j,0.);
+
698 copy->SetBinError(j,0.);
+
+
+
701 obs_[i]->set_shape(std::move(copy),
true);
-
703 if (systs_[i]->shape_u() && systs_[i]->shape_d()) {
-
-
705 std::unique_ptr<TH1> copy_u(systs_[i]->ClonedShapeU());
-
706 std::unique_ptr<TH1> copy_d(systs_[i]->ClonedShapeD());
-
-
-
-
-
711 copy_u->Scale(systs_[i]->value_u() * prev_rate);
-
712 copy_d->Scale(systs_[i]->value_d() * prev_rate);
-
-
714 for (
int j = 1; j <= copy_u->GetNbinsX();j ++){
-
715 if(copy_u->GetBinLowEdge(j) >= min && copy_u->GetBinLowEdge(j+1) <= max){
-
716 copy_u->SetBinContent(j,0.);
-
717 copy_u->SetBinError(j,0.);
-
-
719 if(copy_d->GetBinLowEdge(j) >= min && copy_d->GetBinLowEdge(j+1) <= max){
-
720 copy_d->SetBinContent(j,0.);
-
721 copy_d->SetBinError(j,0.);
-
+
+
704 for (
unsigned i = 0; i < systs_.size(); ++i) {
+
705 TH1
const* proc_hist =
nullptr;
+
706 double prev_rate = 0.;
+
707 for (
unsigned j = 0; j < procs_.size(); ++j) {
+
+
709 proc_hist = scaled_procs[j].get();
+
710 prev_rate = prev_proc_rates[j];
+
+
+
713 if (systs_[i]->shape_u() && systs_[i]->shape_d()) {
+
+
715 std::unique_ptr<TH1> copy_u(systs_[i]->ClonedShapeU());
+
716 std::unique_ptr<TH1> copy_d(systs_[i]->ClonedShapeD());
+
+
+
+
+
721 copy_u->Scale(systs_[i]->value_u() * prev_rate);
+
722 copy_d->Scale(systs_[i]->value_d() * prev_rate);
-
-
-
726 systs_[i]->set_shapes(std::move(copy_u), std::move(copy_d), proc_hist);
-
-
-
-
-
-
-
733 for (
unsigned i = 0; i < procs_.size(); ++i) {
-
734 std::set<std::string> binning_vars;
-
735 if (procs_[i]->pdf()) {
-
736 RooAbsData
const* data_obj = FindMatchingData(procs_[i].get());
-
737 std::string var_name =
"CMS_th1x";
-
738 if (data_obj) var_name = data_obj->get()->first()->GetName();
-
739 binning_vars.insert(var_name);
-
-
741 for (
auto & it : wspaces_) {
-
742 for (
auto & var : binning_vars) {
-
-
744 dynamic_cast<RooRealVar*>(it.second->var(var.c_str()));
-
745 if (avar) avar->setBins(nbins);
-
-
-
-
-
-
-
-
753 double CombineHarvester::logKappaForX(
double x,
double k_low,
-
754 double k_high)
const {
-
755 if (k_high == 0. || k_low == 0.) {
-
756 if (verbosity_ >= 1) {
-
757 LOGLINE(log(),
"Have kappa=0.0 (scaling ill-defined), returning 1.0");
-
-
-
-
-
762 return (x >= 0 ? std::pow(k_high, x) : std::pow(k_low, -1.0 * x));
-
-
-
-
-
-
-
-
-
771 double logKhi = std::log(k_high);
-
772 double logKlo = -std::log(k_low);
-
773 double avg = 0.5 * (logKhi + logKlo), halfdiff = 0.5 * (logKhi - logKlo);
-
774 double twox = x + x, twox2 = twox * twox;
-
775 double alpha = 0.125 * twox * (twox2 * (3 * twox2 - 10.) + 15.);
-
776 double ret = avg + alpha * halfdiff;
-
777 return std::exp(ret * x);
-
-
-
-
781 std::vector<std::string>
const& patterns) {
-
782 std::vector<boost::regex> rgx;
-
783 for (
auto const& pt : patterns) rgx.emplace_back(pt);
-
784 for (
auto it = params_.begin(); it != params_.end(); ++it) {
-
785 std::string par = it->first;
-
786 auto & groups = it->second->groups();
-
787 if (groups.count(name))
continue;
-
-
-
-
-
-
-
-
795 std::vector<std::string>
const& patterns) {
-
796 std::vector<boost::regex> rgx;
-
797 for (
auto const& pt : patterns) rgx.emplace_back(pt);
-
798 for (
auto it = params_.begin(); it != params_.end(); ++it) {
-
799 std::string par = it->first;
-
800 auto & groups = it->second->groups();
-
801 if (!groups.count(name))
continue;
-
-
-
-
-
-
-
-
809 std::string
const& newname) {
-
810 for (
auto it = params_.begin(); it != params_.end(); ++it) {
-
811 auto & groups = it->second->groups();
-
812 if (groups.count(oldname)) {
-
813 groups.erase(oldname);
-
814 groups.insert(newname);
-
-
-
-
-
-
820 post_lines_.push_back(line);
-
-
-
-
-
-
-
@@ -55,7 +52,7 @@
onblur="searchBox.OnSearchFieldFocus(false)"
onkeyup="searchBox.OnSearchFieldChange(event)"/>
- 
+ 
Inheritance diagram for ch::Object:
-
-
-