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
#	Content
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	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	else if(name == DYB)
35	return SF[1]*l/lumi();
36	else if(name == TTbar)
37	return SF[2]*l/lumi();
38	else
39	return scale(l);
40	}
41	else
42	return 0; }
43
44	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	// return ((TH1F*)file()->Get("Count"))->GetBinContent(1);
86	return ((TH1F*)file()->Get("CountWithPU"))->GetBinContent(1);
87	}
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