ConfigurableHisto.h

Go to the documentation of this file.

#ifndef ConfigurableAnalysis_ConfigurableHisto_H
#define ConfigurableAnalysis_ConfigurableHisto_H
 
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "CommonTools/Utils/interface/TFileDirectory.h"
 
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "THStack.h"
 
#include "PhysicsTools/UtilAlgos/interface/VariableHelper.h"
#include "PhysicsTools/UtilAlgos/interface/CachingVariable.h"
 
class ConfigurableAxis {
public:
  ConfigurableAxis() {}
  ConfigurableAxis(const edm::ParameterSet& par) : nBin_(0), Min_(0), Max_(0), Label_("") {
    Label_ = par.getParameter<std::string>("Label");
 
    if (par.exists("nBins")) {
      nBin_ = par.getParameter<unsigned int>("nBins");
      Min_ = par.getParameter<double>("Min");
      Max_ = par.getParameter<double>("Max");
    } else {
      if (par.exists("vBins"))
        vBins_ = par.getParameter<std::vector<double> >("vBins");
      else {
        Min_ = par.getParameter<double>("Min");
        Max_ = par.getParameter<double>("Max");
      }
    }
  }
  bool variableSize() { return !vBins_.empty(); }
  unsigned int nBin() {
    if (variableSize())
      return vBins_.size() - 1;
    else
      return nBin_;
  }
  double Min() {
    if (variableSize())
      return vBins_.front();
    else
      return Min_;
  }
  double Max() {
    if (variableSize())
      return vBins_.back();
    else
      return Max_;
  }
  const std::string& Label() { return Label_; }
  const double* xBins() {
    if (!vBins_.empty())
      return &(vBins_.front());
    else
      return nullptr;
  }
 
private:
  std::vector<double> vBins_;
  unsigned int nBin_;
  double Min_;
  double Max_;
  std::string Label_;
};
 
class ConfigurableHisto {
public:
  enum HType { h1, h2, prof };
  ConfigurableHisto(HType t, std::string name, edm::ParameterSet& iConfig)
      : type_(t), h_(nullptr), name_(name), conf_(iConfig), x_(nullptr), y_(nullptr), z_(nullptr), w_(nullptr) {}
 
  virtual ~ConfigurableHisto() {}
 
  virtual ConfigurableHisto* clone() const { return new ConfigurableHisto(*this); }
 
  virtual void book(TFileDirectory* dir) {
    std::string title = conf_.getParameter<std::string>("title");
    edm::ParameterSet xAxisPSet = conf_.getParameter<edm::ParameterSet>("xAxis");
    ConfigurableAxis xAxis(xAxisPSet);
    x_ = edm::Service<VariableHelperService>()->get().variable(xAxisPSet.getParameter<std::string>("var"));
 
    std::string yLabel = "";
    bool yVBin = false;
    ConfigurableAxis yAxis;
    if (conf_.exists("yAxis")) {
      edm::ParameterSet yAxisPSet = conf_.getParameter<edm::ParameterSet>("yAxis");
      yAxis = ConfigurableAxis(yAxisPSet);
      yLabel = yAxis.Label();
      //at least TH2 or TProfile
      if (yAxisPSet.exists("var"))
        y_ = edm::Service<VariableHelperService>()->get().variable(yAxisPSet.getParameter<std::string>("var"));
      yVBin = yAxis.variableSize();
    }
 
    if (conf_.exists("zAxis")) {
      throw;
    }
 
    bool xVBin = xAxis.variableSize();
 
    if (type() == h1) {
      //make TH1F
      if (xVBin)
        h_ = dir->make<TH1F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins());
      else
        h_ = dir->make<TH1F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.Min(), xAxis.Max());
    } else if (type() == h2) {
      //make TH2F
      if (xVBin) {
        if (yVBin)
          h_ = dir->make<TH2F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.nBin(), yAxis.xBins());
        else
          h_ = dir->make<TH2F>(
              name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.nBin(), yAxis.Min(), yAxis.Max());
      } else {
        if (yVBin)
          h_ = dir->make<TH2F>(
              name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.Min(), xAxis.Max(), yAxis.nBin(), yAxis.xBins());
        else
          h_ = dir->make<TH2F>(name_.c_str(),
                               title.c_str(),
                               xAxis.nBin(),
                               xAxis.Min(),
                               xAxis.Max(),
                               yAxis.nBin(),
                               yAxis.Min(),
                               yAxis.Max());
      }
    } else if (type() == prof) {
      //make TProfile
      std::string pFopt = "";
      if (conf_.exists("Option"))
        pFopt = conf_.getParameter<std::string>("Option");
 
      if (xVBin)
        h_ = dir->make<TProfile>(
            name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.Min(), yAxis.Max(), pFopt.c_str());
      else
        h_ = dir->make<TProfile>(name_.c_str(),
                                 title.c_str(),
                                 xAxis.nBin(),
                                 xAxis.Min(),
                                 xAxis.Max(),
                                 yAxis.Min(),
                                 yAxis.Max(),
                                 pFopt.c_str());
 
    } else {
      edm::LogError("ConfigurableHisto") << "cannot book: " << name_ << "\n" << conf_.dump();
      throw;
    }
 
    //cosmetics
    h_->GetXaxis()->SetTitle(xAxis.Label().c_str());
    h_->SetYTitle(yLabel.c_str());
 
    if (conf_.exists("weight")) {
      w_ = edm::Service<VariableHelperService>()->get().variable(conf_.getParameter<std::string>("weight"));
    }
  }
 
  virtual void fill(const edm::Event& iEvent) {
    if (!h_)
      return;
 
    double weight = 1.0;
    if (w_) {
      if (w_->compute(iEvent))
        weight = (*w_)(iEvent);
      else {
        edm::LogInfo("ConfigurableHisto") << "could not compute the weight for: " << name_ << " with config:\n"
                                          << conf_.dump() << " default to 1.0";
      }
    }
 
    TProfile* pcast(nullptr);
    TH2* h2cast(nullptr);
    switch (type_) {
      case h1:
        if (!h_)
          throw;
        if (x_->compute(iEvent))
          h_->Fill((*x_)(iEvent), weight);
        else {
          edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
        }
        break;
      case prof:
        pcast = dynamic_cast<TProfile*>(h_);
        if (!pcast)
          throw;
        if (x_->compute(iEvent) && y_->compute(iEvent))
          pcast->Fill((*x_)(iEvent), (*y_)(iEvent), weight);
        else {
          edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
        }
        break;
      case h2:
        h2cast = dynamic_cast<TH2*>(h_);
        if (!h2cast)
          throw;
        if (x_->compute(iEvent) && y_->compute(iEvent))
          h2cast->Fill((*x_)(iEvent), (*y_)(iEvent), weight);
        else {
          edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
        }
        break;
    }
  }
  const HType& type() { return type_; }
 
  void complete() {}
  TH1* h() { return h_; }
 
