Merge pull request #638 from kpedro88/track_errors

add option to track parameter errors

Author: adewit

docs/part3/runningthetool.md

@@ -67,6 +67,11 @@ There are a number of useful command line options which can be used to alter the
     - the name of the branch will be **trackedParam_*name***.
     - the exact behaviour depends on the method. For example, when using `MultiDimFit` with the `--algo scan`, the value of the parameter at each point in the scan will be saved while for `FitDiagnostics`, only the value at the end of the method will be saved.
 
+-   `--trackErrors name1[,name2,...]` will add a branch to the output tree for the error of each of the named parameters. This option supports the use of regexp via by replacing `name` with `rgx{some regular expression}`
+
+    - the name of the branch will be **trackedError_*name***.
+    - the behaviour is the same as `--trackParameters` above.
+
 #### Generic Minimizer Options
 
 Combine uses its own minimizer class which is used to steer Minuit (via RooMinimizer) named the `CascadeMinimizer`. This allows for sequential minimization which can help in case a particular setting/algo fails. Also, the `CascadeMinimizer` knows about extra features of Combine such as *discrete* nuisance parameters.

interface/Combine.h

@@ -6,6 +6,7 @@
 #include <boost/filesystem.hpp>
 #include "RooArgSet.h"
 #include "RooAbsReal.h"
+#include "RooRealVar.h"
 #include <boost/algorithm/string/split.hpp>
 #include <boost/algorithm/string/classification.hpp>
 
@@ -73,6 +74,8 @@ class Combine {
   void addNuisances(const RooArgSet *);
   void addFloatingParameters(const RooArgSet &);
   void addPOI(const RooArgSet *);
+  template <class Var>
+  void addBranches(const std::string&, RooWorkspace*, std::vector<std::pair<Var*,float>>&, const std::string&);
 
   boost::program_options::options_description statOptions_, ioOptions_, miscOptions_;
 
@@ -119,7 +122,9 @@ class Combine {
   static TTree *tree_;
 
   static std::vector<std::pair<RooAbsReal*,float> > trackedParametersMap_;
+  static std::vector<std::pair<RooRealVar*,float> > trackedErrorsMap_;
   static std::string  trackParametersNameString_;
+  static std::string  trackErrorsNameString_;
   static std::string  textToWorkspaceString_;
 
 };

src/Combine.cc

@@ -13,6 +13,7 @@
 #include <algorithm>
 #include <unistd.h>
 #include <errno.h>
+#include <sstream>
 
 #include <TCanvas.h>
 #include <TFile.h>
@@ -92,9 +93,11 @@ std::string setPhysicsModelParameterRangeExpression_ = "";
 std::string defineBackgroundOnlyModelParameterExpression_ = "";
 
 std::string Combine::trackParametersNameString_="";
+std::string Combine::trackErrorsNameString_="";
 std::string Combine::textToWorkspaceString_="";
 
 std::vector<std::pair<RooAbsReal*,float> > Combine::trackedParametersMap_;
+std::vector<std::pair<RooRealVar*,float> > Combine::trackedErrorsMap_;
 
 Combine::Combine() :
     statOptions_("Common statistics options"),
@@ -155,6 +158,7 @@ Combine::Combine() :
       ("genUnbinnedChannels", po::value<std::string>(&genAsUnbinned_)->default_value(genAsUnbinned_), "Flag the given channels to be generated unbinned (irrespectively of how they were flagged at workspace creation)") 
       ("text2workspace",   boost::program_options::value<std::string>(&textToWorkspaceString_)->default_value(""), "Pass along options to text2workspace (default = none)")
       ("trackParameters",   boost::program_options::value<std::string>(&trackParametersNameString_)->default_value(""), "Keep track of parameters in workspace, also accepts regexp with syntax 'rgx{<my regexp>}' (default = none)")
+      ("trackErrors",   boost::program_options::value<std::string>(&trackErrorsNameString_)->default_value(""), "Keep track of errors on parameters in workspace, also accepts regexp with syntax 'rgx{<my regexp>}' (default = none)")
       ; 
 }
 
