CombineHarvester.h

Go to the documentation of this file.

#ifndef CombineTools_CombineHarvester_h
#define CombineTools_CombineHarvester_h
#include <string>
#include <utility>
#include <vector>
#include <memory>
#include <map>
#include <unordered_map>
#include <list>
#include <cmath>
#include <set>
#include <functional>
#include "boost/range/algorithm_ext/erase.hpp"
#include "TFile.h"
#include "TH1.h"
#include "RooWorkspace.h"
#include "CombineHarvester/CombineTools/interface/Process.h"
#include "CombineHarvester/CombineTools/interface/Systematic.h"
#include "CombineHarvester/CombineTools/interface/Parameter.h"
#include "CombineHarvester/CombineTools/interface/Observation.h"
#include "CombineHarvester/CombineTools/interface/Utilities.h"
#include "CombineHarvester/CombineTools/interface/HistMapping.h"
 
 
namespace ch {
 
// Define some useful CombineHarvester-specific typedefs
typedef std::vector<std::pair<int, std::string>> Categories;
 
class CombineHarvester {
 public:
  CombineHarvester();
  ~CombineHarvester();
 
  CombineHarvester(CombineHarvester const& other);
  CombineHarvester(CombineHarvester&& other);
  CombineHarvester& operator=(CombineHarvester other);
 
  void SetFlag(std::string const& flag, bool const& value);
  bool GetFlag(std::string const& flag) const;
 
  CombineHarvester cp();
 
  CombineHarvester deep();
  CombineHarvester& PrintAll();
  CombineHarvester& PrintObs();
  CombineHarvester& PrintProcs();
  CombineHarvester& PrintSysts();
  CombineHarvester& PrintParams();
  inline void SetVerbosity(unsigned verbosity) { verbosity_ = verbosity; }
  inline unsigned Verbosity() { return verbosity_; }
  int ParseDatacard(std::string const& filename,
      std::string const& analysis,
      std::string const& era,
      std::string const& channel,
      int bin_id,
      std::string const& mass);
  int ParseDatacard(std::string const& filename,
      std::string parse_rule = "");
 