protected:
  ConfigurableHisto(const ConfigurableHisto& master) {
    type_ = master.type_;
    h_ = nullptr;  //no histogram attached in copy constructor
    name_ = master.name_;
    conf_ = master.conf_;
    x_ = master.x_;
    y_ = master.y_;
    z_ = master.z_;
    w_ = master.w_;
  }
  HType type_;
  TH1* h_;
  std::string name_;
  edm::ParameterSet conf_;
 
  const CachingVariable* x_;
  const CachingVariable* y_;
  const CachingVariable* z_;
  const CachingVariable* w_;
};
 
class SplittingConfigurableHisto : public ConfigurableHisto {
public:
  SplittingConfigurableHisto(HType t, std::string name, edm::ParameterSet& pset)
      : ConfigurableHisto(t, name, pset), splitter_(nullptr) {
    std::string title = pset.getParameter<std::string>("title");
 
    //allow for many splitters ...
    if (pset.exists("splitters")) {
      //you want more than one splitter
      std::vector<std::string> splitters = pset.getParameter<std::vector<std::string> >("splitters");
      for (unsigned int s = 0; s != splitters.size(); ++s) {
        const CachingVariable* v = edm::Service<VariableHelperService>()->get().variable(splitters[s]);
        const Splitter* splitter = dynamic_cast<const Splitter*>(v);
        if (!splitter) {
          edm::LogError("SplittingConfigurableHisto")
              << "for: " << name_ << " the splitting variable: " << splitters[s] << " is not a Splitter";
          continue;
        }
 
        //insert in the map
        std::vector<ConfigurableHisto*>& insertedHisto = subHistoMap_[splitter];
        //now configure the histograms
        unsigned int mSlots = splitter->maxSlots();
        for (unsigned int i = 0; i != mSlots; ++i) {
          //---   std::cout<<" slot: "<<i<<std::endl;
          const std::string& slabel = splitter->shortLabel(i);
          const std::string& label = splitter->label(i);
          edm::ParameterSet mPset = pset;
          edm::Entry e("string", title + " for " + label, true);
          mPset.insert(true, "title", e);
          insertedHisto.push_back(new ConfigurableHisto(t, name + slabel, mPset));
        }  //loop on slots
      }    //loop on splitters
 
    }  //if splitters exists
    else {
      //single splitter
      //get the splitting variable
      const CachingVariable* v =
          edm::Service<VariableHelperService>()->get().variable(pset.getParameter<std::string>("splitter"));
      splitter_ = dynamic_cast<const Splitter*>(v);
      if (!splitter_) {
        edm::LogError("SplittingConfigurableHisto")
            << "for: " << name_ << " the splitting variable: " << v->name() << " is not a Splitter";
      } else {
        //configure the splitted plots
        unsigned int mSlots = splitter_->maxSlots();
        for (unsigned int i = 0; i != mSlots; i++) {
          const std::string& slabel = splitter_->shortLabel(i);
          const std::string& label = splitter_->label(i);
          edm::ParameterSet mPset = pset;
          edm::Entry e("string", title + " for " + label, true);
          mPset.insert(true, "title", e);
          subHistos_.push_back(new ConfigurableHisto(t, name + slabel, mPset));
        }
      }
    }
  }  //end of ctor
 
