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