VHbb/interface/controlRegions.h

#ifndef CONTROLREGION_H
#define CONTROLREGION_H

#include "samples.hpp"
#include "ntupleReader.hpp"
#include "CutsAndHistos.h"
#include "TH1.h"
#include <string>

class controlRegion{

 public:
  controlRegion(){ init(); };
 controlRegion( double lumi_, double f_a, double f_b )
   :lumi(lumi_), fa(f_a), fb(f_b){ init(); };
  ~controlRegion() {};

  void setName(std::string name_){ name = name_; }
  std::string getName(){ return name; }

  double sqrtSum( double ex, double ey ){
    if(ex != 0 || ey != 0)
      return sqrt(ex*ex + ey*ey) ;
    else 
      return 0;
  }

  double cData(){return data;}  
  double cSignal(){return signal;}
  double cDYL(){return dyl;}
  double cDYC(){return dyc;}
  double cDYB(){return dyb;}
  double cTTbar(){return ttbar;}
  double cVV(){return vv;}
  double cST(){return st;}
  double cWJETS(){return wjets;}
  double cOthers(){return others;}
  double cRest(){ return ( data-(vv+st+wjets+others) );}
  double cTotal(){return dyl+dyc+dyb+ttbar+vv+st+wjets+others;}

  TH1F * hData(){ hdata = new TH1F(*(TH1F*)(list_data->At(0))); hdata->Reset(); hdata->Merge(list_data); return hdata;}  
  TH1F * hSignal(){ hsignal = new TH1F(*(TH1F*)(list_signal->At(0))); hsignal->Reset(); hsignal->Merge(list_signal); return hsignal;}
  TH1F * hDYL(){ hdyl = new TH1F(*(TH1F*)(list_dyl->At(0))); hdyl->Reset(); hdyl->Merge(list_dyl); return hdyl;}
  TH1F * hDYC(){ hdyc = new TH1F(*(TH1F*)(list_dyc->At(0))); hdyc->Reset(); hdyc->Merge(list_dyc); return hdyc;}
  TH1F * hDYB(){ hdyb = new TH1F(*(TH1F*)(list_dyb->At(0))); hdyb->Reset(); hdyb->Merge(list_dyb); return hdyb;}
  TH1F * hTTbar(){ httbar = new TH1F(*(TH1F*)(list_ttbar->At(0))); httbar->Reset(); httbar->Merge(list_ttbar); return httbar;}
  TH1F * hVV(){ hvv = new TH1F(*(TH1F*)(list_vv->At(0))); hvv->Reset(); hvv->Merge(list_vv); return hvv;}
  TH1F * hST(){ hst = new TH1F(*(TH1F*)(list_st->At(0))); hst->Reset(); hst->Merge(list_st);  return hst;}
  TH1F * hWJETS(){hwjets = new TH1F(*(TH1F*)(list_wjets->At(0))); hwjets->Reset(); hwjets->Merge(list_wjets); return hwjets;}
  TH1F * hOthers(){ hothers = new TH1F(*(TH1F*)(list_others->At(0))); hothers->Reset(); hothers->Merge(list_others); return hothers;}
  TH1F * hTotal(){ htotal = new TH1F(*(TH1F*)(list_total->At(0))); htotal->Reset(); htotal->Merge(list_total); return htotal;}

  double eData(){if(data > 0) err_data = sqrt(data); return err_data;}  
  double eSignal(){return err_signal;}
  double eDYL(){return err_dyl;}
  double eDYC(){return err_dyc;}
  double eDYB(){return err_dyb;}
  double eTTbar(){return err_ttbar;}
  double eVV(){return err_vv = sqrtSum( err_vv, 0.3*vv );}
  double eST(){return err_st = sqrtSum( err_st, 0.3*st );}
  double eWJETS(){return err_wjets;}
  double eOthers(){return err_others;}
  double eRest(){ return sqrt( eData()*eData() + eVV()*eVV() + eST()*eST() + eWJETS()*eWJETS() + eOthers()*eOthers() );}
  double eTotal(){ return sqrt( eDYL()*eDYL() + eDYC()*eDYC() + eDYB()*eDYB() + eTTbar()*eTTbar() + eVV()*eVV() + eST()*eST() + eWJETS()*eWJETS() + eOthers()*eOthers() );}

  double count( std::string type ){
    if( ((int)type.find("+")) > 0 ) return combinedCount( type );
    if( type == Run2011A_s || type == Run2011B_s ){ return cData(); }
    else{
      if(type == Signal_s){ return cSignal(); }
      else if(type == DYL_s ){ return cDYL(); }
      else if(type == DYC_s ){ return cDYC(); }
      else if(type == DYB_s ){  return cDYB(); }
      else if(type == TTbar_s){  return cTTbar(); }
      else if(type == VV_s){ return cVV(); }
      else if(type == ST_s){ return cST();  }
      else if(type == WJETS_s){ return cWJETS();  }
      else{  return cOthers();  }
    }
  }
  double combinedCount( std::string type ){
    count_comb = 0;
    std::vector<std::string> types;
    size_t pos = type.size();
    while ( (int)pos > 0){
      std::cout << type << " has size = " << (int)pos << std::endl; 
      pos = type.rfind("+");
      std::cout << type << " has size = " << (int)pos << std::endl; 
      if( (int)pos > 0){
        types.push_back( type.substr(pos+1) ); 
        type.resize(pos);
      }
      else
        types.push_back( type );
    }
    for(size_t t=0; t<types.size(); ++t){
      std::cout << types.at(t) << std::endl;
      count_comb += count( types.at(t) );
    }
    return count_comb;
  }

