CAF: Stat.cpp Source File

00001 #include "cafe/Stat.hpp"
00002 #include "cafe/Config.hpp"
00003 #include "cafe/StatSample.hpp"
00004 #include <stdexcept>
00005 #include <fstream>
00006 #include <math.h>
00007 #include <algorithm>
00008 
00009 using namespace std;
00010 
00011 cafe::Stat* STAT = 0 ; 
00012 
00013 namespace cafe {
00014 
00015   StatPointer::StatPointer(Stat* stat) : _stat(stat) { }
00016   StatPointer::~StatPointer() {}
00017 
00018   void StatPointer::EventSelected(const string& selection_name) {
00019     if (_stat) _stat->EventSelected(selection_name)  ;
00020     return ;
00021   }
00022 
00023   void StatPointer::applyWeight(const std::string& weight_name, double weight, 
00024                                 double weight_pos, double weight_neg) {
00025     if (_stat) _stat->applyWeight(weight_name, weight, weight_pos, weight_neg) ;
00026     return ;
00027   }
00028 
00029   //================================================================================
00030 
00031   Stat::Stat(const char *name) : Processor(name), _event(0), _parent_stat(0)  { 
00032 
00037 
00038         if (!STAT) {
00039           STAT = this ;
00040         } else {
00041           throw std::runtime_error ("You may have only a single STAT object around at a time! First Stat object was named " + STAT->name() + " and this second one is named " + name + ".");
00042         }
00043   }
00044 
00045   Stat::Stat (const char *name, bool ignore_duplicate)
00046         : Processor (name), _event(0), _parent_stat(0)
00047   {
00048         if (!STAT) {
00049           STAT = this ;
00050         } else {
00051           if (!ignore_duplicate) {
00052                 throw std::runtime_error ("You may have only a single STAT object around at a time! First Stat object was named " + STAT->name() + " and this second one is named " + name + ".");
00053           }
00054         }
00055   }
00056 
00057   Stat::~Stat() {} ;
00058 
00059   void Stat::chain(void)
00060   {
00062         assert (_parent_stat == 0);
00063         assert (this != STAT);
00064 
00065         _parent_stat = STAT;
00066         STAT = this;
00067   }
00068 
00069   void Stat::unchain(void)
00070   {
00072         assert (_parent_stat != 0);
00073         assert (this == STAT);
00074 
00075         STAT = _parent_stat;
00076         _parent_stat = 0;
00077   }
00078 
00082   void Stat::inheritWeights (void)
00083   {
00084         assert (_parent_stat != 0);
00085         assert (this == STAT);
00086 
00090 
00091         const StatSample *stat_sample = _parent_stat->get_samples().front();
00092 
00093         int num_selections = stat_sample->size();
00094         for (int i = 0; i < num_selections; i++) {
00095           const StatSelection *stat = stat_sample->eventSelection(i);
00096           if (stat->isWeight()) {
00097                 const StatWeight *stat_wt = dynamic_cast<const StatWeight*> (stat);
00098                 double wt = stat_wt->weight();
00099                 double wt_neg = stat_wt->weight_neg();
00100                 double wt_pos = stat_wt->weight_pos();
00101 
00102                 applyWeight (stat_wt->name(), wt, wt_neg, wt_pos);
00103           } else {
00104                 EventSelected (stat->name());
00105           }
00106         }
00107   }
00108 
00109   void Stat::begin() 
00110   {
00111       Config config(name());
00112         
00113       // Get .AddSample entry
00114       vector<string> sample = config.getVString("Sample"," ,");
00115       for (vector<string>::const_iterator it = sample.begin() ;
00116            it != sample.end() ; it++) {
00117 
00118           string s = *it ;
00119           while (s.size() > 0 && s.substr(0,1) == " ") s.erase(0,1) ;
00120           while (s.size() > 0 && s.substr(s.size()-1,1) == " ") s.erase(s.size()-1,1) ;      
00121           if (s.size() ==0 ) continue ;
00122 
00123           _samples.push_back(s) ;
00124 
00125           Config sample_config(s);
00126 
00127           vector<string> tags = sample_config.getVString("UseTag"," ,");
00128           if (tags.size()>0) _samples.back().AddTags(tags) ;
00129 
00130           vector<string> tags_and = sample_config.getVString("AndTag"," ,");
00131           if (tags_and.size()>0) _samples.back().AddAndTags(tags_and) ;
00132       }
00133     
00135       vector<string> syst = config.getVString("Systematics"," ,");      
00136       for( vector<string>::const_iterator it = syst.begin() ;
00137            it != syst.end(); it++) add_syst(*it) ;    
00138 
00139       if (_samples.size() == 0)
00140           _samples.push_back(StatSample(name())) ;
00141 
00142       tag(config.getVString("InitialTag",",")) ;
00143     
00144       _ignoreauto = config.get("IgnoreAutoRecords",false) ;
00145       _precision = config.get("Precision",3) ;
00146       _update = config.get("Update",0) ;
00147 
00148       _output_name = config.get("Output", "efficiencies");
00149       _title = config.get("TexTitle","");
00150 
00151       // read this variable first here just to avoid the error message "not used" in cafe
00152       config.get("PrintSystematics", false) ;
00153   }
00154 
00155   void Stat::finish(){
00156 
00157     ofstream html( (_output_name + ".html").c_str() );
00158     html.setf(ios::fixed);
00159     html.precision(_precision) ;
00160     print_html(html) ;
00161 
00162     if (_title != "") {
00163       ofstream tex( (_output_name + ".tex").c_str() );
00164       tex.setf(ios::fixed);
00165       tex.precision(_precision) ;
00166       print_tex(tex, _title) ;
00167       tex.close() ;
00168     }
00169   }
00170 
00171   void Stat::eventEnd() {
00172     _event = 0 ;
00173 
00174     for (vector<StatSample>::iterator it = _samples.begin() ;
00175          it != _samples.end() ; it++) it->Clear() ;
00176 
00177     for (vector<Syst>::iterator it = _syst.begin() ;
00178          it != _syst.end() ; it++) {
00179       it->pos.Clear() ;
00180       it->neg.Clear() ;
00181     }    
00182   }
00183 
00184   bool Stat::processEvent(cafe::Event& event) {
00185 
00186     StatPointer stat ;
00187     if (event.get("StatPointer", stat)) {
00188       err() << "ERROR! Stat: [" << name() << "] An another instance of the Stat proccesor is running! "
00189             << "Only one instance of the Stat processor must be used!" << endl;
00190       throw runtime_error("Stat:: Only one instance allow for the cafe::Stat processor!") ;
00191     }
00192     stat = this ;
00193     event.put("StatPointer", stat) ;
00194 
00195 
00196     //--------------------------------------------------
00197 
00198     _event = &event ;
00199 
00200     // tag event
00201     for (vector<string>::const_iterator it = _tags.begin() ;
00202          it!= _tags.end(); it++) _event->tag(*it)  ;
00203 
00204     for (vector<StatSample>::iterator it = _samples.begin() ;
00205          it != _samples.end() ; it++)  
00206       it->add(_event) ;
00207 
00208     for (vector<Syst>::iterator it = _syst.begin() ;
00209          it != _syst.end() ; it++) {
00210       it->pos.add(_event) ;
00211       it->neg.add(_event) ;
00212     }
00213 
00214     if (_update > 0 && eventCount()%_update == 0) finish() ;
00215 
00216     return true ;
00217   }
00218 
00219   void Stat::EventSelected(const string& selection_name) {
00220     if (!_event)  return ;
00221     for (vector<StatSample>::iterator it = _samples.begin() ;
00222          it != _samples.end() ; it++) {
00223       if (!_ignoreauto || selection_name.substr(0,9) != "PROCESSOR")
00224         it->add(_event, selection_name) ;
00225     }
00226     for (vector<Syst>::iterator it = _syst.begin() ;
00227          it != _syst.end() ; it++) {
00228       if (!_ignoreauto || selection_name.substr(0,9) != "PROCESSOR") {
00229         it->pos.add(_event, selection_name) ;
00230         it->neg.add(_event, selection_name) ;
00231       }
00232     }
00233     return ; 
00234   }
00235 
00236 
00237   void Stat::applyWeight(const std::string& weight_name, 
00238                          double weight, double weight_pos, double weight_neg) {
00239     for (vector<StatSample>::iterator it = _samples.begin() ;
00240          it != _samples.end() ; it++) 
00241       it->applyWeight(_event, weight_name, weight, weight_pos, weight_neg) ;
00242 
00243     for (vector<Syst>::iterator it = _syst.begin() ;
00244          it != _syst.end() ; it++) {
00245       it->pos.applyWeight(_event, weight_name, weight, weight_pos, weight_neg) ;
00246       it->neg.applyWeight(_event, weight_name, weight, weight_pos, weight_neg) ;
00247     }
00248 
00249     return ; 
00250   }
00251 
00252   //================================================================================
00253 
00255   ostream& Stat::print_html (ostream& os) const {
00256 
00257     // head of the page ===========================================================
00258     os << "<!DOCTYPE HTML PUBLIC \"-//W3C//DTD HTML 4.0 Transitional//EN\">" << endl ;
00259     os << "<HTML><HEAD><TITLE>ANALYSIS EFFICIENCIES</TITLE>" << endl ;
00260     os << "<META content=\"text/html; charset=windows-1252\" http-equiv=Content-Type>" << endl ;
00261     os << "</HEAD>" << endl ;
00262     os << "<BODY aLink=#006600 bgColor=#fffaf0" << endl ;
00263     os << "link=\"blue\" text=\"black\" vLink=\"blue\">" << endl ;
00264     os << " <H2 ALIGN=CENTER> " << endl ;
00265     os << "ANALYSIS EFFICIENCIES" << endl ;
00266     os << "</H2>" << endl ;
00267     os << "" << endl ;
00268     
00269     os << "<TABLE align=right border=0 cellPadding=10 width=\"95%\">" << endl ;
00270     os << "  <TBODY>" << endl ;
00271     os << "  <TR>" << endl ;
00272     os << "    <TD WIDTH=\"60%\">" << endl ;
00273     os << "  </TD>" << endl ;
00274     os << "  <TD>last updated:" << endl ;
00275     os << "      <SCRIPT language=JavaScript>" << endl ;
00276     os << " <!-- " << endl ;
00277     os << "document.writeln(document.lastModified);" << endl ;
00278     os << "// -->" << endl ;
00279     os << "</SCRIPT>" << endl ;
00280     os << "       </TD></TR></TBODY></TABLE> <BR>" << endl ;
00281     os << " " << endl ;
00282     
00283     os << "<HR noShade>" << endl ;
00284     
00285     //--------------- Results table ------------------------------
00286 
00287     bool only_one = true ;
00288     
00289     if (_samples.size() > 1) {
00290       for (vector<StatSample>::const_iterator it = _samples.begin()+1 ;
00291            it != _samples.end() ; it++) 
00292         if (!_samples.begin()->compareNames(*it)) {
00293           only_one = false ;
00294           break ;
00295         }
00296     }
00297   
00298     // number of tabels per page
00299     size_t ntables = 3 ;
00300     if(!only_one) ntables = 2 ;       
00301 
00302     for (size_t ns = 0 ; ns < _samples.size(); ns += ntables) {
00303       
00304       vector<StatSample>::const_iterator itbegin = _samples.begin() + ns;
00305       vector<StatSample>::const_iterator itend  ;
00306       if (ns+ntables < _samples.size()) itend = _samples.begin() + ns + ntables ;
00307       else itend = _samples.end() ;
00308       
00309       // table width in pixels
00310       int column_width = 450 ;
00311       int width = column_width ;
00312       if (_samples.size() - ns == 1) width = 2*width  ;
00313       else if (only_one) width = (1+min(_samples.size() - ns, ntables))*width ;
00314       else width = min(_samples.size() - ns, ntables)*2*width ;
00315       
00316       os << "<TABLE BORDERCOLORDARK=\"996633\"  BORDERCOLOR=\"CC9966\" BORDERCOLORLIGHT=\"FFCC99\"   border=3 cellPadding=5 width=" << width << ">" << endl ;
00317       os << "<TBODY align=right>" << endl ;
00318       
00319       os << "<TR>" << endl ;
00320       for (vector<StatSample>::const_iterator jt = itbegin ;
00321            jt != itend ; jt++) {
00322         
00323         bool weightMode = jt->eventWeight()->nevents() > 0 ?  true : false ;
00324 
00325         if (!only_one || jt == itbegin) {
00326           os << "<TH align=left>" << endl ;
00327           os << "SELECTION"   ;
00328           os << " </TH>" << endl ;
00329         }
00330         
00331         if (weightMode) 
00332           os << "<TH align=center colspan=4>" << endl ; 
00333         else 
00334           os << "<TH align=center colspan=3>" << endl ; 
00335 
00336         os <<  jt->name() << endl ;
00337         os << "</TH>" << endl ; 
00338       }
00339 
00340       os << "</TR>" << endl ;
00341       
00342       bool color = false ;
00343       
00344       unsigned int selsize = 1 ;
00345       for (vector<StatSample>::const_iterator jt = itbegin ;
00346              jt != itend ; jt++) 
00347         selsize = max(jt->size(),selsize) ;
00348 
00349       for (unsigned int n = 0; n < selsize ; n++) {
00350         
00351         if (color) {
00352           os << "<TR BGCOLOR=#bfefff>" << endl ; 
00353           color = !color ;
00354         } else { 
00355           os << "<TR>" << endl ; 
00356           color = !color ;
00357         }      
00358      
00359       for (vector<StatSample>::const_iterator jt = itbegin ;
00360              jt != itend ; jt++) {
00361 
00362         bool weightMode = jt->eventWeight()->nevents() > 0 ?  true : false ;
00363         
00364         if (n >= jt->size()) {
00365           if (!only_one || jt == itbegin) os << "<TD>&#160;</TD> " ;        
00366           os << "<TD>&#160;</TD> <TD>&#160;</TD> <TD>&#160;</TD>" ;
00367           if (weightMode) os << "<TD>&#160;</TD> " ;
00368           os << endl ; 
00369           continue ;
00370         }
00371 
00372         if (!only_one || jt == itbegin) {
00373           os << "<TD align=left>" << endl ; 
00374           os <<  jt->eventSelection(n)->name() << endl ;
00375           os << "</TD>" << endl ; 
00376         }
00377 
00378           
00379         os << "<TD TITLE=\"Number of events\">" << endl ; 
00380         os <<  jt->nevents(n) << endl ;
00381         if (n==0) {
00382           os << "<TD>&#160;</TD> <TD>&#160;</TD>" ; 
00383           if (weightMode) os << "<TD>&#160;</TD> " ;
00384           os << endl ; 
00385           continue ;
00386         }
00387         
00388         if (jt->eventSelection(n)->isWeight()) {
00389           if (jt->eventSelection(n)->name() == jt->eventWeight()->name()) 
00390             os << "</TD> <TD TITLE=\"Global event weight\">" << endl ; 
00391           else 
00392             os << "</TD> <TD TITLE=\"Average weight\">" << endl ; 
00393           os << jt->eff(n) << "&#177;" << jt->effErr(n)<< endl ;
00394           const StatWeight* w = jt->eventWeight(n) ;
00395           if (w->weight_average_pos()>=0 || 
00396               w->weight_average_neg()>=0 )
00397             os << " +" << w->weight_average_pos() - w->weight_average()
00398                << " -" << w->weight_average() - w->weight_average_neg() ;
00399           os << "</TD> <TD> &#160;</TD>" ;
00400           if (weightMode) os << "<TD>&#160;</TD> " ;
00401           os << endl ; 
00402           continue ;
00403         } 
00404         
00405         os << "</TD><TD TITLE=\"Selection efficiency\">"<< endl ;
00406         os << 100.0*jt->eff(n) << "&#177;" << 100.0*jt->effErr(n) << "%"
00407            << endl ;
00408         os << "</TD><TD TITLE=\"Global efficiency\">"<< endl ;
00409         os << 100.0*jt->effGlob(n) << "&#177;" << 100.0*jt->effErrGlob(n) << "%" 
00410            << endl ;
00411         os << "</TD>" << endl ; 
00412         if (weightMode) {
00413           os << "</TD><TD TITLE=\"Global efficiency corrected by the event weight\">"<< endl ;
00414           os << 100.0*jt->correctedEfficiency(n) << "&#177;" << 100.0*jt->correctedEffErr(n) << "%" 
00415              << endl ;
00416           os << "</TD>" << endl ;               
00417         }
00418       }
00419       os << "</TR>" << endl ; 
00420       }
00421 
00422       os << "</TBODY></TABLE>" << endl ;
00423     }
00424 
00425     os << "<BR><BR>" << endl ;
00426 
00427     //--------------- Systematics table ------------------------------
00428 
00429     // table width in pixels
00430     int width = 800 ;
00431     os << " <H3 ALIGN=CENTER> " << endl ;
00432     os << "WEIGHTS AVERAGE VALUES AND SYSTEMATICS AFTER FINAL SELECTION" << endl ;
00433     os << "</H3>" << endl ;
00434 
00435     for (vector<StatSample>::const_iterator jt = _samples.begin() ;
00436          jt != _samples.end() ; jt++) {      
00437 
00438       if(jt->eventWeight()->weight_average()  <= 0 ) continue ;
00439       //      if (jt->eventWeight()->weight_average_pos()  <= 0 &&
00440       //          jt->eventWeight()->weight_average_neg() <=  0 ) continue ;    
00441       
00442       os << "<TABLE BORDERCOLORDARK=\"996633\"  BORDERCOLOR=\"CC9966\" BORDERCOLORLIGHT=\"FFCC99\"   border=3 cellPadding=5 width=" << width << ">" << endl ;
00443       os << "<TBODY align=right>" << endl ;
00444       
00445       os << "<TR>" << endl ;
00446       os << "<TH align=left>" << endl ;
00447       os << "WEIGHT"   ;
00448       os << " </TH>" << endl ;
00449       
00450       os << "<TH align=center colspan=3>" << endl ; 
00451       os <<  jt->name() << endl ;
00452       os << "</TH>" << endl ; 
00453     
00454       os << "</TR>" << endl ;
00455       
00456 
00457       bool color = false ;
00458       for (unsigned int n = 0; n < jt->size() ; n++) {
00459         const StatWeight* w = jt->eventWeight(n) ;
00460         if (!w->isWeight()) continue ;
00461         if (w->weight_average() <0 ) continue ;
00462         //      if (w->weight_average_pos() <0 &&
00463         //          w->weight_average_neg() <0 ) continue ;
00464         if (w->name()=="Initial") continue ;
00465         
00466         if (color) {
00467           os << "<TR BGCOLOR=#bfefff>" << endl ; 
00468           color = !color ;
00469         } else { 
00470           os << "<TR>" << endl ; 
00471           color = !color ;
00472         }      
00473         
00474         os << "<TD align=left>" << endl ; 
00475         os <<  w->name() << endl ;
00476         os << "</TD>" << endl ; 
00477 
00478         double pos = w->weight_average_pos() / w->weight_average() - 1.0 ;
00479         double neg = w->weight_average_neg() / w->weight_average() - 1.0 ;
00480 
00481         os << "<TD TITLE=\"Average central value\">" << endl ; 
00482         os <<   w->weight_average() << endl ;
00483         os << "</TD>" << endl ; 
00484 
00485         os << "<TD TITLE=\"Relative positive shift\">" << endl ; 
00486         if (w->weight_average_pos() < 0)  os <<  " " << endl ;
00487         else os <<  100.0*pos << " % " << endl ;
00488         os << " </TD>" << endl ; 
00489         
00490         os << "<TD TITLE=\"Relative negative shift\">" << endl ; 
00491         if (w->weight_average_neg() < 0)  os <<  " " << endl ;
00492         else os <<  100.0*neg  << " % " << endl ;
00493         os << " </TD>" << endl ; 
00494         
00495         os << "</TR>" << endl ; 
00496       }
00497       os << "</TBODY></TABLE>" << endl ;
00498     }
00499 
00500     // systematics samples   
00501 
00502     if (_syst.size()>0) {
00503     os << "<TABLE BORDERCOLORDARK=\"996633\"  BORDERCOLOR=\"CC9966\" BORDERCOLORLIGHT=\"FFCC99\"   border=3 cellPadding=5 width=1200>" << endl ;
00504     os << "<TBODY align=right>" << endl ;
00505     
00506     os << "<TR>" << endl ;
00507     os << "<TH align=left>SYSTEMATICS</TH>" << endl ;
00508     os << "<TH align=left>SAMPLE</TH>" << endl ;
00509     os << "<TH>POSITIVE</TH>" << endl ;
00510     os << "<TH>NEGATIVE</TH>" << endl ;
00511     os << "</TR>" << endl ;
00512     
00513     for (vector<Syst>::const_iterator jt = _syst.begin() ;
00514          jt != _syst.end() ; jt++) {
00515       os << "<TR>" << endl ;
00516       os << "<TD align=left>" ;
00517       os << jt->name ;
00518       os << "</TD>" << endl ;      
00519       os << "<TD align=left>" ;
00520       os << jt->reference->name() ;
00521       os << "</TD>" << endl ;      
00522 
00523       os << "<TD align>" ;
00524       os << 100*syst_pos(jt) << "&#177;" << 100*systerr_pos(jt) << " %";
00525       os << "</TD>" << endl ;      
00526 
00527       os << "<TD align>" ;
00528       os << 100*syst_neg(jt) << "&#177;" << 100*systerr_neg(jt) << " %";
00529       os << "</TD>" << endl ;      
00530       os << "</TR>" << endl ;
00531     }         
00532 
00533     os << "</TBODY></TABLE>" << endl ;
00534     }
00535     
00536     // print systematics sample selections 
00537       Config config(name());
00538       if (config.get("PrintSystematics", false) ) {
00539         for (vector<Syst>::const_iterator jt = _syst.begin() ;
00540              jt != _syst.end() ; jt++) {
00541           jt->pos.HtmlTable(os) ;
00542           jt->neg.HtmlTable(os) ;
00543         }
00544       }
00545 
00546     os << "</BODY>" << endl ;
00547 
00548     return os ;
00549     
00550   }
00551 
00552   //================================================================================
00553 
00555   ostream& Stat::print_tex (ostream& os, const string title) const {
00556 
00557     for (vector<StatSample>::const_iterator jt = _samples.begin() ;
00558          jt != _samples.end() ; jt++) jt->print_tex(os, title) ;
00559 
00560 
00561     //------ Systematics tables ----------------------------------------
00562 
00563     for (vector<Syst>::const_iterator jt = _syst.begin() ;
00564          jt != _syst.end() ; jt++) {
00565       
00566       os << "" << endl;
00567       os << "% Systematics table %" << endl;
00568       os << "" << endl;
00569       os << "\\begin{table}[p]" << endl;
00570       os << "\\begin{center}" << endl;
00571       os << "\\begin{tabular}{llrr}" << endl;
00572       os << "  \\hline  \\hline" << endl;
00573       os << "  Systematics & Sample & Positive [\\%] & Negative  [\\%] \\\\ \\hline" << endl;      
00574       os << StatSample::tex(jt->name) << " &  " ;
00575       os << StatSample::tex(jt->reference->name()) << " &  " ;
00576       os << " $ " << 100.0*syst_pos(jt)<< "\\pm " << 100.0*systerr_pos(jt) << " $ & $ " 
00577          << 100.0*syst_neg(jt)<< "\\pm " << 100.0*systerr_neg(jt) << " $ \\\\ " << endl ;
00578       os << "  \\hline " << endl;
00579       
00580       os << " \\hline" << endl;
00581       os << " \\end{tabular}" << endl;
00582       os << " \\caption{ Systematics }" << endl;
00583       os << " \\label{Table:syst:}" << endl;
00584       os << "\\end{center}" << endl;
00585       os << "\\end{table}" << endl; 
00586       os << endl;      
00587     }
00588 
00589     return os ;
00590       
00591   }
00592 
00593   double Stat::efficiencyCorrected(const std::string& sample_name) const {
00594     if (sample_name == "not_specified") return _samples.front().correctedEfficiency() ;
00595     for (vector<StatSample>::const_iterator it = _samples.begin() ;
00596          it != _samples.end() ; it++)
00597       if (it->name() == sample_name) return it->correctedEfficiency() ;
00598     cerr << "Stat::efficiencyCorrected ERROR: no sample with name \"" << sample_name
00599           << "\" has been found. " << endl ;
00600     return 0 ;
00601   }
00602                                                         
00603   double Stat::efficiency(const std::string& sample_name) const {
00604     if (sample_name == "not_specified") return _samples.front().effGlob() ;
00605     for (vector<StatSample>::const_iterator it = _samples.begin() ;
00606          it != _samples.end() ; it++)
00607       if (it->name() == sample_name) return it->effGlob() ;
00608     cerr << "Stat::efficiency ERROR: no sample with name \"" << sample_name
00609           << "\" has been found. " << endl ;
00610     return 0 ;
00611   } 
00612 
00613   double Stat::eventWeight(const std::string& name, const std::string sample_name) const {
00614     if (sample_name == "not_specified") 
00615       return _samples.front().eventWeight(name)->weight() ;
00616 
00617     for (vector<StatSample>::const_iterator it = _samples.begin() ;
00618          it != _samples.end() ; it++) 
00619       if (it->name() == sample_name )
00620         return it -> eventWeight(name)->weight() ;
00621 
00622     return -1.0 ;
00623   }
00624   
00625   Collection<EventWeight> Stat::ListEventWeights(const std::string sample_name) const {
00626     if (sample_name == "not_specified") 
00627       return _samples.front().ListEventWeights() ;
00628 
00629     for (vector<StatSample>::const_iterator it = _samples.begin() ;
00630          it != _samples.end() ; it++) 
00631       if (it->name() == sample_name ) 
00632         return it -> ListEventWeights();
00633 
00634     Collection<EventWeight> useless_collection;
00635     return useless_collection;
00636   }
00637   
00638   double Stat::nevents(const string& name, const string& sample_name) const {
00639     if (sample_name == "not_specified") return _samples.front().nevents(name) ;
00640     for (vector<StatSample>::const_iterator it = _samples.begin() ;
00641          it != _samples.end() ; it++) 
00642       if (it->name() == sample_name) return it->nevents(name) ;
00643     cerr << "Stat::nevents ERROR: no sample with name \"" << sample_name
00644           << "\" has been found. " << endl ;
00645     return 0 ;
00646   }
00647 
00648   void Stat::tag(const vector<string>& tags) {
00649     for(vector<string>::const_iterator it = tags.begin();
00650         it != tags.end(); it++) _tags.push_back(*it) ;
00651   }
00652   
00653 
00654   void Stat::add_syst(const string& name, StatSample* sample) {
00655     if(sample) {
00656       _syst.push_back(Syst(name, sample)) ;
00657       sample->AddAndTags(name) ;
00658       return ;
00659     }
00660     
00661     for (vector<StatSample>::iterator it = _samples.begin();
00662          it != _samples.end(); it++) {
00663       _syst.push_back(Syst(name,&(*it))) ;
00664       it->AddAndTags(name) ;
00665     }
00666   }
00667 
00668 
00669   double Stat::syst_pos(std::vector<Syst>::const_iterator jt) const {
00670     if (jt->pos.nevents() <= 0  || jt->reference->nevents() <= 0 || jt->reference->effGlob() <= 0 ) return 0 ;
00671     return jt->pos.effGlob() / jt->reference->effGlob() - 1.0  ;
00672   }
00673 
00674   double Stat::systerr_pos(std::vector<Syst>::const_iterator jt) const {
00675     double s1 = jt->pos.nevents(jt->pos.size()-1) ;
00676     double s2 = jt->reference->nevents(jt->reference->size()-1) ;
00677     if (s1 == s2)  return 0 ;
00678     if (s1 > s2) return sqrt(s2/s1*(s1-s2))/s1 ;
00679     return sqrt(s1/s2*(s2-s1))/s2 ;
00680   }
00681 
00682   double Stat::syst_pos(const std::string& name, const std::string& sample_name) const {
00683     vector<Syst>::const_iterator jt = _syst.begin() ;
00684     for (; jt != _syst.end(); jt++) 
00685       if (jt->name == name && (sample_name == "not_specified" || jt->reference->name() == sample_name))
00686         break ;
00687     if (jt != _syst.end()) return syst_pos(jt) ;
00688 
00689     // if systematics sample with name "name"  was not found, 
00690     // try find weight with name "name"
00691 
00692     vector<StatSample> ::const_iterator it = _samples.begin();    
00693     if (sample_name != "not_specified") {
00694       for ( ; it != _samples.end() ; it++) 
00695         if (it->name() == sample_name) break ;
00696       if (it == _samples.end()) return 0 ;
00697     }
00698     
00699     const StatWeight* w = it->eventWeight(name) ;
00700     if (!w->isWeight()) return 0 ;
00701     
00702     if (w->weight_average_pos() < 0 ) return 0 ;
00703     float err = w->weight_average_pos() - w->weight_average() ;
00704 
00705     // add also  statistical fluctuation  on weight
00706     float sign = err < 0 ? -1.0 : 1.0 ;
00707     return sign*sqrt(err*err +pow(w->err(),2)) ;
00708   }
00709 
00710   double Stat::systerr_pos(const std::string& name, const std::string& sample_name) const {
00711     vector<Syst>::const_iterator jt = _syst.begin() ;
00712     for (; jt != _syst.end(); jt++) 
00713       if (jt->name == name && (sample_name == "not_specified" || jt->reference->name() == sample_name))
00714         break ;
00715     if (jt != _syst.end()) return systerr_pos(jt) ;
00716     return -1.0 ;
00717   }
00718 
00719 
00720   double Stat::syst_neg(std::vector<Syst>::const_iterator jt) const {
00721     if (jt->neg.nevents() <= 0  || jt->reference->nevents() <= 0 || jt->reference->effGlob() <= 0 ) return 0 ;
00722     return jt->neg.effGlob() / jt->reference->effGlob() - 1.0  ;
00723   }
00724 
00725   double Stat::systerr_neg(std::vector<Syst>::const_iterator jt) const {
00726     double s1 = jt->neg.nevents(jt->neg.size()-1) ;
00727     double s2 = jt->reference->nevents(jt->reference->size()-1) ;
00728     if (s1 == s2)  return 0 ;
00729     if (s1 > s2) return sqrt(s2/s1*(s1-s2))/s1 ;
00730     return sqrt(s1/s2*(s2-s1))/s2 ;
00731   }
00732 
00733   double Stat::syst_neg(const std::string& name, const std::string& sample_name) const {
00734     vector<Syst>::const_iterator jt = _syst.begin() ;
00735     for (; jt != _syst.end(); jt++) 
00736       if (jt->name == name && (sample_name == "not_specified" || jt->reference->name() == sample_name))
00737         break ;
00738     if (jt != _syst.end()) return syst_neg(jt) ;
00739 
00740     // if systematics sample with name "name"  was not found, 
00741     // try find weight with name "name"
00742 
00743     vector<StatSample> ::const_iterator it = _samples.begin();    
00744     if (sample_name != "not_specified") {
00745       for ( ; it != _samples.end() ; it++) 
00746         if (it->name() == sample_name) break ;
00747       if (it == _samples.end()) return 0 ;
00748     }
00749     
00750     const StatWeight* w = it->eventWeight(name) ;
00751     if (!w->isWeight()) return 0 ;
00752     
00753     if (w->weight_average_neg() < 0 ) return 0 ;
00754     float err = w->weight_average_neg() - w->weight_average() ;
00755 
00756     // add also  statistical fluctuation  on weight
00757     float sign = err < 0 ? -1.0 : 1.0 ;
00758     return sign*sqrt(err*err +pow(w->err(),2)) ;
00759 
00760   }
00761 
00762   double Stat::systerr_neg(const std::string& name, const std::string& sample_name) const {
00763     vector<Syst>::const_iterator jt = _syst.begin() ;
00764     for (; jt != _syst.end(); jt++) 
00765       if (jt->name == name && (sample_name == "not_specified" || jt->reference->name() == sample_name))
00766         break ;
00767     if (jt != _syst.end()) return systerr_neg(jt) ;
00768     return -1.0 ;
00769   }
00770 
00771   vector<const StatSample*> Stat::get_samples() const {
00772     vector<const StatSample*> s ;
00773     for (vector<StatSample>::const_iterator it = _samples.begin();
00774          it != _samples.end(); it++) s.push_back(&(*it)) ;
00775     return s ;
00776   }
00777 
00778   vector<const Syst*> Stat::get_syst() const {
00779     vector<const Syst*> s ;
00780     for (vector<Syst>::const_iterator it = _syst.begin();
00781          it != _syst.end(); it++) s.push_back(&(*it)) ;
00782     return s ;
00783   }
00784 
00785 }
00786 
00787 ClassImp(cafe::Stat);
00788 ClassImp(cafe::StatPointer);