@@ -812,44 +816,10 @@ void Combine::run(TString hlfFile, const std::string &dataset, double &limit, do
 
   tree_ = tree;
 
-  // Set up additional branches 
-  if(trackParametersNameString_!=""){
-    char tmp[10240] ;
-    strlcpy(tmp,trackParametersNameString_.c_str(),10240) ;
-    char* token = strtok(tmp,",") ;
-    while(token) {
-      if (boost::starts_with(token, "rgx{") && boost::ends_with(token, "}")) {
-          std::string tokenstr(token);
-          std::string reg_esp = tokenstr.substr(4, tokenstr.size()-5);
-          std::cout<<"interpreting "<<reg_esp<<" as regex "<<std::endl;
-          std::regex rgx( reg_esp, std::regex::ECMAScript);
+  // Set up additional branches
+  addBranches(trackParametersNameString_,w,trackedParametersMap_,"Param");
+  addBranches(trackErrorsNameString_,w,trackedErrorsMap_,"Error");
 
-          RooArgSet allParams(w->allVars());
-          std::auto_ptr<TIterator> iter(allParams.createIterator());
-          for (RooAbsArg *a = (RooAbsArg*) iter->Next(); a != 0; a = (RooAbsArg*) iter->Next()) {
-              RooAbsReal *tmp = dynamic_cast<RooAbsReal *>(a);
-              const std::string &target = tmp->GetName();
-              std::smatch match;
-              if (std::regex_match(target, match, rgx)) {
-                  if (tmp->isConstant()) continue;
-                  Combine::trackedParametersMap_.push_back(std::pair<RooAbsReal*,float>(tmp,tmp->getVal()));
-              }
-          }
-          token = strtok(0,",") ;
-      } else {
-          RooAbsReal *a =(RooAbsReal*)w->obj(token); 
-          if (a == 0) throw std::invalid_argument(std::string("Parameter ")+(token)+" not in model.");
-          Combine::trackedParametersMap_.push_back(std::pair<RooAbsReal*,float>(a,a->getVal()));
-          token = strtok(0,",") ;
-      } 
-
-    }
-  }
-  
-  for (std::vector<std::pair<RooAbsReal*,float> >::iterator it = Combine::trackedParametersMap_.begin(); it!=Combine::trackedParametersMap_.end();it++){
-    const char * token = (it->first)->GetName();
-    addBranch((std::string("trackedParam_")+token).c_str(), &(it->second), (std::string("trackedParam_")+token+std::string("/F")).c_str()); 
-  }
   // Should have the PDF at this point, if not something is really odd?
   if (!(mc->GetPdf())){
 	std::cerr << " FATAL ERROR! PDF not found in ModelConfig. \n Try to build the workspace first with text2workspace.py and run with the binary output." << std::endl;
@@ -1190,8 +1160,11 @@ void Combine::commitPoint(bool expected, float quantile) {
     Float_t saveQuantile =  g_quantileExpected_;
     g_quantileExpected_ = quantile;
 
-    for (std::vector<std::pair<RooAbsReal*,float> >::iterator it = Combine::trackedParametersMap_.begin(); it!=Combine::trackedParametersMap_.end();it++){
-	it->second = (it->first)->getVal();
+    for (auto& it : trackedParametersMap_){
+      it.second = (it.first)->getVal();
+    }
+    for (auto& it : trackedErrorsMap_){
+      it.second = (it.first)->getError();
     }
 
     if (g_fillTree_) tree_->Fill();
@@ -1283,3 +1256,41 @@ void Combine::addDiscreteNuisances(RooWorkspace *w){
       }
     }
 }
+
+template <class Var>
+void Combine::addBranches(const std::string& trackString, RooWorkspace* w, std::vector<std::pair<Var*,float>>& trackMap, const std::string& vtype) {
+  if(trackString!=""){
+    std::stringstream ss(trackString);
+    std::string token;
+    while(std::getline(ss,token,',')) {
+      if (boost::starts_with(token, "rgx{") && boost::ends_with(token, "}")) {
+          std::string reg_esp = token.substr(4, token.size()-5);
+          std::cout<<"interpreting "<<reg_esp<<" as regex "<<std::endl;
+          std::regex rgx( reg_esp, std::regex::ECMAScript);
+
+          RooArgSet allParams(w->allVars());
+          std::auto_ptr<TIterator> iter(allParams.createIterator());
+          for (RooAbsArg *a = (RooAbsArg*) iter->Next(); a != 0; a = (RooAbsArg*) iter->Next()) {
+              Var *tmp = dynamic_cast<Var *>(a);
+              if(tmp==nullptr) continue;
+              const std::string &target = tmp->GetName();
+              std::smatch match;
+              if (std::regex_match(target, match, rgx)) {
+                  if (tmp->isConstant()) continue;
+                  trackMap.emplace_back(tmp,0.f);
+              }
+          }
+      } else {
+          Var *a =(Var*)w->obj(token.c_str());
+          if (a == 0) throw std::invalid_argument(vtype+" "+(token)+" not in model.");
+          trackMap.emplace_back(a,0.f);
+      }
+
+    }
+  }
+
+  for (auto& it : trackMap){
+    const char * token = (it.first)->GetName();
+    addBranch((std::string("tracked")+vtype+"_"+token).c_str(), &(it.second), (std::string("tracked")+vtype+"_"+token+std::string("/F")).c_str());
+  }
+}