  TH1F * histo( std::string type ){
    if( ((int)type.find("+")) > 0 ) return combinedHistos( type );    
    if( type == Run2011A_s || type == Run2011B_s ){ return hData(); }
    else{
      if(type == Signal_s){ return hSignal(); }
      else if(type == DYL_s ){ return hDYL(); }
      else if(type == DYC_s ){ return hDYC(); }
      else if(type == DYB_s ){  return hDYB(); }
      else if(type == TTbar_s){  return hTTbar(); }
      else if(type == VV_s){ return hVV(); }
      else if(type == ST_s){ return hST();  }
      else if(type == WJETS_s){ return hWJETS();  }
      else{  return hOthers();  }
    }
  }

  TH1F * combinedHistos( std::string type ){
    std::cout << "Combine histos " << std::endl;
    std::vector<std::string> types;
    size_t pos = type.size();
    while ((int)pos > 0){
      std::cout << type << " has size = " << (int)pos << std::endl; 
      pos = type.rfind("+");
      std::cout << type << " has size = " << (int)pos << std::endl; 
      if( (int)pos > 0 ){
        types.push_back( type.substr(pos+1) ); 
        type.resize(pos);
      }
      else
        types.push_back( type );
    }
    for(size_t t=0; t<types.size(); ++t)
      list_comb->Add( histo(types.at(t)) );
    h_comb = new TH1F(*(TH1F*)(list_comb->At(0))); h_comb->Reset(); h_comb->Merge(list_comb);  return h_comb;
  }

  void dump(){
    std::cout << "Data            = " << cData() << std::endl;
    std::cout << "Signal          = " << cSignal() << std::endl;
    std::cout << "DrellYann light = " << cDYL() << std::endl;
    std::cout << "DrellYann charm = " << cDYC() << std::endl;
    std::cout << "DrellYann b     = " << cDYB() << std::endl;
    std::cout << "TTbar           = " << cTTbar() << std::endl;
    std::cout << "Diboson         = " << cVV() << std::endl;
    std::cout << "Single Top      = " << cST() << std::endl;
    std::cout << "Wjets           = " << cWJETS() << std::endl;
    std::cout << "Others          = " << cOthers() << std::endl;
    std::cout << "Total MC        = " << cTotal() << std::endl;
  }


  void fillFromHisto( Sample & sample, TH1F & histo, double min=-1e10, double max=1e10 ){
    if(sample.data){
      data += histo.Integral(min,max);
      list_data->Add(&histo);
    }
    else{
      if(sample.name == Signal_s){
        signal = histo.Integral(min,max);
        list_signal->Add(&histo);
        bcsignal=0;
        err_signal = sqrtSum( err_signal, TMath::Sqrt(signal) );
      }
      else{
        list_total->Add(&histo);
        if(sample.name == DYL_s ){
          dyl += histo.Integral(min,max);
          list_dyl->Add(&histo);
          bcdyl=0;
          err_dyl = sqrtSum( err_dyl, TMath::Sqrt(dyl) );
        }
        else if(sample.name == DYC_s ){
          dyc += histo.Integral(min,max);
          list_dyc->Add(&histo);
          bcdyc=0;
          err_dyc = sqrtSum( err_dyl, TMath::Sqrt(dyc) );
        }
        else if(sample.name == DYB_s ){
          dyb += histo.Integral(min,max);
          list_dyb->Add(&histo);
          bcdyb=0;
          err_dyb = sqrtSum( err_dyl, TMath::Sqrt(dyb) );
        }
        else if(sample.name == TTbar_s){
          ttbar += histo.Integral(min,max);
          list_ttbar->Add(&histo);
          bcttbar=0;
          err_ttbar = sqrtSum( err_ttbar, TMath::Sqrt(ttbar) );
        }
        else if(sample.name == VV_s){
          vv +=  histo.Integral(min,max);
          list_vv->Add(&histo);
          bcvv=0;
          err_vv = sqrtSum( err_vv, TMath::Sqrt(vv) );
        }
        else if(sample.name == ST_s){
          st +=  histo.Integral(min,max);
          list_st->Add(&histo);
          bcst=0;
          err_st = sqrtSum( err_st, TMath::Sqrt(st) );
        }
        else if(sample.name == WJETS_s){
          wjets +=  histo.Integral(min,max);
          list_wjets->Add(&histo);
          bcwjets=0;
          err_wjets = sqrtSum( err_wjets, TMath::Sqrt(wjets) );
        }
        else{
          std::cout << "WARNING: Others counter non zero!" << std::endl;
          others +=  histo.Integral(min,max);
          list_others->Add(&histo);
          bcothers=0;
          err_others = sqrtSum( err_others, TMath::Sqrt(others) );
        }
      }
    }
  }


