CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/PhysicsTools/UtilAlgos/interface/ConfigurableHisto.h

Go to the documentation of this file.
00001 #ifndef ConfigurableAnalysis_ConfigurableHisto_H
00002 #define ConfigurableAnalysis_ConfigurableHisto_H
00003 
00004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00005 #include "CommonTools/Utils/interface/TFileDirectory.h"
00006 
00007 #include "TH1.h"
00008 #include "TH1F.h"
00009 #include "TH2F.h"
00010 #include "TProfile.h"
00011 #include "THStack.h"
00012 
00013 #include "PhysicsTools/UtilAlgos/interface/VariableHelper.h"
00014 #include "PhysicsTools/UtilAlgos/interface/CachingVariable.h"
00015 
00016 
00017 class ConfigurableAxis {
00018  public:
00019   ConfigurableAxis(){}
00020   ConfigurableAxis(edm::ParameterSet par) :
00021     nBin_(0),Min_(0),Max_(0),Label_(""){
00022     Label_=par.getParameter<std::string>("Label");
00023 
00024     if (par.exists("nBins")){
00025       nBin_=par.getParameter<unsigned int>("nBins");
00026       Min_=par.getParameter<double>("Min");
00027       Max_=par.getParameter<double>("Max");
00028     }else{
00029       if (par.exists("vBins"))
00030         vBins_=par.getParameter<std::vector<double> >("vBins");
00031       else{
00032         Min_=par.getParameter<double>("Min");
00033         Max_=par.getParameter<double>("Max");
00034       }
00035     }
00036   }
00037   bool variableSize(){return vBins_.size()!=0;}
00038   unsigned int nBin(){if (variableSize()) return vBins_.size()-1;else return nBin_;}
00039   double Min(){if (variableSize()) return vBins_.front(); else return Min_;}
00040   double Max(){if (variableSize()) return vBins_.back(); else return Max_;}
00041   const std::string & Label(){ return Label_;}
00042   const double * xBins(){ if (vBins_.size()!=0) return &(vBins_.front()); else return 0;}
00043 
00044  private:
00045   std::vector<double> vBins_;
00046   unsigned int nBin_;
00047   double Min_;
00048   double Max_;
00049   std::string Label_;
00050 };
00051 
00052 
00053 class ConfigurableHisto {
00054  public:
00055   enum HType { h1 ,h2, prof };
00056   ConfigurableHisto(HType t, std::string name, edm::ParameterSet & iConfig) :
00057     type_(t),h_(0),name_(name), conf_(iConfig),x_(0),y_(0),z_(0),w_(0){}
00058 
00059   virtual ~ConfigurableHisto(){}
00060 
00061   virtual ConfigurableHisto * clone() const { return new ConfigurableHisto(*this);}
00062 
00063   virtual void book(TFileDirectory* dir){
00064     std::string title=conf_.getParameter<std::string>("title");
00065     edm::ParameterSet xAxisPSet=conf_.getParameter<edm::ParameterSet>("xAxis");
00066     ConfigurableAxis xAxis(xAxisPSet);
00067     x_=edm::Service<VariableHelperService>()->get().variable(xAxisPSet.getParameter<std::string>("var"));
00068 
00069     std::string yLabel="";    
00070     bool yVBin=false;
00071     ConfigurableAxis yAxis;
00072     if (conf_.exists("yAxis")){
00073       edm::ParameterSet yAxisPSet=conf_.getParameter<edm::ParameterSet>("yAxis");
00074       yAxis=ConfigurableAxis(yAxisPSet);
00075       yLabel=yAxis.Label();
00076       //at least TH2 or TProfile
00077       if (yAxisPSet.exists("var"))
00078         y_=edm::Service<VariableHelperService>()->get().variable(yAxisPSet.getParameter<std::string>("var"));
00079       yVBin=yAxis.variableSize();
00080     }
00081     
00082     if (conf_.exists("zAxis")){
00083       throw;
00084     }
00085 
00086     bool xVBin=xAxis.variableSize();
00087 
00088 
00089     if (type()==h1){
00090       //make TH1F
00091       if (xVBin)
00092         h_=dir->make<TH1F>(name_.c_str(),title.c_str(),
00093                            xAxis.nBin(), xAxis.xBins());
00094       else
00095         h_=dir->make<TH1F>(name_.c_str(),title.c_str(),
00096                            xAxis.nBin(),xAxis.Min(),xAxis.Max());
00097     }
00098     else if (type()==h2){
00099       //make TH2F
00100       if (xVBin){
00101         if (yVBin)
00102           h_=dir->make<TH2F>(name_.c_str(),title.c_str(),
00103                              xAxis.nBin(),xAxis.xBins(),
00104                              yAxis.nBin(),yAxis.xBins());
00105         else
00106           h_=dir->make<TH2F>(name_.c_str(),title.c_str(),
00107                              xAxis.nBin(),xAxis.xBins(),
00108                              yAxis.nBin(),yAxis.Min(),yAxis.Max());
00109       }else{
00110         if (yVBin)
00111           h_=dir->make<TH2F>(name_.c_str(),title.c_str(),
00112                              xAxis.nBin(),xAxis.Min(),xAxis.Max(),
00113                              yAxis.nBin(),yAxis.xBins());
00114         else
00115           h_=dir->make<TH2F>(name_.c_str(),title.c_str(),
00116                              xAxis.nBin(),xAxis.Min(),xAxis.Max(),
00117                              yAxis.nBin(),yAxis.Min(),yAxis.Max());
00118       }
00119     }
00120     else if (type()==prof){
00121       //make TProfile
00122       std::string pFopt="";
00123       if (conf_.exists("Option"))
00124         pFopt=conf_.getParameter<std::string>("Option");
00125       
00126       if (xVBin)
00127       h_=dir->make<TProfile>(name_.c_str(),title.c_str(),
00128                              xAxis.nBin(), xAxis.xBins(),
00129                              yAxis.Min(),yAxis.Max(),
00130                              pFopt.c_str());
00131       else
00132         h_=dir->make<TProfile>(name_.c_str(),title.c_str(),
00133                                xAxis.nBin(),xAxis.Min(),xAxis.Max(),
00134                                yAxis.Min(),yAxis.Max(),
00135                                pFopt.c_str());
00136       
00137     }
00138     else {
00139       edm::LogError("ConfigurableHisto")<<"cannot book: "<<name_<<"\n"<<conf_.dump();
00140       throw;
00141     }
00142 
00143     //cosmetics
00144     h_->GetXaxis()->SetTitle(xAxis.Label().c_str());
00145     h_->SetYTitle(yLabel.c_str());
00146     
00147     if (conf_.exists("weight"))
00148       {
00149         w_=edm::Service<VariableHelperService>()->get().variable(conf_.getParameter<std::string>("weight"));
00150       }
00151   }
00152   
00153   virtual void fill(const edm::Event &iEvent) {
00154     if (!h_)      return;
00155 
00156     double weight=1.0;
00157     if (w_){
00158       if (w_->compute(iEvent))
00159         weight=(*w_)(iEvent);
00160       else{
00161         edm::LogInfo("ConfigurableHisto")<<"could not compute the weight for: "<<name_
00162                                          <<" with config:\n"<<conf_.dump()
00163                                          <<" default to 1.0";
00164       }
00165     }
00166     
00167     TProfile * pcast(0);
00168     TH2 * h2cast(0);
00169     switch(type_){
00170     case h1:
00171       if (!h_) throw;
00172       if (x_->compute(iEvent)) h_->Fill((*x_)(iEvent),weight);
00173       else{
00174         edm::LogInfo("ConfigurableHisto")<<"could not fill: "<<name_
00175                                          <<" with config:\n"<<conf_.dump();
00176       }
00177       break;
00178     case prof:
00179       pcast=dynamic_cast<TProfile*>(h_);
00180       if (!pcast) throw;
00181       if (x_->compute(iEvent) && y_->compute(iEvent)) pcast->Fill((*x_)(iEvent),(*y_)(iEvent),weight);
00182       else{
00183         edm::LogInfo("ConfigurableHisto")<<"could not fill: "<<name_
00184                                          <<" with config:\n"<<conf_.dump();
00185       }
00186       break;
00187     case h2:
00188       h2cast=dynamic_cast<TH2*>(h_);
00189       if (!h2cast) throw;
00190       if (x_->compute(iEvent) && y_->compute(iEvent)) h2cast->Fill((*x_)(iEvent),(*y_)(iEvent),weight);
00191       else{
00192         edm::LogInfo("ConfigurableHisto")<<"could not fill: "<<name_
00193                                          <<" with config:\n"<<conf_.dump();
00194       }
00195       break;
00196     }
00197   }
00198   const HType & type() { return type_;}  
00199 
00200   void complete() {}
00201   TH1 * h() {return h_;}
00202 
00203  protected:
00204   ConfigurableHisto(const ConfigurableHisto & master){
00205     type_=master.type_;
00206     h_=0; //no histogram attached in copy constructor
00207     name_=master.name_;
00208     conf_=master.conf_;
00209     x_=master.x_;
00210     y_=master.y_;
00211     z_=master.z_;
00212     w_=master.w_;
00213   }
00214   HType type_;
00215   TH1 * h_;
00216   std::string name_;
00217   edm::ParameterSet conf_;
00218   
00219   const CachingVariable * x_;
00220   const CachingVariable * y_;
00221   const CachingVariable * z_;
00222   const CachingVariable * w_;
00223 };
00224 
00225 class SplittingConfigurableHisto : public ConfigurableHisto {
00226  public:
00227   SplittingConfigurableHisto(HType t, std::string name, edm::ParameterSet & pset) :
00228     ConfigurableHisto(t,name,pset) , splitter_(0) {
00229     std::string title=pset.getParameter<std::string>("title");
00230 
00231     //allow for many splitters ...
00232     if (pset.exists("splitters")){
00233       //you want more than one splitter
00234       std::vector<std::string> splitters = pset.getParameter<std::vector<std::string> >("splitters");
00235       for (unsigned int s=0;s!=splitters.size();++s){
00236         const CachingVariable * v=edm::Service<VariableHelperService>()->get().variable(splitters[s]);
00237         const Splitter * splitter = dynamic_cast<const Splitter*>(v);
00238         if (!splitter){
00239           edm::LogError("SplittingConfigurableHisto")<<"for: "<<name_<<" the splitting variable: "<<splitters[s]<<" is not a Splitter";
00240           continue;
00241         }
00242 
00243         //insert in the map
00244         std::vector<ConfigurableHisto*> & insertedHisto=subHistoMap_[splitter];
00245         //now configure the histograms
00246         unsigned int mSlots=splitter->maxSlots();
00247         for (unsigned int i=0;i!=mSlots;++i){
00248           //---   std::cout<<" slot: "<<i<<std::endl;
00249           const std::string & slabel=splitter->shortLabel(i);
00250           const std::string & label=splitter->label(i);
00251           edm::ParameterSet mPset=pset;
00252           edm::Entry e("string",title+" for "+label,true);
00253           mPset.insert(true,"title",e);
00254           insertedHisto.push_back(new ConfigurableHisto(t,name+slabel ,mPset));
00255         }//loop on slots
00256       }//loop on splitters
00257 
00258     }//if splitters exists
00259     else{
00260       //single splitter
00261       //get the splitting variable
00262       const CachingVariable * v=edm::Service<VariableHelperService>()->get().variable(pset.getParameter<std::string>("splitter"));
00263       splitter_ = dynamic_cast<const Splitter*>(v);
00264       if (!splitter_){
00265         edm::LogError("SplittingConfigurableHisto")<<"for: "<<name_<<" the splitting variable: "<<v->name()<<" is not a Splitter";
00266       }
00267       else{
00268         //configure the splitted plots
00269         unsigned int mSlots=splitter_->maxSlots();
00270         for (unsigned int i=0;i!=mSlots;i++){
00271           const std::string & slabel=splitter_->shortLabel(i);
00272           const std::string & label=splitter_->label(i);
00273           edm::ParameterSet mPset=pset;
00274           edm::Entry e("string",title+" for "+label,true);
00275           mPset.insert(true,"title",e);
00276           subHistos_.push_back(new ConfigurableHisto(t,name+slabel, mPset));
00277         }
00278       }
00279     }
00280   }//end of ctor
00281     
00282   void book(TFileDirectory* dir){
00283     //book the base histogram
00284     ConfigurableHisto::book(dir);
00285 
00286     if (subHistoMap_.size()!=0){
00287       SubHistoMap::iterator i=subHistoMap_.begin();
00288       SubHistoMap::iterator i_end=subHistoMap_.end();
00289       for (;i!=i_end;++i){for (unsigned int h=0;h!=i->second.size();++h){ 
00290           i->second[h]->book(dir);}
00291         //book the THStack
00292         std::string sName= name_+"_"+i->first->name();
00293         std::string sTitle="Stack histogram of "+name_+" for splitter "+i->first->name();
00294         subHistoStacks_[i->first]= dir->make<THStack>(sName.c_str(),sTitle.c_str());
00295       }
00296     }else{
00297       for (unsigned int h=0;h!=subHistos_.size();h++){subHistos_[h]->book(dir);}
00298       //book a THStack
00299       std::string sName= name_+"_"+splitter_->name();
00300       std::string sTitle="Stack histogram of "+name_+" for splitter "+splitter_->name();
00301       stack_ = dir->make<THStack>(sName.c_str(),sTitle.c_str());
00302     }
00303     
00304   }
00305 
00306   ConfigurableHisto * clone() const {    return new SplittingConfigurableHisto(*this);  }
00307 
00308   void fill(const edm::Event & e) {
00309     //fill the base histogram
00310     ConfigurableHisto::fill(e);
00311     
00312     if (subHistoMap_.size()!=0){
00313       SubHistoMap::iterator i=subHistoMap_.begin();
00314       SubHistoMap::iterator i_end=subHistoMap_.end();
00315       for (;i!=i_end;++i){
00316         const Splitter * splitter=i->first;
00317         if (!splitter) continue;
00318         if (!splitter->compute(e)) continue;
00319         unsigned int slot=(unsigned int)  (*splitter)(e);
00320         if (slot>=i->second.size()){
00321           edm::LogError("SplittingConfigurableHisto")<<"slot index: "<<slot<<" is bigger than slots size: "<<i->second.size()<<" from variable value: "<<(*splitter)(e);
00322           continue;}
00323         //fill in the proper slot
00324         i->second[slot]->fill(e);
00325       }
00326     }
00327     else{
00328       //fill the component histograms
00329       if (!splitter_) return;
00330       if (!splitter_->compute(e)){
00331         return;}
00332       unsigned int slot=(unsigned int) (*splitter_)(e);
00333       if (slot>=subHistos_.size()){
00334         edm::LogError("SplittingConfigurableHisto")<<"slot index: "<<slot<<" is bigger than slots size: "<< subHistos_.size() <<" from variable value: "<<(*splitter_)(e);
00335         return;}
00336       subHistos_[slot]->fill(e);
00337     }
00338   }
00339 
00340   void complete(){
00341     if (subHistoMap_.size()!=0){
00342       //fill up the stacks
00343       SubHistoMap::iterator i=subHistoMap_.begin();
00344       SubHistoMap::iterator i_end=subHistoMap_.end();
00345       for (;i!=i_end;++i){
00346         for (unsigned int h=0;h!=i->second.size();h++){
00347           //      if (i->second[h]->h()->Integral==0) continue;// do not add empty histograms. NO, because it will be tough to merge two THStack
00348           subHistoStacks_[i->first]->Add(i->second[h]->h(), i->first->label(h).c_str());
00349         }
00350       }
00351 
00352     }else{
00353       //fill up the only stack
00354       for (unsigned int i=0;i!=subHistos_.size();i++){  stack_->Add(subHistos_[i]->h(), splitter_->label(i).c_str());      }
00355     }
00356 
00357 
00358   }
00359  private:
00360   SplittingConfigurableHisto(const SplittingConfigurableHisto & master) : ConfigurableHisto(master){
00361     splitter_ = master.splitter_;
00362     if (master.subHistoMap_.size()!=0){
00363       SubHistoMap::const_iterator i=master.subHistoMap_.begin();
00364       SubHistoMap::const_iterator i_end=master.subHistoMap_.end();
00365       for (;i!=i_end;++i){
00366         const std::vector<ConfigurableHisto*> & masterHistos=i->second;
00367         std::vector<ConfigurableHisto*> & clonedHistos=subHistoMap_[i->first];
00368         for (unsigned int i=0;i!=masterHistos.size();i++){clonedHistos.push_back(masterHistos[i]->clone());}
00369       }
00370     }else{
00371       for (unsigned int i=0;i!=master.subHistos_.size();i++){subHistos_.push_back(master.subHistos_[i]->clone());}
00372     }
00373   }
00374 
00375   typedef std::map<const Splitter *, std::vector<ConfigurableHisto* > > SubHistoMap;
00376   typedef std::map<const Splitter *, THStack *> SubHistoStacks;
00377   SubHistoStacks subHistoStacks_;
00378   SubHistoMap subHistoMap_;
00379   
00380   const Splitter * splitter_;
00381   std::vector<ConfigurableHisto* > subHistos_;
00382   // do you want to have a stack already made from the splitted histos?
00383   THStack * stack_;
00384 };
00385 
00386 #endif