  void book(TFileDirectory* dir) override {
    //book the base histogram
    ConfigurableHisto::book(dir);
 
    if (!subHistoMap_.empty()) {
      SubHistoMap::iterator i = subHistoMap_.begin();
      SubHistoMap::iterator i_end = subHistoMap_.end();
      for (; i != i_end; ++i) {
        for (unsigned int h = 0; h != i->second.size(); ++h) {
          i->second[h]->book(dir);
        }
        //book the THStack
        std::string sName = name_ + "_" + i->first->name();
        std::string sTitle = "Stack histogram of " + name_ + " for splitter " + i->first->name();
        subHistoStacks_[i->first] = dir->make<THStack>(sName.c_str(), sTitle.c_str());
      }
    } else {
      for (unsigned int h = 0; h != subHistos_.size(); h++) {
        subHistos_[h]->book(dir);
      }
      //book a THStack
      std::string sName = name_ + "_" + splitter_->name();
      std::string sTitle = "Stack histogram of " + name_ + " for splitter " + splitter_->name();
      stack_ = dir->make<THStack>(sName.c_str(), sTitle.c_str());
    }
  }
 
  ConfigurableHisto* clone() const override { return new SplittingConfigurableHisto(*this); }
 
  void fill(const edm::Event& e) override {
    //fill the base histogram
    ConfigurableHisto::fill(e);
 
    if (!subHistoMap_.empty()) {
      SubHistoMap::iterator i = subHistoMap_.begin();
      SubHistoMap::iterator i_end = subHistoMap_.end();
      for (; i != i_end; ++i) {
        const Splitter* splitter = i->first;
        if (!splitter)
          continue;
        if (!splitter->compute(e))
          continue;
        unsigned int slot = (unsigned int)(*splitter)(e);
        if (slot >= i->second.size()) {
          edm::LogError("SplittingConfigurableHisto")
              << "slot index: " << slot << " is bigger than slots size: " << i->second.size()
              << " from variable value: " << (*splitter)(e);
          continue;
        }
        //fill in the proper slot
        i->second[slot]->fill(e);
      }
    } else {
      //fill the component histograms
      if (!splitter_)
        return;
      if (!splitter_->compute(e)) {
        return;
      }
      unsigned int slot = (unsigned int)(*splitter_)(e);
      if (slot >= subHistos_.size()) {
        edm::LogError("SplittingConfigurableHisto")
            << "slot index: " << slot << " is bigger than slots size: " << subHistos_.size()
            << " from variable value: " << (*splitter_)(e);
        return;
      }
      subHistos_[slot]->fill(e);
    }
  }
 
  void complete() {
    if (!subHistoMap_.empty()) {
      //fill up the stacks
      SubHistoMap::iterator i = subHistoMap_.begin();
      SubHistoMap::iterator i_end = subHistoMap_.end();
      for (; i != i_end; ++i) {
        for (unsigned int h = 0; h != i->second.size(); h++) {
          //      if (i->second[h]->h()->Integral==0) continue;// do not add empty histograms. NO, because it will be tough to merge two THStack
          subHistoStacks_[i->first]->Add(i->second[h]->h(), i->first->label(h).c_str());
        }
      }
 
    } else {
      //fill up the only stack
      for (unsigned int i = 0; i != subHistos_.size(); i++) {
        stack_->Add(subHistos_[i]->h(), splitter_->label(i).c_str());
      }
    }
  }
 
private:
  SplittingConfigurableHisto(const SplittingConfigurableHisto& master) : ConfigurableHisto(master) {
    splitter_ = master.splitter_;
    if (!master.subHistoMap_.empty()) {
      SubHistoMap::const_iterator i = master.subHistoMap_.begin();
      SubHistoMap::const_iterator i_end = master.subHistoMap_.end();
      for (; i != i_end; ++i) {
        const std::vector<ConfigurableHisto*>& masterHistos = i->second;
        std::vector<ConfigurableHisto*>& clonedHistos = subHistoMap_[i->first];
        for (unsigned int i = 0; i != masterHistos.size(); i++) {
          clonedHistos.push_back(masterHistos[i]->clone());
        }
      }
    } else {
      for (unsigned int i = 0; i != master.subHistos_.size(); i++) {
        subHistos_.push_back(master.subHistos_[i]->clone());
      }
    }
  }
 
  typedef std::map<const Splitter*, std::vector<ConfigurableHisto*> > SubHistoMap;
  typedef std::map<const Splitter*, THStack*> SubHistoStacks;
  SubHistoStacks subHistoStacks_;
  SubHistoMap subHistoMap_;
 
  const Splitter* splitter_;
  std::vector<ConfigurableHisto*> subHistos_;
  // do you want to have a stack already made from the splitted histos?
  THStack* stack_;
};
 
#endif