  void fill( Sample &sample, CutSample &cut, PCutSet &addCut, int n, ntupleReader &ev){
    bool ok = cut.pass( ev, sample );
    bool addOk = addCut.pass( ev, n );
    double scale = sample.scale(lumi,fa,fb);
    double weight = sample.scale(lumi,fa,fb) * cut.weight( ev , sample );
    if(!ok || !addOk)
      return;
    else{
      if(sample.data)
        data++;
      else{
        if(sample.name == Signal_s){
          signal += weight;
          bcsignal++;
          err_signal = sqrtSum( err_signal, weight );
        }
        else if(sample.name == DY_s || sample.name == DYBOOSTED_s){
          if( TMath::Abs(ev.eventFlav) < 4){
            dyl += weight;
            bcdyl++;
            err_dyl = sqrtSum( err_dyl, weight );
          }
          else if( TMath::Abs(ev.eventFlav) == 4){
            dyc += weight;
            bcdyc++;
            err_dyc = sqrtSum( err_dyc, weight );
          }
          else if( TMath::Abs(ev.eventFlav) == 5){
            dyb += weight;
            bcdyb++;
            err_dyb = sqrtSum( err_dyb, weight );
          }
        }
        else if(sample.name == TTbar_s){
          ttbar += weight;
          bcttbar++;
          err_ttbar = sqrtSum( err_ttbar, weight );
        }
        else if(sample.name == VV_s){
          vv += weight;
          bcvv++;
          err_vv = sqrtSum( err_vv, weight );
        }
        else if(sample.name == ST_s){
          st += weight;
          bcst++;
          err_st = sqrtSum( err_st, weight );
        }
        else if(sample.name == WJETS_s){
          wjets += weight;
          bcwjets++;
          err_wjets = sqrtSum( err_wjets, weight );
        }
        else{
          std::cout << "WARNING: Others counter non zero!" << std::endl;
          others += weight;
          bcothers++;
          err_others = sqrtSum( err_others, weight );
        }
      }
    }
  }

  void fill( Sample &sample, CutSample &cut,  ntupleReader &ev){
    bool ok = cut.pass( ev, sample );
    double scale = sample.scale(lumi,fa,fb);
    double weight = sample.scale(lumi,fa,fb) * cut.weight( ev , sample );
    if(!ok)
      return;
    else{
      if(sample.data)
        data++;
      else{
        if(sample.name == Signal_s){
          signal += weight;
          bcsignal++;
          err_signal = sqrtSum( err_signal, weight );
        }
        else if(sample.name == DY_s || sample.name == DYBOOSTED_s){
          if( TMath::Abs(ev.eventFlav) < 4){
            dyl += weight;
            bcdyl++;
            err_dyl = sqrtSum( err_dyl, weight );
          }
          else if( TMath::Abs(ev.eventFlav) == 4){
            dyc += weight;
            bcdyc++;
            err_dyc = sqrtSum( err_dyc, weight );
          }
          else if( TMath::Abs(ev.eventFlav) == 5){
            dyb += weight;
            bcdyb++;
            err_dyb = sqrtSum( err_dyb, weight );
          }
        }
        else if(sample.name == TTbar_s){
          ttbar += weight;
          bcttbar++;
          err_ttbar = sqrtSum( err_ttbar, weight );
        }
        else if(sample.name == VV_s){
          vv += weight;
          bcvv++;
          err_vv = sqrtSum( err_vv, weight );
        }
        else if(sample.name == ST_s){
          st += weight;
          bcst++;
          err_st = sqrtSum( err_st, weight );
        }
        else if(sample.name == WJETS_s){
          wjets += weight;
          bcwjets++;
          err_wjets = sqrtSum( err_wjets, weight );
        }
        else{
          std::cout << "WARNING: Others counter non zero!" << std::endl;
          others += weight;
          bcothers++;
          err_others = sqrtSum( err_others, weight );
        }
      }
    }
  }


