VHbb/interface/samples.hpp

#ifndef samples_h
#define samples_h

#include <string>
#include <vector>
#include <TFile.h>
#include <TH1F.h>
#include <TTree.h>
#include <iostream>
#include <sstream>
#include "TMath.h"
#include <TCut.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

struct Sample {
  Sample();
Sample( float xs, std::string n, std::string f, int c, bool isdata, float datalumi=-1., float nEvents=-1 )
: xsec(xs),luminosity(datalumi),name(n),filename(f),color(c),data(isdata),f(0),nevents(nEvents) {}

  float lumi() {  if(data) return luminosity; else return numberOfEvents()/xsec; }
  float lumi(double fA, double fB) {  if(data) return luminosity; else return numberOfEvents(fA,fB)/xsec; }
  float scale(float l) { if(lumi()>0) return l/lumi(); else return 0;}
  float scale(float l, double *SF) { 
    std::string DYL("DYL");
    std::string DYC("DYC");
    std::string DYNoB("DYNoB");
    std::string DYB("DYB");
    std::string TTbar("TTbar");
    if(lumi()>0){
      if(name == DYL || name == DYC || name == DYNoB)
        return SF[0]*l/lumi();
      else if(name == DYB)
        return SF[1]*l/lumi();
      else if(name == TTbar) 
        return SF[2]*l/lumi(); 
      else
        return scale(l);
    }
    else 
      return 0; }

  float scale(float l, double fA, double fB) { if(lumi(fA,fB)>0) return l/lumi(fA,fB); else return 0;}
  float scale(float l, double fA, double fB, double *SF) { 
    std::string DYL("DYL");
    std::string DYC("DYC");
    std::string DYNoB("DYNoB");
    std::string DYB("DYB");
    std::string TTbar("TTbar");
    if(lumi(fA,fB)>0){
      if(name == DYL || name == DYC || name == DYNoB)
        return SF[0]*l/lumi(fA,fB);
      else if(name == TTbar) 
        return SF[1]*l/lumi(fA,fB); 
      else if(name == DYB)
        return SF[2]*l/lumi(fA,fB);
      else
        return scale(l,fA,fB);
    }
    else 
      return 0; }
  TFile * file() { if(f) return f; else return f=TFile::Open(filename.c_str());}
  float numberOfEvents( double fA, double fB ) 
  {
    if(data) return luminosity;
    double CountWithPU,CountWithPU2011B;
    double nOfEvents;
    if(nevents !=-1) return nevents;
    else
      {
       CountWithPU = ((TH1F*)file()->Get("CountWithPU"))->GetBinContent(1);
       CountWithPU2011B = ((TH1F*)file()->Get("CountWithPU2011B"))->GetBinContent(1);
       nOfEvents = fA*CountWithPU + fB*CountWithPU2011B;
       return nOfEvents;
      }
  }  

  float numberOfEvents() 
  {
    if(data) return luminosity;
    if(nevents !=-1) return nevents;
    else
      {
        //      return ((TH1F*)file()->Get("Count"))->GetBinContent(1);
        return  ((TH1F*)file()->Get("CountWithPU"))->GetBinContent(1);
      }
  }  

  void dump(float l)
  {
    std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi()/1000 << "/fb \t& " << scale(l) << "\t& " << numberOfEvents() << std::endl;
    //     std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi()/1000 << "/fb \t& " << scale(l) << std::endl;
    //     std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi()/1000 << std::endl;
  }

  void dump(float l, double fa, double fb)
  {
    std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi(fa,fb)/1000 << "/fb \t& " << scale(l,fa,fb) << "\t& " << numberOfEvents(fa,fb) << std::endl;
    //     std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi()/1000 << "/fb \t& " << scale(l) << std::endl;
    //     std::cout << name <<  "\t& " << xsec << "\t& " <<  lumi()/1000 << std::endl;
  }
  
  float bareCount(TCut iCut){
    float NpassedEvents = 0;
    file()->cd();
    NpassedEvents = ((TTree*)file()->Get("tree"))->Draw("H.dPhi", iCut, "goff");
    return NpassedEvents;
  }
  
  float bareCount(TCut iWeight, TCut iCut){
    float NpassedEvents = 0;
    file()->cd();
    NpassedEvents = ((TTree*)file()->Get("tree"))->Draw("H.dPhi", iWeight * iCut, "goff");
    return NpassedEvents;
  }