  void WriteDatacard(std::string const& name, std::string const& root_file);
  void WriteDatacard(std::string const& name, TFile & root_file);
  void WriteDatacard(std::string const& name);
  CombineHarvester& bin(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& bin_id(std::vector<int> const& vec, bool cond = true);
  CombineHarvester& process(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& analysis(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& era(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& channel(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& mass(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& attr(std::vector<std::string> const& vec,std::string attr_label, bool cond = true);
  CombineHarvester& syst_name(std::vector<std::string> const& vec, bool cond = true);
  CombineHarvester& syst_type(std::vector<std::string> const& vec, bool cond = true);
 
  CombineHarvester& process_rgx(std::vector<std::string> const& vec, bool cond = true);
 
  CombineHarvester& signals();
  CombineHarvester& backgrounds();
  CombineHarvester& histograms();
  CombineHarvester& pdfs();
  CombineHarvester& data();
 
  template<typename Function>
  CombineHarvester& FilterAll(Function func);
  template<typename Function>
  CombineHarvester& FilterObs(Function func);
  template<typename Function>
  CombineHarvester& FilterProcs(Function func);
  template<typename Function>
  CombineHarvester& FilterSysts(Function func);
  // Set generation
  std::set<std::string> bin_set();
  std::set<int> bin_id_set();
  std::set<std::string> process_set();
  std::set<std::string> analysis_set();
  std::set<std::string> era_set();
  std::set<std::string> channel_set();
  std::set<std::string> mass_set();
  std::set<std::string> syst_name_set();
  std::set<std::string> syst_type_set();
 
  template <typename T,
            typename R = typename std::decay<
                typename std::result_of<T(Object const*)>::type>::type>
  std::set<R> SetFromAll(T func);
 
  template <typename T,
            typename R = typename std::decay<
                typename std::result_of<T(Observation const*)>::type>::type>
  std::set<R> SetFromObs(T func);
 
  template <typename T,
            typename R = typename std::decay<
                typename std::result_of<T(Process const*)>::type>::type>
  std::set<R> SetFromProcs(T func);
 
  template <typename T,
            typename R = typename std::decay<
                typename std::result_of<T(Systematic const*)>::type>::type>
  std::set<R> SetFromSysts(T func);
  ch::Parameter const* GetParameter(std::string const& name) const;
  ch::Parameter* GetParameter(std::string const& name);
 
  void UpdateParameters(std::vector<ch::Parameter> const& params);
 
  void UpdateParameters(RooFitResult const* fit);
  void UpdateParameters(RooFitResult const& fit);
 
  std::vector<ch::Parameter> GetParameters() const;
  void RenameParameter(std::string const& oldname, std::string const& newname);
 
  template<typename Function>
  void ForEachObj(Function func);
 
  template<typename Function>
  void ForEachProc(Function func);
 
  template<typename Function>
  void ForEachObs(Function func);
 
  template<typename Function>
  void ForEachSyst(Function func);
 
  void VariableRebin(std::vector<double> bins);
  void ZeroBins(double min, double max);
  void SetPdfBins(unsigned nbins);
 
  // 
  double getParFromWs(const std::string name);
  void setParInWs(const std::string name,double value) ;
  void renameParInWs(const std::string& name, const std::string& newName,const std::string& wsName="");
 
  void SetGroup(std::string const& name, std::vector<std::string> const& patterns);
 
  void RemoveGroup(std::string const& name, std::vector<std::string> const& patterns);
 
  void RenameGroup(std::string const& oldname, std::string const& newname);
 
  void AddDatacardLineAtEnd(std::string const& line);
 
  void ClearDatacardLinesAtEnd();
  double GetRate();
  double GetObservedRate();
  double GetUncertainty();
 
  double GetUncertainty(RooFitResult const* fit, unsigned n_samples);
  double GetUncertainty(RooFitResult const& fit, unsigned n_samples);
  TH1F GetShape();
  TH1F GetShapeWithUncertainty();
 
  TH1F GetShapeWithUncertainty(RooFitResult const* fit, unsigned n_samples);
  TH1F GetShapeWithUncertainty(RooFitResult const& fit, unsigned n_samples);
  TH1F GetObservedShape();
 
  TH2F GetRateCovariance(RooFitResult const& fit, unsigned n_samples);
  TH2F GetRateCorrelation(RooFitResult const& fit, unsigned n_samples);
  void AddObservations(std::vector<std::string> mass,
                       std::vector<std::string> analysis,
                       std::vector<std::string> era,
                       std::vector<std::string> channel,
                       ch::Categories bin);
 
  void AddProcesses(std::vector<std::string> mass,
                    std::vector<std::string> analysis,
                    std::vector<std::string> era,
                    std::vector<std::string> channel,
                    std::vector<std::string> procs,
                    ch::Categories bin, bool signal);
  void AddSystVar(std::string const& name, double val, double err);
  void AddSystFromProc(Process const& proc, std::string const& name,
                       std::string const& type, bool asymm, double val_u,
                       double val_d, std::string const& formula,
                       std::string const& args);
 
  template <class Map>
  void AddSyst(CombineHarvester & target, std::string const& name,
               std::string const& type, Map const& valmap);
 
  void ExtractShapes(std::string const& file, std::string const& rule,
                     std::string const& syst_rule);
  void ExtractPdfs(CombineHarvester& target, std::string const& ws_name,
                   std::string const& rule, std::string norm_rule = "");
  void ExtractData(std::string const& ws_name, std::string const& rule);
 
  void AddWorkspace(RooWorkspace const& ws, bool can_rename = false);
 
  void InsertObservation(ch::Observation const& obs);
  void InsertProcess(ch::Process const& proc);
  void InsertSystematic(ch::Systematic const& sys);
 
  void RenameSystematic(CombineHarvester& target, std::string const& old_name, std::string const& new_name);
 
  void CreateParameterIfEmpty(std::string const& name);
 
  void AddBinByBin(double threshold, bool fixed_norm, CombineHarvester* other);
 
 
  void AddBinByBin(double threshold, bool fixed_norm, CombineHarvester & other);
 
 
  void MergeBinErrors(double bbb_threshold, double merge_threshold);
  void SetAutoMCStats(CombineHarvester &target, double thresh, bool sig=false, int mode=1);
  void RenameAutoMCStatsBin(std::string const& oldname, std::string const& newname);
  std::set<std::string> GetAutoMCStatsBins() const;
 
  void AddExtArgValue(std::string const& name, double const& value);
 private:
  friend void swap(CombineHarvester& first, CombineHarvester& second);
 
  // ---------------------------------------------------------------
  // Main data members
  // ---------------------------------------------------------------
  std::vector<std::shared_ptr<Observation>> obs_;
  std::vector<std::shared_ptr<Process>> procs_;
  std::vector<std::shared_ptr<Systematic>> systs_;
  std::map<std::string, std::shared_ptr<Parameter>> params_;
  std::map<std::string, std::shared_ptr<RooWorkspace>> wspaces_;
 
  std::unordered_map<std::string, bool> flags_;
 
  struct AutoMCStatsSettings {
    double event_threshold;
    bool include_signal;
    int hist_mode;
 
    AutoMCStatsSettings(double thresh, bool sig=false, int mode=1) {
      event_threshold = thresh;
      include_signal = sig;
      hist_mode = mode;
    }
 
    AutoMCStatsSettings() : AutoMCStatsSettings(0.) {}
  };
 
  std::map<std::string, AutoMCStatsSettings> auto_stats_settings_;
  std::vector<std::string> post_lines_;
 
  // ---------------------------------------------------------------
  // typedefs
  // ---------------------------------------------------------------
  typedef std::pair<std::string, std::string> StrPair;
  typedef std::vector<StrPair> StrPairVec;
  typedef std::vector<std::string> StrVec;
 
 
  // ---------------------------------------------------------------
  // Logging
  // ---------------------------------------------------------------
  unsigned verbosity_;
  std::ostream * log_;
  std::ostream& log() const { return *log_; }
 
  // ---------------------------------------------------------------
  // Private methods for the shape extraction routines
  // --> implementation in src/CombineHarvester.cc
  // ---------------------------------------------------------------
  void LoadShapes(Observation* entry,
                     std::vector<HistMapping> const& mappings);
  void LoadShapes(Process* entry,
                     std::vector<HistMapping> const& mappings);
  void LoadShapes(Systematic* entry,
                     std::vector<HistMapping> const& mappings);
 
  HistMapping const& ResolveMapping(std::string const& process,
                                    std::string const& bin,
                                    std::vector<HistMapping> const& mappings);
 
  StrPairVec GenerateShapeMapAttempts(std::string process,
      std::string category);
 
  std::shared_ptr<RooWorkspace> SetupWorkspace(RooWorkspace const& ws,
                                    bool can_rename = false);
 
  void ImportParameters(RooArgSet *vars);
 
  RooAbsData const* FindMatchingData(Process const* proc);
 
  ch::Parameter * SetupRateParamVar(std::string const& name, double val, bool is_ext_arg = false);
  void SetupRateParamFunc(std::string const& name, std::string const& formula,
                          std::string const& pars);
  void SetupRateParamWspObj(std::string const& name, std::string const& obj, bool is_ext_arg = false);
  bool SetupRateParamWspObjFromWsStore(std::string const& name, std::string const& obj,                                                                                               std::map<std::string, std::shared_ptr<RooWorkspace>> const& ws_store);
  // ---------------------------------------------------------------
  // Private methods for the shape writing routines
  // ---------------------------------------------------------------
  void WriteHistToFile(
      TH1 * hist,
      TFile * file,
      std::vector<HistMapping> const& mappings,
      std::string const& bin,
      std::string const& process,
      std::string const& mass,
      std::string const& nuisance,
      unsigned type);
 
void FillHistMappings(std::vector<HistMapping> & mappings);
 
  // ---------------------------------------------------------------
  // Private methods for shape/yield evaluation
  // --> implementation in src/CombineHarvester_Evaluate.cc
  // ---------------------------------------------------------------
  typedef std::vector<std::vector<Systematic const*>> ProcSystMap;
  ProcSystMap GenerateProcSystMap();
 
  double GetRateInternal(ProcSystMap const& lookup,
    std::string const& single_sys = "");
 
  TH1F GetShapeInternal(ProcSystMap const& lookup,
    std::string const& single_sys = "");
 
  inline double smoothStepFunc(double x) const {
    if (std::fabs(x) >= 1.0/*_smoothRegion*/) return x > 0 ? +1 : -1;
    double xnorm = x/1.0;/*_smoothRegion*/
    double xnorm2 = xnorm*xnorm;
    return 0.125 * xnorm * (xnorm2 * (3.*xnorm2 - 10.) + 15);
  }
 
  double logKappaForX(double x, double k_low, double k_high) const;
 
  void ShapeDiff(double x, TH1F* target, TH1 const* nom, TH1 const* low,
                 TH1 const* high, bool linear);
  void ShapeDiff(double x, TH1F* target, RooDataHist const* nom,
                 RooDataHist const* low, RooDataHist const* high);
 
  // bug fix for RooConstVar compatibility between ROOT626 and workspace created with earlier versions
  RooWorkspace* fixRooConstVar(RooWorkspace *win, bool useRooRealVar=true, bool clean=true);
};
 
 
// ---------------------------------------------------------------
// Template method implementation
// ---------------------------------------------------------------
template <typename T, typename R>
std::set<R> CombineHarvester::SetFromAll(T func) {
  std::set<R> ret;
  for (auto const& item : obs_) ret.insert(func(item.get()));
  for (auto const& item : procs_) ret.insert(func(item.get()));
  for (auto const& item : systs_) ret.insert(func(item.get()));
  return ret;
};
 
template <typename T, typename R>
std::set<R> CombineHarvester::SetFromObs(T func) {
  std::set<R> ret;
  for (auto const& item : obs_) ret.insert(func(item.get()));
  return ret;
};
 
template <typename T, typename R>
std::set<R> CombineHarvester::SetFromProcs(T func) {
  std::set<R> ret;
  for (auto const& item : procs_) ret.insert(func(item.get()));
  return ret;
};
 
template <typename T, typename R>
std::set<R> CombineHarvester::SetFromSysts(T func) {
  std::set<R> ret;
  for (auto const& item : systs_) ret.insert(func(item.get()));
  return ret;
};
 
template<typename Function>
void CombineHarvester::ForEachObj(Function func) {
  ForEachObs(func);
  ForEachProc(func);
  ForEachSyst(func);
}
 
template<typename Function>
void CombineHarvester::ForEachProc(Function func) {
  for (auto & item: procs_) func(item.get());
}
 
template<typename Function>
void CombineHarvester::ForEachObs(Function func) {
  for (auto & item: obs_) func(item.get());
}
 
template<typename Function>
void CombineHarvester::ForEachSyst(Function func) {
  for (auto & item: systs_) func(item.get());
}
 
template<typename Function>
CombineHarvester& CombineHarvester::FilterAll(Function func) {
  FilterObs(func);
  FilterProcs(func);
  FilterSysts(func);
  return *this;
}
 
template<typename Function>
CombineHarvester& CombineHarvester::FilterObs(Function func) {
  boost::remove_erase_if(
      obs_, [&](std::shared_ptr<Observation> ptr) { return func(ptr.get());
  });
  return *this;
}
 
template<typename Function>
CombineHarvester& CombineHarvester::FilterProcs(Function func) {
  boost::remove_erase_if(
      procs_, [&](std::shared_ptr<Process> ptr) { return func(ptr.get());
  });
  return *this;
}
template<typename Function>
CombineHarvester& CombineHarvester::FilterSysts(Function func) {
  boost::remove_erase_if(
      systs_, [&](std::shared_ptr<Systematic> ptr) { return func(ptr.get());
  });
  return *this;
}
 
template <class Map>
void CombineHarvester::AddSyst(CombineHarvester& target,
                               std::string const& name, std::string const& type,
                               Map const& valmap) {
  // Keep track of which Process entries get a Systematic assigned and which
  // don't. If verbosity is on we'll print lists of these processes at the end.
  std::vector<ch::Process *> not_added_procs;
  std::vector<ch::Process *> added_procs;
  // Also track which tuples in the map did not get used. Do this by getting the
  // full map here and then removing elements as they are used to create a
  // Systematic.
  auto tuples = valmap.GetTupleSet();
  if (verbosity_ >= 1) {
    log() << (name + ":" + type) << "\n";
  }
  for (unsigned i = 0; i < procs_.size(); ++i) {
    if (!valmap.Contains(procs_[i].get())) {
      not_added_procs.push_back(procs_[i].get());
      continue;
    }
    tuples.erase(valmap.GetTuple(procs_[i].get()));
    added_procs.push_back(procs_[i].get());
    double val_u = valmap.ValU(procs_[i].get());
    double val_d = valmap.ValD(procs_[i].get());
    std::string formula = valmap.Formula(procs_[i].get());
    std::string args = valmap.Args(procs_[i].get());
    target.AddSystFromProc(*(procs_[i]), name, type, valmap.IsAsymm(),
                           val_u, val_d, formula, args);
  }
  if (tuples.size() && verbosity_ >= 1) {
    log() << ">> Map keys that were not used to create a Systematic:\n";
    for (auto s : tuples) {
      log() << ch::Tuple2String(s) << "\n";
    }
  }
  if (verbosity_ >= 2) {
    Process::PrintHeader(log());
    log() << ">> Process entries that did not get a Systematic:\n";
    for (auto p : not_added_procs) {
      log() << *p << "\n";
    }
    log() << ">> Process entries that did get a Systematic:\n";
    for (auto p : added_procs) {
      log() << *p << "\n";
    }
  }
}
}
 
#endif