Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:24:17

0001 #ifndef ConfigurableAnalysis_ConfigurableHisto_H
0002 #define ConfigurableAnalysis_ConfigurableHisto_H
0003 
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "CommonTools/Utils/interface/TFileDirectory.h"
0006 
0007 #include "TH1.h"
0008 #include "TH1F.h"
0009 #include "TH2F.h"
0010 #include "TProfile.h"
0011 #include "THStack.h"
0012 
0013 #include "PhysicsTools/UtilAlgos/interface/VariableHelper.h"
0014 #include "PhysicsTools/UtilAlgos/interface/CachingVariable.h"
0015 
0016 class ConfigurableAxis {
0017 public:
0018   ConfigurableAxis() {}
0019   ConfigurableAxis(const edm::ParameterSet& par) : nBin_(0), Min_(0), Max_(0), Label_("") {
0020     Label_ = par.getParameter<std::string>("Label");
0021 
0022     if (par.exists("nBins")) {
0023       nBin_ = par.getParameter<unsigned int>("nBins");
0024       Min_ = par.getParameter<double>("Min");
0025       Max_ = par.getParameter<double>("Max");
0026     } else {
0027       if (par.exists("vBins"))
0028         vBins_ = par.getParameter<std::vector<double> >("vBins");
0029       else {
0030         Min_ = par.getParameter<double>("Min");
0031         Max_ = par.getParameter<double>("Max");
0032       }
0033     }
0034   }
0035   bool variableSize() { return !vBins_.empty(); }
0036   unsigned int nBin() {
0037     if (variableSize())
0038       return vBins_.size() - 1;
0039     else
0040       return nBin_;
0041   }
0042   double Min() {
0043     if (variableSize())
0044       return vBins_.front();
0045     else
0046       return Min_;
0047   }
0048   double Max() {
0049     if (variableSize())
0050       return vBins_.back();
0051     else
0052       return Max_;
0053   }
0054   const std::string& Label() { return Label_; }
0055   const double* xBins() {
0056     if (!vBins_.empty())
0057       return &(vBins_.front());
0058     else
0059       return nullptr;
0060   }
0061 
0062 private:
0063   std::vector<double> vBins_;
0064   unsigned int nBin_;
0065   double Min_;
0066   double Max_;
0067   std::string Label_;
0068 };
0069 
0070 class ConfigurableHisto {
0071 public:
0072   enum HType { h1, h2, prof };
0073   ConfigurableHisto(HType t, std::string name, edm::ParameterSet& iConfig)
0074       : type_(t), h_(nullptr), name_(name), conf_(iConfig), x_(nullptr), y_(nullptr), z_(nullptr), w_(nullptr) {}
0075 
0076   virtual ~ConfigurableHisto() {}
0077 
0078   virtual ConfigurableHisto* clone() const { return new ConfigurableHisto(*this); }
0079 
0080   virtual void book(TFileDirectory* dir) {
0081     std::string title = conf_.getParameter<std::string>("title");
0082     edm::ParameterSet xAxisPSet = conf_.getParameter<edm::ParameterSet>("xAxis");
0083     ConfigurableAxis xAxis(xAxisPSet);
0084     x_ = edm::Service<VariableHelperService>()->get().variable(xAxisPSet.getParameter<std::string>("var"));
0085 
0086     std::string yLabel = "";
0087     bool yVBin = false;
0088     ConfigurableAxis yAxis;
0089     if (conf_.exists("yAxis")) {
0090       edm::ParameterSet yAxisPSet = conf_.getParameter<edm::ParameterSet>("yAxis");
0091       yAxis = ConfigurableAxis(yAxisPSet);
0092       yLabel = yAxis.Label();
0093       //at least TH2 or TProfile
0094       if (yAxisPSet.exists("var"))
0095         y_ = edm::Service<VariableHelperService>()->get().variable(yAxisPSet.getParameter<std::string>("var"));
0096       yVBin = yAxis.variableSize();
0097     }
0098 
0099     if (conf_.exists("zAxis")) {
0100       throw;
0101     }
0102 
0103     bool xVBin = xAxis.variableSize();
0104 
0105     if (type() == h1) {
0106       //make TH1F
0107       if (xVBin)
0108         h_ = dir->make<TH1F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins());
0109       else
0110         h_ = dir->make<TH1F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.Min(), xAxis.Max());
0111     } else if (type() == h2) {
0112       //make TH2F
0113       if (xVBin) {
0114         if (yVBin)
0115           h_ = dir->make<TH2F>(name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.nBin(), yAxis.xBins());
0116         else
0117           h_ = dir->make<TH2F>(
0118               name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.nBin(), yAxis.Min(), yAxis.Max());
0119       } else {
0120         if (yVBin)
0121           h_ = dir->make<TH2F>(
0122               name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.Min(), xAxis.Max(), yAxis.nBin(), yAxis.xBins());
0123         else
0124           h_ = dir->make<TH2F>(name_.c_str(),
0125                                title.c_str(),
0126                                xAxis.nBin(),
0127                                xAxis.Min(),
0128                                xAxis.Max(),
0129                                yAxis.nBin(),
0130                                yAxis.Min(),
0131                                yAxis.Max());
0132       }
0133     } else if (type() == prof) {
0134       //make TProfile
0135       std::string pFopt = "";
0136       if (conf_.exists("Option"))
0137         pFopt = conf_.getParameter<std::string>("Option");
0138 
0139       if (xVBin)
0140         h_ = dir->make<TProfile>(
0141             name_.c_str(), title.c_str(), xAxis.nBin(), xAxis.xBins(), yAxis.Min(), yAxis.Max(), pFopt.c_str());
0142       else
0143         h_ = dir->make<TProfile>(name_.c_str(),
0144                                  title.c_str(),
0145                                  xAxis.nBin(),
0146                                  xAxis.Min(),
0147                                  xAxis.Max(),
0148                                  yAxis.Min(),
0149                                  yAxis.Max(),
0150                                  pFopt.c_str());
0151 
0152     } else {
0153       edm::LogError("ConfigurableHisto") << "cannot book: " << name_ << "\n" << conf_.dump();
0154       throw;
0155     }
0156 
0157     //cosmetics
0158     h_->GetXaxis()->SetTitle(xAxis.Label().c_str());
0159     h_->SetYTitle(yLabel.c_str());
0160 
0161     if (conf_.exists("weight")) {
0162       w_ = edm::Service<VariableHelperService>()->get().variable(conf_.getParameter<std::string>("weight"));
0163     }
0164   }
0165 
0166   virtual void fill(const edm::Event& iEvent) {
0167     if (!h_)
0168       return;
0169 
0170     double weight = 1.0;
0171     if (w_) {
0172       if (w_->compute(iEvent))
0173         weight = (*w_)(iEvent);
0174       else {
0175         edm::LogInfo("ConfigurableHisto") << "could not compute the weight for: " << name_ << " with config:\n"
0176                                           << conf_.dump() << " default to 1.0";
0177       }
0178     }
0179 
0180     TProfile* pcast(nullptr);
0181     TH2* h2cast(nullptr);
0182     switch (type_) {
0183       case h1:
0184         if (!h_)
0185           throw;
0186         if (x_->compute(iEvent))
0187           h_->Fill((*x_)(iEvent), weight);
0188         else {
0189           edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
0190         }
0191         break;
0192       case prof:
0193         pcast = dynamic_cast<TProfile*>(h_);
0194         if (!pcast)
0195           throw;
0196         if (x_->compute(iEvent) && y_->compute(iEvent))
0197           pcast->Fill((*x_)(iEvent), (*y_)(iEvent), weight);
0198         else {
0199           edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
0200         }
0201         break;
0202       case h2:
0203         h2cast = dynamic_cast<TH2*>(h_);
0204         if (!h2cast)
0205           throw;
0206         if (x_->compute(iEvent) && y_->compute(iEvent))
0207           h2cast->Fill((*x_)(iEvent), (*y_)(iEvent), weight);
0208         else {
0209           edm::LogInfo("ConfigurableHisto") << "could not fill: " << name_ << " with config:\n" << conf_.dump();
0210         }
0211         break;
0212     }
0213   }
0214   const HType& type() { return type_; }
0215 
0216   void complete() {}
0217   TH1* h() { return h_; }
0218 
0219 protected:
0220   ConfigurableHisto(const ConfigurableHisto& master) {
0221     type_ = master.type_;
0222     h_ = nullptr;  //no histogram attached in copy constructor
0223     name_ = master.name_;
0224     conf_ = master.conf_;
0225     x_ = master.x_;
0226     y_ = master.y_;
0227     z_ = master.z_;
0228     w_ = master.w_;
0229   }
0230   HType type_;
0231   TH1* h_;
0232   std::string name_;
0233   edm::ParameterSet conf_;
0234 
0235   const CachingVariable* x_;
0236   const CachingVariable* y_;
0237   const CachingVariable* z_;
0238   const CachingVariable* w_;
0239 };
0240 
0241 class SplittingConfigurableHisto : public ConfigurableHisto {
0242 public:
0243   SplittingConfigurableHisto(HType t, std::string name, edm::ParameterSet& pset)
0244       : ConfigurableHisto(t, name, pset), splitter_(nullptr) {
0245     std::string title = pset.getParameter<std::string>("title");
0246 
0247     //allow for many splitters ...
0248     if (pset.exists("splitters")) {
0249       //you want more than one splitter
0250       std::vector<std::string> splitters = pset.getParameter<std::vector<std::string> >("splitters");
0251       for (unsigned int s = 0; s != splitters.size(); ++s) {
0252         const CachingVariable* v = edm::Service<VariableHelperService>()->get().variable(splitters[s]);
0253         const Splitter* splitter = dynamic_cast<const Splitter*>(v);
0254         if (!splitter) {
0255           edm::LogError("SplittingConfigurableHisto")
0256               << "for: " << name_ << " the splitting variable: " << splitters[s] << " is not a Splitter";
0257           continue;
0258         }
0259 
0260         //insert in the map
0261         std::vector<ConfigurableHisto*>& insertedHisto = subHistoMap_[splitter];
0262         //now configure the histograms
0263         unsigned int mSlots = splitter->maxSlots();
0264         for (unsigned int i = 0; i != mSlots; ++i) {
0265           //---   std::cout<<" slot: "<<i<<std::endl;
0266           const std::string& slabel = splitter->shortLabel(i);
0267           const std::string& label = splitter->label(i);
0268           edm::ParameterSet mPset = pset;
0269           edm::Entry e("string", title + " for " + label, true);
0270           mPset.insert(true, "title", e);
0271           insertedHisto.push_back(new ConfigurableHisto(t, name + slabel, mPset));
0272         }  //loop on slots
0273       }    //loop on splitters
0274 
0275     }  //if splitters exists
0276     else {
0277       //single splitter
0278       //get the splitting variable
0279       const CachingVariable* v =
0280           edm::Service<VariableHelperService>()->get().variable(pset.getParameter<std::string>("splitter"));
0281       splitter_ = dynamic_cast<const Splitter*>(v);
0282       if (!splitter_) {
0283         edm::LogError("SplittingConfigurableHisto")
0284             << "for: " << name_ << " the splitting variable: " << v->name() << " is not a Splitter";
0285       } else {
0286         //configure the splitted plots
0287         unsigned int mSlots = splitter_->maxSlots();
0288         for (unsigned int i = 0; i != mSlots; i++) {
0289           const std::string& slabel = splitter_->shortLabel(i);
0290           const std::string& label = splitter_->label(i);
0291           edm::ParameterSet mPset = pset;
0292           edm::Entry e("string", title + " for " + label, true);
0293           mPset.insert(true, "title", e);
0294           subHistos_.push_back(new ConfigurableHisto(t, name + slabel, mPset));
0295         }
0296       }
0297     }
0298   }  //end of ctor
0299 
0300   void book(TFileDirectory* dir) override {
0301     //book the base histogram
0302     ConfigurableHisto::book(dir);
0303 
0304     if (!subHistoMap_.empty()) {
0305       SubHistoMap::iterator i = subHistoMap_.begin();
0306       SubHistoMap::iterator i_end = subHistoMap_.end();
0307       for (; i != i_end; ++i) {
0308         for (unsigned int h = 0; h != i->second.size(); ++h) {
0309           i->second[h]->book(dir);
0310         }
0311         //book the THStack
0312         std::string sName = name_ + "_" + i->first->name();
0313         std::string sTitle = "Stack histogram of " + name_ + " for splitter " + i->first->name();
0314         subHistoStacks_[i->first] = dir->make<THStack>(sName.c_str(), sTitle.c_str());
0315       }
0316     } else {
0317       for (unsigned int h = 0; h != subHistos_.size(); h++) {
0318         subHistos_[h]->book(dir);
0319       }
0320       //book a THStack
0321       std::string sName = name_ + "_" + splitter_->name();
0322       std::string sTitle = "Stack histogram of " + name_ + " for splitter " + splitter_->name();
0323       stack_ = dir->make<THStack>(sName.c_str(), sTitle.c_str());
0324     }
0325   }
0326 
0327   ConfigurableHisto* clone() const override { return new SplittingConfigurableHisto(*this); }
0328 
0329   void fill(const edm::Event& e) override {
0330     //fill the base histogram
0331     ConfigurableHisto::fill(e);
0332 
0333     if (!subHistoMap_.empty()) {
0334       SubHistoMap::iterator i = subHistoMap_.begin();
0335       SubHistoMap::iterator i_end = subHistoMap_.end();
0336       for (; i != i_end; ++i) {
0337         const Splitter* splitter = i->first;
0338         if (!splitter)
0339           continue;
0340         if (!splitter->compute(e))
0341           continue;
0342         unsigned int slot = (unsigned int)(*splitter)(e);
0343         if (slot >= i->second.size()) {
0344           edm::LogError("SplittingConfigurableHisto")
0345               << "slot index: " << slot << " is bigger than slots size: " << i->second.size()
0346               << " from variable value: " << (*splitter)(e);
0347           continue;
0348         }
0349         //fill in the proper slot
0350         i->second[slot]->fill(e);
0351       }
0352     } else {
0353       //fill the component histograms
0354       if (!splitter_)
0355         return;
0356       if (!splitter_->compute(e)) {
0357         return;
0358       }
0359       unsigned int slot = (unsigned int)(*splitter_)(e);
0360       if (slot >= subHistos_.size()) {
0361         edm::LogError("SplittingConfigurableHisto")
0362             << "slot index: " << slot << " is bigger than slots size: " << subHistos_.size()
0363             << " from variable value: " << (*splitter_)(e);
0364         return;
0365       }
0366       subHistos_[slot]->fill(e);
0367     }
0368   }
0369 
0370   void complete() {
0371     if (!subHistoMap_.empty()) {
0372       //fill up the stacks
0373       SubHistoMap::iterator i = subHistoMap_.begin();
0374       SubHistoMap::iterator i_end = subHistoMap_.end();
0375       for (; i != i_end; ++i) {
0376         for (unsigned int h = 0; h != i->second.size(); h++) {
0377           //      if (i->second[h]->h()->Integral==0) continue;// do not add empty histograms. NO, because it will be tough to merge two THStack
0378           subHistoStacks_[i->first]->Add(i->second[h]->h(), i->first->label(h).c_str());
0379         }
0380       }
0381 
0382     } else {
0383       //fill up the only stack
0384       for (unsigned int i = 0; i != subHistos_.size(); i++) {
0385         stack_->Add(subHistos_[i]->h(), splitter_->label(i).c_str());
0386       }
0387     }
0388   }
0389 
0390 private:
0391   SplittingConfigurableHisto(const SplittingConfigurableHisto& master) : ConfigurableHisto(master) {
0392     splitter_ = master.splitter_;
0393     if (!master.subHistoMap_.empty()) {
0394       SubHistoMap::const_iterator i = master.subHistoMap_.begin();
0395       SubHistoMap::const_iterator i_end = master.subHistoMap_.end();
0396       for (; i != i_end; ++i) {
0397         const std::vector<ConfigurableHisto*>& masterHistos = i->second;
0398         std::vector<ConfigurableHisto*>& clonedHistos = subHistoMap_[i->first];
0399         for (unsigned int i = 0; i != masterHistos.size(); i++) {
0400           clonedHistos.push_back(masterHistos[i]->clone());
0401         }
0402       }
0403     } else {
0404       for (unsigned int i = 0; i != master.subHistos_.size(); i++) {
0405         subHistos_.push_back(master.subHistos_[i]->clone());
0406       }
0407     }
0408   }
0409 
0410   typedef std::map<const Splitter*, std::vector<ConfigurableHisto*> > SubHistoMap;
0411   typedef std::map<const Splitter*, THStack*> SubHistoStacks;
0412   SubHistoStacks subHistoStacks_;
0413   SubHistoMap subHistoMap_;
0414 
0415   const Splitter* splitter_;
0416   std::vector<ConfigurableHisto*> subHistos_;
0417   // do you want to have a stack already made from the splitted histos?
0418   THStack* stack_;
0419 };
0420 
0421 #endif