  float bareCountError(TCut iCut){ bareCountError_ = 0; if(bareCount(iCut)>0) bareCountError_ = TMath::Sqrt(bareCount(iCut)); return bareCountError_; }
  
  float count(TCut iCut){
    float NweightedEvents = 0;
    if(data)
      return bareCount(iCut);
    else{
      TH1F* histo;
      ((TTree*)file()->Get("tree"))->Draw("H.dPhi>>histo", iCut, "goff");
      histo = (TH1F*)gDirectory->Get("histo");
      NweightedEvents = histo->Integral()/lumi();
      return NweightedEvents;
    }
  }

  float count(TCut iWeight, TCut iCut){
    float NweightedEvents = 0;
    if(data)
      return bareCount(iCut);
    else{
      TH1F* histo;
      ((TTree*)file()->Get("tree"))->Draw("H.dPhi>>histo", iWeight*iCut, "goff");
      histo = (TH1F*)gDirectory->Get("histo");
      NweightedEvents = histo->Integral()/lumi();
      return NweightedEvents;
    }
  }

  float count(TCut iCut, double fA, double fB){
    float NweightedEvents = 0;
    if(data)
      return bareCount(iCut);
    else{
      TH1F* histo;
      ((TTree*)file()->Get("tree"))->Draw("H.dPhi>>histo", iCut, "goff");
      histo = (TH1F*)gDirectory->Get("histo");
      NweightedEvents = histo->Integral()/lumi(fA, fB);
      return NweightedEvents;
    }
  }

  float count(TCut iWeight, TCut iCut, double fA, double fB){
    float NweightedEvents = 0;
    if(data)
      return bareCount(iCut);
    else{
      TH1F* histo;
      std::stringstream oss;
      srand( time(NULL) );
      int num = rand() % 10000 + 1;
      oss << num;
      std::string var("H.dPhi>>");
      std::string histoName("histo");
      histoName=histoName+name+oss.str();
      ((TTree*)file()->Get("tree"))->Draw((var+histoName).c_str(), iWeight*iCut, "goff");
      histo = (TH1F*)gDirectory->Get(histoName.c_str());
      NweightedEvents = histo->Integral()/lumi(fA, fB);
      return NweightedEvents;
    }
  }

  float countError(TCut iCut){ countError_ = count(iCut)/bareCountError(iCut); return countError_; }

  float countError(TCut iWeight, TCut iCut, double fA, double fB){ countError_ = 0; if(bareCountError(iCut)>0) countError_ = count(iWeight,iCut,fA, fB)/bareCountError(iCut); return countError_; }

  float countError(TCut iCut, float systematicError ){
    float totalError=0;
    if(systematicError < 1)
      totalError = TMath::Sqrt( TMath::Power(count(iCut)/bareCountError(iCut),2) + TMath::Power( systematicError*count(iCut),2) );
    else
      std::cout << "Systematic error has to be espressed in %" << std::endl;
    return totalError;
  }
  
  float addSysError( float error, float systematicError ){ return TMath::Sqrt( TMath::Power(error,2) + TMath::Power( systematicError,2) ); }
  
  void setCountError( float newCountError ){
    countError_ = newCountError;
  }
  

  float countError_;
  float bareCountError_;
  float nevents;
  float xsec;
  float luminosity;
  std::string name;
  std::string filename;
  int color;
  bool data;
  TFile * f;
  
}; 

