StatSample.cpp

Go to the documentation of this file.

#include "cafe/StatSample.hpp"
#include <string>

using namespace std ;

namespace cafe {

  StatSample::StatSample(const string& sample) : _name(sample)
  {
    _events.push_back(new StatSelection("Initial",0)) ;
    _weight = new StatWeight("Efficiency correction") ;
  }

  StatSample::~StatSample() {
  }

  void StatSample::AddTags(const vector<string>& tags) {
    for (vector<string>::const_iterator it = tags.begin() ;
         it != tags.end(); it++) {
      string s = *it ;
      while (s.size() > 0 && s.substr(0,1) == " ") s.erase(0,1) ;
      while (s.size() > 0 && s.substr(s.size()-1,1) == " ") s.erase(s.size()-1,1) ;      
      if (s.size() ==0 ) continue ;      
      _tags.push_back(s) ;
    }
    return ;
  }


  void StatSample::AddAndTags(const vector<string>& tags) {
    for (vector<string>::const_iterator it = tags.begin() ;
         it != tags.end(); it++) {
      string s = *it ;
      while (s.size() > 0 && s.substr(0,1) == " ") s.erase(0,1) ;
      while (s.size() > 0 && s.substr(s.size()-1,1) == " ") s.erase(s.size()-1,1) ;      
      if (s.size() ==0 ) continue ;      
      _tags_and.push_back(s) ;
    }
    return ;
  }

  unsigned long StatSample::add(cafe::Event* event, const std::string& name) {
    if (!event) {
      cerr << "StatSample ERROR!. Event pointer == 0." << endl ;
      return 0;
    }    

    if (!tagged(event)) return 0 ;

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      _events.push_back(new StatSelection(name)) ;
      it = _events.end()-1 ;
    } 