  void fill( Sample& sample, double entries){
    double weight = entries;
    if(sample.data)
      data++;
    else{
      if(sample.name == Signal_s){
        signal += weight;
        bcsignal++;
        err_signal = sqrtSum( err_signal, weight );
      }
      else if(sample.name == DYL_s ){
        dyl += weight;
        bcdyl++;
        err_dyl = sqrtSum( err_dyl, weight );
      }
      else if(sample.name == DYC_s ){
        dyc += weight;
        bcdyc++;
        err_dyc = sqrtSum( err_dyc, weight );
      }
      else if(sample.name == DYB_s ){
        dyb += weight;
        bcdyb++;
        err_dyb = sqrtSum( err_dyb, weight );
      }
      else if(sample.name == TTbar_s){
        ttbar += weight;
        bcttbar++;
        err_ttbar = sqrtSum( err_ttbar, weight );
      }
      else if(sample.name == VV_s){
        vv += weight;
        bcvv++;
        err_vv = sqrtSum( err_vv, weight );
      }
      else if(sample.name == ST_s){
        st += weight;
        bcst++;
        err_st = sqrtSum( err_st, weight );
      }
      else if(sample.name == WJETS_s){
        wjets += weight;
        bcwjets++;
        err_wjets = sqrtSum( err_wjets, weight );
      }
      else{
        std::cout << "WARNING: Others counter non zero!" << std::endl;
        others += weight;
        bcothers++;
        err_others = sqrtSum( err_others, weight );
      }
    }
  }

  
  void init(){

    list_data = new TList;
    list_signal = new TList;
    list_dyl = new TList;
    list_dyc = new TList;
    list_dyb = new TList;
    list_ttbar = new TList;
    list_vv = new TList;
    list_st = new TList;
    list_wjets = new TList;
    list_others = new TList;
    list_total = new TList;
    list_comb = new TList;

   data=0;  
   signal=0;
   dyl=0;
   dyc=0;
   dyb=0;
   ttbar=0;
   vv=0;
   st=0;
   wjets=0;
   others=0;

   bcdata=0;  
   bcsignal=0;
   bcdyl=0;
   bcdyc=0;
   bcdyb=0;
   bcttbar=0;
   bcvv=0;
   bcst=0;
   bcwjets=0;
   bcothers=0;

   err_data=0;  
   err_signal=0;
   err_dyl=0;
   err_dyc=0;
   err_dyb=0;
   err_ttbar=0;
   err_vv=0;
   err_st=0;
   err_wjets=0;
   err_others=0;
   
   Run2011A_s = "Run2011A";
   Run2011B_s = "Run2011B";
   Signal_s = "ZH";
   DYL_s = "DYL";
   DYC_s = "DYC";
   DYB_s = "DYB";
   DY_s = "DY";
   DYBOOSTED_s = "DYBOOSTED";
   TTbar_s = "TTbar";
   VV_s = "VV";
   ST_s = "ST";
   WJETS_s = "WJ";

  }
  
 private:

  double lumi;
  double fa;
  double fb;
  std::string name;

  std::string Run2011A_s;
  std::string Run2011B_s;
  std::string Signal_s;
  std::string DYL_s;
  std::string DYC_s;
  std::string DYB_s;
  std::string DY_s;
  std::string DYBOOSTED_s;
  std::string TTbar_s;
  std::string VV_s;
  std::string ST_s;
  std::string WJETS_s;

  double data;  
  double signal;
  double dyl;
  double dyc;
  double dyb;
  double ttbar;
  double vv;
  double st;
  double wjets;
  double others;
  double count_comb;

  double bcdata;  
  double bcsignal;
  double bcdyl;
  double bcdyc;
  double bcdyb;
  double bcttbar;
  double bcvv;
  double bcst;
  double bcwjets;
  double bcothers;

  TH1F * hdata;  
  TH1F * hsignal;
  TH1F * hdyl;
  TH1F * hdyc;
  TH1F * hdyb;
  TH1F * httbar;
  TH1F * hvv;
  TH1F * hst;
  TH1F * hwjets;
  TH1F * hothers;
  TH1F * htotal;
  TH1F * h_comb;

  TList * list_data;  
  TList * list_signal;
  TList * list_dyl;
  TList * list_dyc;
  TList * list_dyb;
  TList * list_ttbar;
  TList * list_vv;
  TList * list_st;
  TList * list_wjets;
  TList * list_others;
  TList * list_total;
  TList * list_comb;

  double err_data;  
  double err_signal;
  double err_dyl;
  double err_dyc;
  double err_dyb;
  double err_ttbar;
  double err_vv;
  double err_st;
  double err_wjets;
  double err_others;


};