#endif
Revision:	1.3
Committed:	Thu Oct 4 09:54:36 2012 UTC (12 years, 7 months ago) by bortigno
Branch:	MAIN
CVS Tags:	lhcp_UnblindFix, hcp_Unblind, lhcp_11April, LHCP_PreAppFixAfterFreeze, LHCP_PreAppFreeze, workingVersionAfterHCP, hcpApproval, hcpPreApp, HEAD
Changes since 1.2:	+3 -3 lines
Log Message:	correct SF order in scale(lumi,SF) function
#	User	Rev	Content
1	bortigno	1.1	#ifndef samples_h
2			#define samples_h
3
4			#include <string>
5			#include <vector>
6			#include <TFile.h>
7			#include <TH1F.h>
8			#include <TTree.h>
9			#include <iostream>
10			#include <sstream>
11			#include "TMath.h"
12			#include <TCut.h>
13			#include <stdio.h>
14			#include <stdlib.h>
15			#include <time.h>
16
17			struct Sample {
18			Sample();
19			Sample( float xs, std::string n, std::string f, int c, bool isdata, float datalumi=-1., float nEvents=-1 )
20			: xsec(xs),luminosity(datalumi),name(n),filename(f),color(c),data(isdata),f(0),nevents(nEvents) {}
21
22			float lumi() { if(data) return luminosity; else return numberOfEvents()/xsec; }
23			float lumi(double fA, double fB) { if(data) return luminosity; else return numberOfEvents(fA,fB)/xsec; }
24			float scale(float l) { if(lumi()>0) return l/lumi(); else return 0;}
25	bortigno	1.2	float scale(float l, double *SF) {
26			std::string DYL("DYL");
27			std::string DYC("DYC");
28			std::string DYNoB("DYNoB");
29			std::string DYB("DYB");
30			std::string TTbar("TTbar");
31			if(lumi()>0){
32			if(name == DYL \|\| name == DYC \|\| name == DYNoB)
33			return SF[0]*l/lumi();
34	bortigno	1.3	else if(name == DYB)
35			return SF[1]*l/lumi();
36	bortigno	1.2	else if(name == TTbar)
37	bortigno	1.3	return SF[2]*l/lumi();
38	bortigno	1.2	else
39			return scale(l);
40			}
41			else
42			return 0; }
43
44	bortigno	1.1	float scale(float l, double fA, double fB) { if(lumi(fA,fB)>0) return l/lumi(fA,fB); else return 0;}
45			float scale(float l, double fA, double fB, double *SF) {
46			std::string DYL("DYL");
47			std::string DYC("DYC");
48			std::string DYNoB("DYNoB");
49			std::string DYB("DYB");
50			std::string TTbar("TTbar");
51			if(lumi(fA,fB)>0){
52			if(name == DYL \|\| name == DYC \|\| name == DYNoB)
53			return SF[0]*l/lumi(fA,fB);
54			else if(name == TTbar)
55			return SF[1]*l/lumi(fA,fB);
56			else if(name == DYB)
57			return SF[2]*l/lumi(fA,fB);
58			else
59			return scale(l,fA,fB);
60			}
61			else
62			return 0; }
63			TFile * file() { if(f) return f; else return f=TFile::Open(filename.c_str());}
64			float numberOfEvents( double fA, double fB )
65			{
66			if(data) return luminosity;
67			double CountWithPU,CountWithPU2011B;
68			double nOfEvents;
69			if(nevents !=-1) return nevents;
70			else
71			{
72			CountWithPU = ((TH1F*)file()->Get("CountWithPU"))->GetBinContent(1);
73			CountWithPU2011B = ((TH1F*)file()->Get("CountWithPU2011B"))->GetBinContent(1);
74			nOfEvents = fACountWithPU + fBCountWithPU2011B;
75			return nOfEvents;
76			}
77			}
78
79			float numberOfEvents()
80			{
81			if(data) return luminosity;
82			if(nevents !=-1) return nevents;
83			else
84			{
85	bortigno	1.2	// return ((TH1F*)file()->Get("Count"))->GetBinContent(1);
86			return ((TH1F*)file()->Get("CountWithPU"))->GetBinContent(1);
87	bortigno	1.1	}
88			}
89
90			void dump(float l)
91			{
92			std::cout << name << "\t& " << xsec << "\t& " << lumi()/1000 << "/fb \t& " << scale(l) << "\t& " << numberOfEvents() << std::endl;
93			// std::cout << name << "\t& " << xsec << "\t& " << lumi()/1000 << "/fb \t& " << scale(l) << std::endl;
94			// std::cout << name << "\t& " << xsec << "\t& " << lumi()/1000 << std::endl;
95			}
96
97			void dump(float l, double fa, double fb)
98			{
99			std::cout << name << "\t& " << xsec << "\t& " << lumi(fa,fb)/1000 << "/fb \t& " << scale(l,fa,fb) << "\t& " << numberOfEvents(fa,fb) << std::endl;
100			// std::cout << name << "\t& " << xsec << "\t& " << lumi()/1000 << "/fb \t& " << scale(l) << std::endl;
101			// std::cout << name << "\t& " << xsec << "\t& " << lumi()/1000 << std::endl;
102			}
103
104			float bareCount(TCut iCut){
105			float NpassedEvents = 0;
106			file()->cd();
107			NpassedEvents = ((TTree*)file()->Get("tree"))->Draw("H.dPhi", iCut, "goff");
108			return NpassedEvents;
109			}
110
111			float bareCount(TCut iWeight, TCut iCut){
112			float NpassedEvents = 0;
113			file()->cd();
114			NpassedEvents = ((TTree)file()->Get("tree"))->Draw("H.dPhi", iWeight iCut, "goff");
115			return NpassedEvents;
116			}
117
118			float bareCountError(TCut iCut){ bareCountError_ = 0; if(bareCount(iCut)>0) bareCountError_ = TMath::Sqrt(bareCount(iCut)); return bareCountError_; }
119
120			float count(TCut iCut){
121			float NweightedEvents = 0;
122			if(data)
123			return bareCount(iCut);
124			else{
125			TH1F* histo;
126			((TTree*)file()->Get("tree"))->Draw("H.dPhi>>histo", iCut, "goff");
127			histo = (TH1F*)gDirectory->Get("histo");
128			NweightedEvents = histo->Integral()/lumi();
129			return NweightedEvents;
130			}
131			}
132
133			float count(TCut iWeight, TCut iCut){
134			float NweightedEvents = 0;
135			if(data)
136			return bareCount(iCut);
137			else{
138			TH1F* histo;
139			((TTree)file()->Get("tree"))->Draw("H.dPhi>>histo", iWeightiCut, "goff");
140			histo = (TH1F*)gDirectory->Get("histo");
141			NweightedEvents = histo->Integral()/lumi();
142			return NweightedEvents;
143			}
144			}
145
146			float count(TCut iCut, double fA, double fB){
147			float NweightedEvents = 0;
148			if(data)
149			return bareCount(iCut);
150			else{
151			TH1F* histo;
152			((TTree*)file()->Get("tree"))->Draw("H.dPhi>>histo", iCut, "goff");
153			histo = (TH1F*)gDirectory->Get("histo");
154			NweightedEvents = histo->Integral()/lumi(fA, fB);
155			return NweightedEvents;
156			}
157			}
158
159			float count(TCut iWeight, TCut iCut, double fA, double fB){
160			float NweightedEvents = 0;
161			if(data)
162			return bareCount(iCut);
163			else{
164			TH1F* histo;
165			std::stringstream oss;
166			srand( time(NULL) );
167			int num = rand() % 10000 + 1;
168			oss << num;
169			std::string var("H.dPhi>>");
170			std::string histoName("histo");
171			histoName=histoName+name+oss.str();
172			((TTree)file()->Get("tree"))->Draw((var+histoName).c_str(), iWeightiCut, "goff");
173			histo = (TH1F*)gDirectory->Get(histoName.c_str());
174			NweightedEvents = histo->Integral()/lumi(fA, fB);
175			return NweightedEvents;
176			}
177			}
178
179			float countError(TCut iCut){ countError_ = count(iCut)/bareCountError(iCut); return countError_; }
180
181			float countError(TCut iWeight, TCut iCut, double fA, double fB){ countError_ = 0; if(bareCountError(iCut)>0) countError_ = count(iWeight,iCut,fA, fB)/bareCountError(iCut); return countError_; }
182
183			float countError(TCut iCut, float systematicError ){
184			float totalError=0;
185			if(systematicError < 1)
186			totalError = TMath::Sqrt( TMath::Power(count(iCut)/bareCountError(iCut),2) + TMath::Power( systematicError*count(iCut),2) );
187			else
188			std::cout << "Systematic error has to be espressed in %" << std::endl;
189			return totalError;
190			}
191
192			float addSysError( float error, float systematicError ){ return TMath::Sqrt( TMath::Power(error,2) + TMath::Power( systematicError,2) ); }
193
194			void setCountError( float newCountError ){
195			countError_ = newCountError;
196			}
197
198
199			float countError_;
200			float bareCountError_;
201			float nevents;
202			float xsec;
203			float luminosity;
204			std::string name;
205			std::string filename;
206			int color;
207			bool data;
208			TFile * f;
209
210			};
211
212			#endif