    return (*it)->addEvent() ;
  }



  double StatSample::applyWeight(cafe::Event* event, const std::string& name, 
                                 double weight, double weight_pos, double weight_neg) {    
    if (!event) {
      cerr << "StatSample ERROR!. Event pointer == 0." << endl ;
      return 0;
    }    

    if (!tagged(event)) return 0.0 ;

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      _events.push_back(new StatWeight(name)) ;
      it = _events.end()-1 ;
    }

    // update global weight
    _weight->applyWeight(weight, weight_pos, weight_neg) ;

    // update specified weight
    return dynamic_cast<StatWeight*>(*it) ->applyWeight(weight, weight_pos, weight_neg) ;
  }

  double StatSample::calculateGlobalWeight(cafe::Event* event) {

    if (!event) {
      cerr << "StatSample ERROR!. Event pointer == 0." << endl ;
      return 0;
    }    

    
    if (!tagged(event)) return -1.0 ;
    
    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == _weight->name() ) break ;

    // create the global weight record, if it is not existing 
    // the global weight itself exist always
    if (it == _events.end()) {
      _events.push_back(new StatWeight(_weight->name())) ;
      it = _events.end()-1 ;
    }



    for(vector<StatSelection*>::const_iterator jt = _events.begin(); 
        jt != it; jt++) {
      if ((*jt)->isWeight()) {
        StatWeight* w = dynamic_cast<StatWeight*>(*jt) ;

        // move global weight calculation to the applyWeight
        //      weight *=  w->weight() ; // global weight calculation

        // systematics
        if (w->weight_average_pos() >= 0 || w->weight_average_neg() >= 0 ) {
          vector<pair<StatWeight*,StatWeight*> >::iterator wt = _syst.begin(); 
          for(; wt != _syst.end(); wt++) 
            if (wt->first->name().find(w->name()) != string::npos) break ;

          if(wt == _syst.end()) {
            _syst.push_back(pair<StatWeight*,StatWeight*>(new StatWeight(w->name()),
                                                          new StatWeight(w->name()))) ;
            wt = _syst.end()-1 ;
          }

          for(vector<StatSelection*>::const_iterator jtt = _events.begin(); 
              jtt != it; jtt++) {
            if ((*jtt)->isWeight()) {
              StatWeight* ww = dynamic_cast<StatWeight*>(*jtt) ;
              if (jtt == jt) {
                wt->first->applyWeight(ww->weight_pos())  ;
                wt->second->applyWeight(ww->weight_neg()) ;
              } else {
                wt->first->applyWeight(ww->weight()) ;
                wt->second->applyWeight(ww->weight()) ;
              }
            }
          }
        }
      }
    }

    return (dynamic_cast<StatWeight*>(*it))->applyWeight(_weight->weight(),_weight->weight_pos(),_weight->weight_neg()) ;
  }


   unsigned long StatSample::nevents(unsigned int n ) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }
    return _events[n]->nevents() ;
  }


   unsigned long StatSample::nevents(const string& name) const {

     vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample: \"" << _name << "\": No weight with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return 0 ;
    }
    return (*it)->nevents() ;
  }

  void StatSample::Clear() {
    for (vector<StatSelection*>::iterator it = _events.begin() ;
         it != _events.end(); it++) 
      (*it)->Clear() ;

    for (vector<pair<StatWeight*,StatWeight*> >::iterator it = _syst.begin() ;
         it != _syst.end(); it++) {
      it->first->Clear() ;
      it->second->Clear() ;
    }

    _weight->Clear() ;
    return ;
  }

  const StatWeight* StatSample::eventWeight(const std::string& name) const {

    if (name == "global") return _weight ;

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample: \"" << _name << "\": No weight with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return 0 ;
    }

    if (!(*it)->isWeight()) {
      cout << "cafe::StatSample: \"" << _name << "\": Selection with name \"" << name 
           << "\" is not the event weight!" << endl ; 
      return 0 ;
    }
    return dynamic_cast<StatWeight*> (*it) ;
  }

  const StatWeight* StatSample::eventWeight(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample: \"" << _name << "\": weight number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }
    if (!_events[n]->isWeight()) {
      cout << "cafe::StatSample \"" << _name << "\": Selection with number \"" << n
           << "\" is not the event weight!" << endl ; 
      return 0 ;
    }
    return dynamic_cast<StatWeight*> (_events[n]);
  }

  const StatSelection* StatSample::eventSelection(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": Selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }
    return  _events[n];
  }

  const pair<StatWeight*,StatWeight*>& StatSample::systematics(unsigned int n) const {
    if (n >= _syst.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": Systematics number " << n 
           << " greater than systematics size " << _syst.size() << endl ;
      exit(1) ;
    }
    return  _syst[n];
  }


  const StatSelection* StatSample::eventSelection(const string& name) const {
    
    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample \"" << _name << "\": No selection with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return 0 ;
    }

    return (*it) ;
  }

  double StatSample::eff(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSamplec\"" << _name << "\": Selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }

    if (_events[n]->isWeight()) return dynamic_cast<StatWeight*> (_events[n])->weight_average() ;

    if ( n > 0) 
      return ((double)_events[n]->nevents())/_events[n-1]->nevents() ;

    return -1.0 ;
  }


  double StatSample::eff(const std::string& name) const {

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample \"" << _name << "\": No selection with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return 0 ;
    }

    if (it == _events.begin()) return 0 ;

    return ((double) (*it)->nevents())/(*(it-1))->nevents() ;
  }


  double StatSample::effErr(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": Selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }

    if (_events[n]->isWeight()) return dynamic_cast<StatWeight*> (_events[n])->err() ;

    if ( n >0) 
      return _events[n]->err(*_events[n-1]) ;

    return -1.0 ;
  }

  double StatSample::effErr(const std::string& name) const {

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample \"" << _name << "\": No selection with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return 0 ;
    }

    if (it == _events.begin()) return 0 ;

    return (*it)->err(**(it-1)) ;
  }


  double StatSample::effGlob(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": Selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }

    if ( _events.front()->nevents() == 0 ) return -1.0 ; 

    if (n == 0) 
      return ((double)_events.back()->nevents())/_events.front()->nevents() ;
    
    return ((double)_events[n]->nevents())/_events.front()->nevents() ;
  }


  double StatSample::effGlob(const std::string& name) const {

    if ( _events.front()->nevents() == 0 ) return -1.0 ; 

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample \"" << _name << "\": No selection with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return -1.0 ;
    }

    if (it == _events.begin()) return -1.0 ;

    return ((double) (*it)->nevents())/_events.front()->nevents() ;
  }


  double StatSample::effErrGlob(unsigned int n) const {
    if (n >= _events.size()) {
      cerr << "ERROR! StatSample \"" << _name << "\": Selection number " << n 
           << " greater than selection size " << _events.size() << endl ;
      exit(1) ;
    }

    if ( _events.front()->nevents() == 0)  return -1.0 ;
    
    if ( n == 0) 
      return _events.back()->err(*_events.front()); 

    return _events[n]->err(*_events.front());
  }

  double StatSample::effErrGlob(const std::string& name) const {

    if ( _events.front()->nevents() == 0 ) return -1.0 ; 

    vector<StatSelection*>::const_iterator it = _events.begin();
    for (; it != _events.end(); it++) if ((*it)->name() == name ) break ;

    if (it == _events.end()) {
      cout << "cafe::StatSample \"" << _name << "\": No selection with name \"" << name 
           << "\" found in the Stat class" << endl ; 
      return -1.0 ;
    }

    if (it == _events.begin()) return -1.0 ;

    return (*it)->err(*_events.front());
  }

  bool StatSample::compareNames(const StatSample& sample) const {
    
    if (_events.size() != sample._events.size()) return false ;

    for (unsigned int n = 0 ; n < _events.size() ; n++)
      if (_events[n]->name() != sample._events[n]->name() ) return false ;

    return true ;
  }

  double StatSample::correctedEff(const StatWeight* w) const {
    if(!w) {
      vector<StatSelection*>::const_iterator it = _events.begin();
      for (; it != _events.end(); it++) if ((*it)->name() == _weight->name() ) break ;
      if (it == _events.end()) {
        cout << "cafe::StatSample ERROR: The global weight processor is not used! Add cafe.Run: StatGWeight. " 
             << "Will return non-corrected efficiency" << endl ;
        return effGlob() ;
      }
      w = dynamic_cast<StatWeight*>(*it) ;
    }
    return effGlob()*w->weight_average() ;
  }

  double StatSample::correctedEffErr(const StatWeight* w) const {
    if(!w) w = _weight ;    
    return sqrt(pow(w->weight_average()*effErrGlob(),2) 
                + pow(effGlob()*w->err(),2)) ;
  }

  double StatSample::totalsyst_pos() const {
    double w = _weight->weight_average() ;
    if (w<=0) return 0.0 ;
    double s = 0.0 ;
    for (std::vector<std::pair<StatWeight*,StatWeight*> >::const_iterator it = _syst.begin() ;
         it != _syst.end(); it++) 
      s += pow(it->first->weight_average() / w - 1.0 , 2) ;    
    return sqrt(s) ;
  }

  double StatSample::totalsyst_neg() const {
    double w = _weight->weight_average() ;
    if (w<=0) return 0.0 ;
    double s = 0.0 ;
    for (std::vector<std::pair<StatWeight*,StatWeight*> >::const_iterator it = _syst.begin() ;
         it != _syst.end(); it++) 
      s += pow(it->second->weight_average() / w - 1.0 , 2) ;    
    return sqrt(s) ;
  }

  bool StatSample::tagged(const cafe::Event* event) const {

    if (_tags.size()  == 0 &&  _tags_and.size()  == 0 ) return true ;

    for (vector<string>::const_iterator it = _tags_and.begin() ; 
         it != _tags_and.end(); it++) {
       if (!event->hasTag(*it)) return false ;
    }

    return _tags.size()  == 0 || event->hasTag(_tags) ; 
  }


  ostream& StatSample::HtmlTable(ostream& os) const {

    os << "<TABLE BORDERCOLORDARK=\"996633\"  BORDERCOLOR=\"CC9966\"    BORDERCOLORLIGHT=\"FFCC99\"   border=3 cellPadding=5 width=\"99%\">" << endl ;
    os << "<TBODY align=right>" << endl ;
    
    os << "<TR>" << endl ;
    os << "<TH align=left>" << endl ;
    os << "SELECTION"   ;
    os << " </TH>" << endl ;
    os << "<TH align=center colspan=3>" << endl ; 
    os <<  name() << endl ;
    os << "</TH>" << endl ; 
    os << "</TR>" << endl ;
      
    bool color = false ;
    
    for (unsigned int n = 0; n < size() ; n++) {
      
      if (color) {
        os << "<TR BGCOLOR=#bfefff>" << endl ; 
        color = !color ;
      } else { 
        os << "<TR>" << endl ; 
        color = !color ;
      }      
     
      os << "<TD align=left>" << endl ; 
      os <<  eventSelection(n)->name() << endl ;
      os << "</TD>" << endl ; 

      os << "<TD TITLE=\"Number of events\">" << endl ; 
      os <<  nevents(n) << endl ;
      if (n==0) {
        os << "<TD>&#160;</TD> <TD>&#160;</TD>" << endl ; 
        continue ;
      }
        
      if (eventSelection(n)->isWeight()) {
          if (eventSelection(n)->name() == eventWeight()->name()) 
            os << "</TD> <TD TITLE=\"Global event weight\">" << endl ; 
          else 
            os << "</TD> <TD TITLE=\"Average weight\">" << endl ; 
          os << eff(n) << "&#177;" << effErr(n)<< endl ;
          os << "</TD> <TD> &#160;</TD>"<< endl ;
          continue ;
      } 
      
      os << "</TD><TD TITLE=\"Selection efficiency\">"<< endl ;
      os << 100.0*eff(n) << "&#177;" << 100.0*effErr(n) << "%"
         << endl ;
      os << "</TD><TD TITLE=\"Global efficiency\">"<< endl ;
      os << 100.0*effGlob(n) << "&#177;" << 100.0*effErrGlob(n) << "%" 
         << endl ;
      os << "</TD>" << endl ;       
      os << "</TR>" << endl ; 
    }
    
    if(eventWeight()->nevents() > 0) {
      
      os << "<TR BGCOLOR=#FFCCCC>" << endl ; 
      os << "<TD align=left>Corrected Efficiency</TD>" << endl ; 
      os << "<TD TITLE=\"Number Of Events\">"<< endl ;
      os <<  nevents(size()-1) << endl ;
      os << "</TD>&#160;<TD>" << endl ; 
      os << "</TD><TD TITLE=\"Global efficiency\">"<< endl ;
      os << 100.0*correctedEff() << "&#177;" 
         << 100.0*correctedEffErr()<< "%" 
         << endl ;
      os << "</TD>"<< endl ;
      os << "</TR>" << endl ; 
    }
    os << "</TBODY></TABLE>" << endl ;

    return os ;
  }


  ostream& StatSample::print_tex (ostream& os, const string title) const {
    
    os << "" << endl;
    os << "% Efficiency table %" << endl;
    os << "" << endl;
    os << "\\begin{table}[p]" << endl;
    os << "\\begin{center}" << endl;
    os << "\\begin{tabular}{lrrr}" << endl;
    os << "  \\hline  \\hline" << endl;
    os << "  Selection & Events & Relative  & Total  \\\\ \\hline" << endl;
    
    os << tex(eventSelection(0)->name()) << " & " ;
    os <<  nevents(0) << " & & \\\\ " << endl ; 

    for (unsigned int n = 1; n < size() ; n++) {
      
      os << tex(eventSelection(n)->name()) << " & " ;
      os <<  nevents(n) << " & ";

      if (eventSelection(n)->isWeight()) {      
        os << " $\\mathit{ " << eff(n) << " \\pm " 
           << effErr(n) << " }$ & " ;
      } else {  
        os << " $ " << 100.0*eff(n) << " \\pm " 
           << 100.0*effErr(n) << " ~\\% $ & " ;
        os << " $ " << 100.0*effGlob(n) << " \\pm " 
           << 100.0*effErrGlob(n) << "~\\% $ " ;
      }

      os << " \\\\ " << endl ; 
    }

    os << "  \\hline " << endl;
    const StatWeight* w = eventWeight() ;
    if (w->nevents() > 0) {
      os << " Corrected efficiency & " ;
      os <<  nevents(size()-1) << " & ";
      
      os << "& $" << 100.0*correctedEff()<< " \\pm " 
         << 100.0*correctedEffErr() << "~\\% $ \\\\"  ;
      
      os << " \\hline"<< endl ;
    }
    os << " \\hline" << endl;
    os << " \\end{tabular}" << endl;
    os << " \\caption{" << title << " Sample " << name() << " }" << endl;
    os << " \\label{Table:" <<  title.c_str() << ":" << name() << "}" << endl;
    os << "\\end{center}" << endl;
    os << "\\end{table}" << endl; 
    os << endl;

    return os ;    
  }

  
  ostream& StatSample::print_tex_syst (ostream& os) const {
    
    if (_weight->weight_average() <= 0) return os;
      
    os << "" << endl;
    os << "% Systematics table %" << endl;
    os << "" << endl;
    os << "\\begin{table}[p]" << endl;
    os << "\\begin{center}" << endl;
    os << "\\begin{tabular}{lrr}" << endl;
    os << "  \\hline  \\hline" << endl;
    os << "  Systematics & Positive [\\%] & Negative  [\\%] \\\\ \\hline" << endl;
    
    
    for (vector<pair<StatWeight*,StatWeight*> >::const_iterator jt = _syst.begin();
         jt != _syst.end() ; jt++) {
      os << tex(jt->first->name()) << " &  " ;
      double pos = jt->first->weight_average() / _weight->weight_average() - 1.0 ;
      double neg = jt->second->weight_average() / _weight->weight_average() - 1.0 ;
      os << " $ " << 100.0*pos << " $ & $ " <<  100.0*neg << " $ \\\\ " << endl ;
    }
    
    os << "  \\hline " << endl;
    
    os << " Total : & $ " << 100.0*totalsyst_pos() << " $ & $ " ;
    os << -100.0*totalsyst_neg() << " $ \\\\ " ;      
    os << " \\hline"<< endl ;
    
    os << " \\hline" << endl;
    os << " \\end{tabular}" << endl;
    os << " \\caption{ Systematics for the sample " << name() << " }" << endl;
    os << " \\label{Table:syst:" << name() << "}" << endl;
    os << "\\end{center}" << endl;
    os << "\\end{table}" << endl; 
    os << endl;

    return os ;
  }

  std::string StatSample::tex(const std::string& init) {
    string name ;
    for (unsigned int i=0; i<init.size(); i++) {
      string ch ;
      if (init.substr(i,1) == "_") ch = "\\_" ;
      else if (init.substr(i,2) == ">=") {
        ch = "$\\geq$" ;
        i++ ;
      }
      else if (init.substr(i,2) == "<=") {
        ch = "$\\leq$" ;
        i++ ;
      }
      else if (init.substr(i,1) == "<") ch = "$<$" ;
      else if (init.substr(i,1) == ">") ch = "$>$" ;
      else ch = init.substr(i,1) ;
      name = name + ch ;
    }
    return name ;
  }

}

ClassImp(cafe::StatSample);

---