#endif
Revision:	1.2
Committed:	Wed May 2 19:20:07 2012 UTC (13 years ago) by bortigno
Content type:	text/plain
Branch:	MAIN
CVS Tags:	lhcp_UnblindFix, hcp_Unblind, lhcp_11April, LHCP_PreAppFixAfterFreeze, LHCP_PreAppFreeze, workingVersionAfterHCP, hcpApproval, hcpPreApp, ICHEP8TeV, ichep8TeV, AN-12-181-7TeV_patch1, AN-12-181-7TeV, HEAD
Changes since 1.1:	+3 -3 lines
Log Message:	working version
#	User	Rev	Content
1	bortigno	1.1	#ifndef CONTROLREGION_H
2			#define CONTROLREGION_H
3
4	bortigno	1.2	#include "samples.hpp"
5			#include "ntupleReader.hpp"
6			#include "CutsAndHistos.h"
7	bortigno	1.1	#include "TH1.h"
8			#include <string>
9
10			class controlRegion{
11
12			public:
13			controlRegion(){ init(); };
14			controlRegion( double lumi_, double f_a, double f_b )
15			:lumi(lumi_), fa(f_a), fb(f_b){ init(); };
16			~controlRegion() {};
17
18			void setName(std::string name_){ name = name_; }
19			std::string getName(){ return name; }
20
21			double sqrtSum( double ex, double ey ){
22			if(ex != 0 \|\| ey != 0)
23			return sqrt(exex + eyey) ;
24			else
25			return 0;
26			}
27
28			double cData(){return data;}
29			double cSignal(){return signal;}
30			double cDYL(){return dyl;}
31			double cDYC(){return dyc;}
32			double cDYB(){return dyb;}
33			double cTTbar(){return ttbar;}
34			double cVV(){return vv;}
35			double cST(){return st;}
36			double cWJETS(){return wjets;}
37			double cOthers(){return others;}
38			double cRest(){ return ( data-(vv+st+wjets+others) );}
39			double cTotal(){return dyl+dyc+dyb+ttbar+vv+st+wjets+others;}
40
41			TH1F * hData(){ hdata = new TH1F((TH1F)(list_data->At(0))); hdata->Reset(); hdata->Merge(list_data); return hdata;}
42			TH1F * hSignal(){ hsignal = new TH1F((TH1F)(list_signal->At(0))); hsignal->Reset(); hsignal->Merge(list_signal); return hsignal;}
43			TH1F * hDYL(){ hdyl = new TH1F((TH1F)(list_dyl->At(0))); hdyl->Reset(); hdyl->Merge(list_dyl); return hdyl;}
44			TH1F * hDYC(){ hdyc = new TH1F((TH1F)(list_dyc->At(0))); hdyc->Reset(); hdyc->Merge(list_dyc); return hdyc;}
45			TH1F * hDYB(){ hdyb = new TH1F((TH1F)(list_dyb->At(0))); hdyb->Reset(); hdyb->Merge(list_dyb); return hdyb;}
46			TH1F * hTTbar(){ httbar = new TH1F((TH1F)(list_ttbar->At(0))); httbar->Reset(); httbar->Merge(list_ttbar); return httbar;}
47			TH1F * hVV(){ hvv = new TH1F((TH1F)(list_vv->At(0))); hvv->Reset(); hvv->Merge(list_vv); return hvv;}
48			TH1F * hST(){ hst = new TH1F((TH1F)(list_st->At(0))); hst->Reset(); hst->Merge(list_st); return hst;}
49			TH1F * hWJETS(){hwjets = new TH1F((TH1F)(list_wjets->At(0))); hwjets->Reset(); hwjets->Merge(list_wjets); return hwjets;}
50			TH1F * hOthers(){ hothers = new TH1F((TH1F)(list_others->At(0))); hothers->Reset(); hothers->Merge(list_others); return hothers;}
51			TH1F * hTotal(){ htotal = new TH1F((TH1F)(list_total->At(0))); htotal->Reset(); htotal->Merge(list_total); return htotal;}
52
53			double eData(){if(data > 0) err_data = sqrt(data); return err_data;}
54			double eSignal(){return err_signal;}
55			double eDYL(){return err_dyl;}
56			double eDYC(){return err_dyc;}
57			double eDYB(){return err_dyb;}
58			double eTTbar(){return err_ttbar;}
59			double eVV(){return err_vv = sqrtSum( err_vv, 0.3*vv );}
60			double eST(){return err_st = sqrtSum( err_st, 0.3*st );}
61			double eWJETS(){return err_wjets;}
62			double eOthers(){return err_others;}
63			double eRest(){ return sqrt( eData()eData() + eVV()eVV() + eST()eST() + eWJETS()eWJETS() + eOthers()*eOthers() );}
64			double eTotal(){ return sqrt( eDYL()eDYL() + eDYC()eDYC() + eDYB()eDYB() + eTTbar()eTTbar() + eVV()eVV() + eST()eST() + eWJETS()eWJETS() + eOthers()eOthers() );}
65
66			double count( std::string type ){
67			if( ((int)type.find("+")) > 0 ) return combinedCount( type );
68			if( type == Run2011A_s \|\| type == Run2011B_s ){ return cData(); }
69			else{
70			if(type == Signal_s){ return cSignal(); }
71			else if(type == DYL_s ){ return cDYL(); }
72			else if(type == DYC_s ){ return cDYC(); }
73			else if(type == DYB_s ){ return cDYB(); }
74			else if(type == TTbar_s){ return cTTbar(); }
75			else if(type == VV_s){ return cVV(); }
76			else if(type == ST_s){ return cST(); }
77			else if(type == WJETS_s){ return cWJETS(); }
78			else{ return cOthers(); }
79			}
80			}
81			double combinedCount( std::string type ){
82			count_comb = 0;
83			std::vector<std::string> types;
84			size_t pos = type.size();
85			while ( (int)pos > 0){
86			std::cout << type << " has size = " << (int)pos << std::endl;
87			pos = type.rfind("+");
88			std::cout << type << " has size = " << (int)pos << std::endl;
89			if( (int)pos > 0){
90			types.push_back( type.substr(pos+1) );
91			type.resize(pos);
92			}
93			else
94			types.push_back( type );
95			}
96			for(size_t t=0; t<types.size(); ++t){
97			std::cout << types.at(t) << std::endl;
98			count_comb += count( types.at(t) );
99			}
100			return count_comb;
101			}
102
103			TH1F * histo( std::string type ){
104			if( ((int)type.find("+")) > 0 ) return combinedHistos( type );
105			if( type == Run2011A_s \|\| type == Run2011B_s ){ return hData(); }
106			else{
107			if(type == Signal_s){ return hSignal(); }
108			else if(type == DYL_s ){ return hDYL(); }
109			else if(type == DYC_s ){ return hDYC(); }
110			else if(type == DYB_s ){ return hDYB(); }
111			else if(type == TTbar_s){ return hTTbar(); }
112			else if(type == VV_s){ return hVV(); }
113			else if(type == ST_s){ return hST(); }
114			else if(type == WJETS_s){ return hWJETS(); }
115			else{ return hOthers(); }
116			}
117			}
118
119			TH1F * combinedHistos( std::string type ){
120			std::cout << "Combine histos " << std::endl;
121			std::vector<std::string> types;
122			size_t pos = type.size();
123			while ((int)pos > 0){
124			std::cout << type << " has size = " << (int)pos << std::endl;
125			pos = type.rfind("+");
126			std::cout << type << " has size = " << (int)pos << std::endl;
127			if( (int)pos > 0 ){
128			types.push_back( type.substr(pos+1) );
129			type.resize(pos);
130			}
131			else
132			types.push_back( type );
133			}
134			for(size_t t=0; t<types.size(); ++t)
135			list_comb->Add( histo(types.at(t)) );
136			h_comb = new TH1F((TH1F)(list_comb->At(0))); h_comb->Reset(); h_comb->Merge(list_comb); return h_comb;
137			}
138
139			void dump(){
140			std::cout << "Data = " << cData() << std::endl;
141			std::cout << "Signal = " << cSignal() << std::endl;
142			std::cout << "DrellYann light = " << cDYL() << std::endl;
143			std::cout << "DrellYann charm = " << cDYC() << std::endl;
144			std::cout << "DrellYann b = " << cDYB() << std::endl;
145			std::cout << "TTbar = " << cTTbar() << std::endl;
146			std::cout << "Diboson = " << cVV() << std::endl;
147			std::cout << "Single Top = " << cST() << std::endl;
148			std::cout << "Wjets = " << cWJETS() << std::endl;
149			std::cout << "Others = " << cOthers() << std::endl;
150			std::cout << "Total MC = " << cTotal() << std::endl;
151			}
152
153
154			void fillFromHisto( Sample & sample, TH1F & histo, double min=-1e10, double max=1e10 ){
155			if(sample.data){
156			data += histo.Integral(min,max);
157			list_data->Add(&histo);
158			}
159			else{
160			if(sample.name == Signal_s){
161			signal = histo.Integral(min,max);
162			list_signal->Add(&histo);
163			bcsignal=0;
164			err_signal = sqrtSum( err_signal, TMath::Sqrt(signal) );
165			}
166			else{
167			list_total->Add(&histo);
168			if(sample.name == DYL_s ){
169			dyl += histo.Integral(min,max);
170			list_dyl->Add(&histo);
171			bcdyl=0;
172			err_dyl = sqrtSum( err_dyl, TMath::Sqrt(dyl) );
173			}
174			else if(sample.name == DYC_s ){
175			dyc += histo.Integral(min,max);
176			list_dyc->Add(&histo);
177			bcdyc=0;
178			err_dyc = sqrtSum( err_dyl, TMath::Sqrt(dyc) );
179			}
180			else if(sample.name == DYB_s ){
181			dyb += histo.Integral(min,max);
182			list_dyb->Add(&histo);
183			bcdyb=0;
184			err_dyb = sqrtSum( err_dyl, TMath::Sqrt(dyb) );
185			}
186			else if(sample.name == TTbar_s){
187			ttbar += histo.Integral(min,max);
188			list_ttbar->Add(&histo);
189			bcttbar=0;
190			err_ttbar = sqrtSum( err_ttbar, TMath::Sqrt(ttbar) );
191			}
192			else if(sample.name == VV_s){
193			vv += histo.Integral(min,max);
194			list_vv->Add(&histo);
195			bcvv=0;
196			err_vv = sqrtSum( err_vv, TMath::Sqrt(vv) );
197			}
198			else if(sample.name == ST_s){
199			st += histo.Integral(min,max);
200			list_st->Add(&histo);
201			bcst=0;
202			err_st = sqrtSum( err_st, TMath::Sqrt(st) );
203			}
204			else if(sample.name == WJETS_s){
205			wjets += histo.Integral(min,max);
206			list_wjets->Add(&histo);
207			bcwjets=0;
208			err_wjets = sqrtSum( err_wjets, TMath::Sqrt(wjets) );
209			}
210			else{
211			std::cout << "WARNING: Others counter non zero!" << std::endl;
212			others += histo.Integral(min,max);
213			list_others->Add(&histo);
214			bcothers=0;
215			err_others = sqrtSum( err_others, TMath::Sqrt(others) );
216			}
217			}
218			}
219			}
220
221
222			void fill( Sample &sample, CutSample &cut, PCutSet &addCut, int n, ntupleReader &ev){
223			bool ok = cut.pass( ev, sample );
224			bool addOk = addCut.pass( ev, n );
225			double scale = sample.scale(lumi,fa,fb);
226			double weight = sample.scale(lumi,fa,fb) * cut.weight( ev , sample );
227			if(!ok \|\| !addOk)
228			return;
229			else{
230			if(sample.data)
231			data++;
232			else{
233			if(sample.name == Signal_s){
234			signal += weight;
235			bcsignal++;
236			err_signal = sqrtSum( err_signal, weight );
237			}
238			else if(sample.name == DY_s \|\| sample.name == DYBOOSTED_s){
239			if( TMath::Abs(ev.eventFlav) < 4){
240			dyl += weight;
241			bcdyl++;
242			err_dyl = sqrtSum( err_dyl, weight );
243			}
244			else if( TMath::Abs(ev.eventFlav) == 4){
245			dyc += weight;
246			bcdyc++;
247			err_dyc = sqrtSum( err_dyc, weight );
248			}
249			else if( TMath::Abs(ev.eventFlav) == 5){
250			dyb += weight;
251			bcdyb++;
252			err_dyb = sqrtSum( err_dyb, weight );
253			}
254			}
255			else if(sample.name == TTbar_s){
256			ttbar += weight;
257			bcttbar++;
258			err_ttbar = sqrtSum( err_ttbar, weight );
259			}
260			else if(sample.name == VV_s){
261			vv += weight;
262			bcvv++;
263			err_vv = sqrtSum( err_vv, weight );
264			}
265			else if(sample.name == ST_s){
266			st += weight;
267			bcst++;
268			err_st = sqrtSum( err_st, weight );
269			}
270			else if(sample.name == WJETS_s){
271			wjets += weight;
272			bcwjets++;
273			err_wjets = sqrtSum( err_wjets, weight );
274			}
275			else{
276			std::cout << "WARNING: Others counter non zero!" << std::endl;
277			others += weight;
278			bcothers++;
279			err_others = sqrtSum( err_others, weight );
280			}
281			}
282			}
283			}
284
285			void fill( Sample &sample, CutSample &cut, ntupleReader &ev){
286			bool ok = cut.pass( ev, sample );
287			double scale = sample.scale(lumi,fa,fb);
288			double weight = sample.scale(lumi,fa,fb) * cut.weight( ev , sample );
289			if(!ok)
290			return;
291			else{
292			if(sample.data)
293			data++;
294			else{
295			if(sample.name == Signal_s){
296			signal += weight;
297			bcsignal++;
298			err_signal = sqrtSum( err_signal, weight );
299			}
300			else if(sample.name == DY_s \|\| sample.name == DYBOOSTED_s){
301			if( TMath::Abs(ev.eventFlav) < 4){
302			dyl += weight;
303			bcdyl++;
304			err_dyl = sqrtSum( err_dyl, weight );
305			}
306			else if( TMath::Abs(ev.eventFlav) == 4){
307			dyc += weight;
308			bcdyc++;
309			err_dyc = sqrtSum( err_dyc, weight );
310			}
311			else if( TMath::Abs(ev.eventFlav) == 5){
312			dyb += weight;
313			bcdyb++;
314			err_dyb = sqrtSum( err_dyb, weight );
315			}
316			}
317			else if(sample.name == TTbar_s){
318			ttbar += weight;
319			bcttbar++;
320			err_ttbar = sqrtSum( err_ttbar, weight );
321			}
322			else if(sample.name == VV_s){
323			vv += weight;
324			bcvv++;
325			err_vv = sqrtSum( err_vv, weight );
326			}
327			else if(sample.name == ST_s){
328			st += weight;
329			bcst++;
330			err_st = sqrtSum( err_st, weight );
331			}
332			else if(sample.name == WJETS_s){
333			wjets += weight;
334			bcwjets++;
335			err_wjets = sqrtSum( err_wjets, weight );
336			}
337			else{
338			std::cout << "WARNING: Others counter non zero!" << std::endl;
339			others += weight;
340			bcothers++;
341			err_others = sqrtSum( err_others, weight );
342			}
343			}
344			}
345			}
346
347
348			void fill( Sample& sample, double entries){
349			double weight = entries;
350			if(sample.data)
351			data++;
352			else{
353			if(sample.name == Signal_s){
354			signal += weight;
355			bcsignal++;
356			err_signal = sqrtSum( err_signal, weight );
357			}
358			else if(sample.name == DYL_s ){
359			dyl += weight;
360			bcdyl++;
361			err_dyl = sqrtSum( err_dyl, weight );
362			}
363			else if(sample.name == DYC_s ){
364			dyc += weight;
365			bcdyc++;
366			err_dyc = sqrtSum( err_dyc, weight );
367			}
368			else if(sample.name == DYB_s ){
369			dyb += weight;
370			bcdyb++;
371			err_dyb = sqrtSum( err_dyb, weight );
372			}
373			else if(sample.name == TTbar_s){
374			ttbar += weight;
375			bcttbar++;
376			err_ttbar = sqrtSum( err_ttbar, weight );
377			}
378			else if(sample.name == VV_s){
379			vv += weight;
380			bcvv++;
381			err_vv = sqrtSum( err_vv, weight );
382			}
383			else if(sample.name == ST_s){
384			st += weight;
385			bcst++;
386			err_st = sqrtSum( err_st, weight );
387			}
388			else if(sample.name == WJETS_s){
389			wjets += weight;
390			bcwjets++;
391			err_wjets = sqrtSum( err_wjets, weight );
392			}
393			else{
394			std::cout << "WARNING: Others counter non zero!" << std::endl;
395			others += weight;
396			bcothers++;
397			err_others = sqrtSum( err_others, weight );
398			}
399			}
400			}
401
402
403			void init(){
404
405			list_data = new TList;
406			list_signal = new TList;
407			list_dyl = new TList;
408			list_dyc = new TList;
409			list_dyb = new TList;
410			list_ttbar = new TList;
411			list_vv = new TList;
412			list_st = new TList;
413			list_wjets = new TList;
414			list_others = new TList;
415			list_total = new TList;
416			list_comb = new TList;
417
418			data=0;
419			signal=0;
420			dyl=0;
421			dyc=0;
422			dyb=0;
423			ttbar=0;
424			vv=0;
425			st=0;
426			wjets=0;
427			others=0;
428
429			bcdata=0;
430			bcsignal=0;
431			bcdyl=0;
432			bcdyc=0;
433			bcdyb=0;
434			bcttbar=0;
435			bcvv=0;
436			bcst=0;
437			bcwjets=0;
438			bcothers=0;
439
440			err_data=0;
441			err_signal=0;
442			err_dyl=0;
443			err_dyc=0;
444			err_dyb=0;
445			err_ttbar=0;
446			err_vv=0;
447			err_st=0;
448			err_wjets=0;
449			err_others=0;
450
451			Run2011A_s = "Run2011A";
452			Run2011B_s = "Run2011B";
453			Signal_s = "ZH";
454			DYL_s = "DYL";
455			DYC_s = "DYC";
456			DYB_s = "DYB";
457			DY_s = "DY";
458			DYBOOSTED_s = "DYBOOSTED";
459			TTbar_s = "TTbar";
460			VV_s = "VV";
461			ST_s = "ST";
462			WJETS_s = "WJ";
463
464			}
465
466			private:
467
468			double lumi;
469			double fa;
470			double fb;
471			std::string name;
472
473			std::string Run2011A_s;
474			std::string Run2011B_s;
475			std::string Signal_s;
476			std::string DYL_s;
477			std::string DYC_s;
478			std::string DYB_s;
479			std::string DY_s;
480			std::string DYBOOSTED_s;
481			std::string TTbar_s;
482			std::string VV_s;
483			std::string ST_s;
484			std::string WJETS_s;
485
486			double data;
487			double signal;
488			double dyl;
489			double dyc;
490			double dyb;
491			double ttbar;
492			double vv;
493			double st;
494			double wjets;
495			double others;
496			double count_comb;
497
498			double bcdata;
499			double bcsignal;
500			double bcdyl;
501			double bcdyc;
502			double bcdyb;
503			double bcttbar;
504			double bcvv;
505			double bcst;
506			double bcwjets;
507			double bcothers;
508
509			TH1F * hdata;
510			TH1F * hsignal;
511			TH1F * hdyl;
512			TH1F * hdyc;
513			TH1F * hdyb;
514			TH1F * httbar;
515			TH1F * hvv;
516			TH1F * hst;
517			TH1F * hwjets;
518			TH1F * hothers;
519			TH1F * htotal;
520			TH1F * h_comb;
521
522			TList * list_data;
523			TList * list_signal;
524			TList * list_dyl;
525			TList * list_dyc;
526			TList * list_dyb;
527			TList * list_ttbar;
528			TList * list_vv;
529			TList * list_st;
530			TList * list_wjets;
531			TList * list_others;
532			TList * list_total;
533			TList * list_comb;
534
535			double err_data;
536			double err_signal;
537			double err_dyl;
538			double err_dyc;
539			double err_dyb;
540			double err_ttbar;
541			double err_vv;
542			double err_st;
543			double err_wjets;
544			double err_others;
545
546
547			};
548
549
550			#endif