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
#	Content
1	#ifndef CONTROLREGION_H
2	#define CONTROLREGION_H
3
4	#include "samples.hpp"
5	#include "ntupleReader.hpp"
6	#include "CutsAndHistos.h"
7	#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