VHbbDataFormats/interface/TriggerWeight.h

#ifndef TRIGGERWEIGHT_H
#define TRIGGERWEIGHT_H

#include "FWCore/ParameterSet/interface/ProcessDesc.h"
#include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerZnunuCurve.h"
#include "VHbbAnalysis/VHbbDataFormats/interface/MultiThresholdEfficiency.h"
#include <TH1F.h>
#include <TF1.h>
#include <TFile.h>
#include <TTree.h>
#include <iostream>

class TriggerWeight
{
public:  
  TriggerWeight(const edm::ParameterSet& ana) : combiner2Thr(2)
  {
   tscaleHLTmu=openFile(ana,"hltMuFileName");
   tscaleIDmu=openFile(ana,"idMuFileName");
   tscaleHLTele1=openFile(ana,"hltEle1FileName");
   tscaleHLTele2=openFile(ana,"hltEle2FileName");
   tscaleID80Ele=openFile(ana,"idEle80FileName");
   tscaleID95Ele=openFile(ana,"idEle95FileName");
   tscaleHLTeleJet1=openFile(ana,"hltJetEle1FileName");
   tscaleHLTeleJet2=openFile(ana,"hltJetEle2FileName");
   tscaleRecoEle=openFile(ana,"recoEleFileName");
//   tscalePFMHTele=openFile(ana,"hltPFMHTEleFileName");
   tscaleSingleEleMay=openFile(ana,"hltSingleEleMayFileName");
   tscaleSingleEleV4=openFile(ana,"hltSingleEleV4FileName");
   tscaleHLTmuOr30=openFile(ana,"hltMuOr30FileName");

   if(tscaleHLTmu == 0 || tscaleIDmu == 0) 
    {
      std::cout << "ERROR: cannot load Muon Trigger efficiencies" << std::endl;
    } 
   

  }
 
  TTree * openFile(const edm::ParameterSet& ana, const char * name)
  {
   TFile *hltMuFile = new TFile (ana.getParameter<std::string> (name).c_str(),"read");
   if(hltMuFile)   return (TTree*) hltMuFile->Get("tree");
   else return 0;
  }

  std::pair<float,float> efficiencyFromPtEta(float pt1, float eta1, TTree *t)
  {
     float s1 = 1.,err=1.;
     std::pair<float,float> r(s1,err);
     if(!t) return r;
     float ptMin,ptMax,etaMin,etaMax,scale,error;
     int count = 0;
     t->SetBranchAddress("ptMin",&ptMin);
     t->SetBranchAddress("ptMax",&ptMax);
     t->SetBranchAddress("etaMin",&etaMin);
     t->SetBranchAddress("etaMax",&etaMax);
     t->SetBranchAddress("scale",&scale);
     t->SetBranchAddress("error",&error);

    for(int jentry = 0; jentry < t->GetEntries(); jentry++)
      {
        t->GetEntry(jentry);
        if((pt1 > ptMin) && (pt1 < ptMax) && (eta1 > etaMin) && (eta1 < etaMax))
          {
            s1 = scale;
            err=error;
            count++;
          }
      }

    if(count == 0 || s1 == 0)
      {
        return r;
      }

     r.first=s1;
     r.second = err;
    return (r);
  }
 
  float scaleMuIsoHLT(float pt1, float eta1)
  {
    return efficiencyFromPtEta(pt1,eta1,tscaleHLTmu).first;
  }
  
  
  float scaleMuID(float pt1, float eta1)
  {
    return efficiencyFromPtEta(pt1,eta1,tscaleIDmu).first;
  }

double  scaleMetHLT( double met){

    float s1 = 1;
    TF1 * f = new TF1 ("f",TriggerZnunuCurve::trigMet, 0,99999, 0, "triggerZnunuCurve"  );
    
    s1 = f->Eval(met);
    
    return (s1);
    
  }
  

double scaleDoubleEle17Ele8( std::vector<float> pt, std::vector<float> eta )
{
   std::vector< std::vector<float> > allEleWithEffs;
for(unsigned int j=0; j< pt.size(); j++)
 {
  std::vector<float> thisEleEffs;
  thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele1).first);
  thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele2).first);
  allEleWithEffs.push_back(thisEleEffs);
 }

  return   combiner2Thr.weight<Trigger1High2Loose>(allEleWithEffs);

}

double scaleSingleEleMay( float pt, float eta){    return efficiencyFromPtEta(pt,eta,tscaleSingleEleMay).first;}
double scaleSingleEleV4( float pt, float eta){   return efficiencyFromPtEta(pt,eta,tscaleSingleEleV4).first; }
double scaleID80Ele( float pt, float eta) {      return efficiencyFromPtEta(pt,eta,tscaleID80Ele).first; }
double scaleID95Ele( float pt, float eta) {      return efficiencyFromPtEta(pt,eta,tscaleID95Ele).first; }
double scaleRecoEle( float pt, float eta){     return efficiencyFromPtEta(pt,eta,tscaleRecoEle).first; }
double scalePFMHTEle( float MetPFPt){
    double weightPFMHTrigger=0.;

    //FIXME: read from file
    if(MetPFPt>0. && MetPFPt<5.) weightPFMHTrigger=0.3834;
    if(MetPFPt>5. && MetPFPt<10.) weightPFMHTrigger=0.4493;
    if(MetPFPt>10. && MetPFPt<15.) weightPFMHTrigger=0.5676;
    if(MetPFPt>15. && MetPFPt<20.) weightPFMHTrigger=0.6474;
    if(MetPFPt>20. && MetPFPt<25.) weightPFMHTrigger=0.7695;
    if(MetPFPt>25. && MetPFPt<30.) weightPFMHTrigger=0.8936;
    if(MetPFPt>30. && MetPFPt<35.) weightPFMHTrigger=0.9304;
    if(MetPFPt>35. && MetPFPt<40.) weightPFMHTrigger=0.9620;
    if(MetPFPt>40. && MetPFPt<45.) weightPFMHTrigger=0.9894;
    if(MetPFPt>45. && MetPFPt<50.) weightPFMHTrigger=0.9863;
    if(MetPFPt>50. && MetPFPt<60.) weightPFMHTrigger=0.9978;
    if(MetPFPt>60.) weightPFMHTrigger=1;
    return weightPFMHTrigger;
}

double scaleJet30Jet25( std::vector<float> pt, std::vector<float> eta)
{

  std::vector< std::vector<float> > allJetsWithEffs;
for(unsigned int j=0; j< pt.size(); j++)
 {
  std::vector<float> thisJetEffs;
  thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet1).first);
  thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet2).first);
  allJetsWithEffs.push_back(thisJetEffs);
 }

  return   combiner2Thr.weight<Trigger1High2Loose>(allJetsWithEffs);
}

float scaleMuOr30IsoHLT(float pt1, float eta1)
{
   return efficiencyFromPtEta(pt1,eta1,tscaleHLTmuOr30).first;
}


private:
  TTree * tscaleHLTele1;
  TTree * tscaleHLTele2;
  TTree * tscaleHLTeleJet1;
  TTree * tscaleHLTeleJet2;
  TTree * tscaleID80Ele;
  TTree * tscaleID95Ele;
  TTree * tscaleRecoEle;
  TTree * tscaleHLTmuOr30;

//  TTree * tscalePFMHTele;
  TTree * tscaleSingleEleMay;
  TTree * tscaleSingleEleV4;

  TTree * tscaleHLTmu;
  TTree * tscaleIDmu;
  MultiThresholdEfficiency combiner2Thr; 
};

#endif
Revision:	1.6
Committed:	Tue Oct 18 10:23:09 2011 UTC (13 years, 6 months ago) by arizzi
Content type:	text/plain
Branch:	MAIN
CVS Tags:	AR_step2_oct19
Changes since 1.5:	+7 -0 lines
Log Message:	changes asked by david. Add jet_e, btagSF for 1Tight+1AntiCustom, new trigger, new PU weights
#	User	Rev	Content
1	arizzi	1.1	#ifndef TRIGGERWEIGHT_H
2			#define TRIGGERWEIGHT_H
3
4			#include "FWCore/ParameterSet/interface/ProcessDesc.h"
5			#include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
6	arizzi	1.4	#include "VHbbAnalysis/VHbbDataFormats/interface/TriggerZnunuCurve.h"
7			#include "VHbbAnalysis/VHbbDataFormats/interface/MultiThresholdEfficiency.h"
8	arizzi	1.1	#include <TH1F.h>
9	arizzi	1.3	#include <TF1.h>
10	arizzi	1.1	#include <TFile.h>
11			#include <TTree.h>
12			#include <iostream>
13
14			class TriggerWeight
15			{
16			public:
17	arizzi	1.5	TriggerWeight(const edm::ParameterSet& ana) : combiner2Thr(2)
18	arizzi	1.1	{
19	arizzi	1.4	tscaleHLTmu=openFile(ana,"hltMuFileName");
20			tscaleIDmu=openFile(ana,"idMuFileName");
21			tscaleHLTele1=openFile(ana,"hltEle1FileName");
22			tscaleHLTele2=openFile(ana,"hltEle2FileName");
23	arizzi	1.5	tscaleID80Ele=openFile(ana,"idEle80FileName");
24			tscaleID95Ele=openFile(ana,"idEle95FileName");
25	arizzi	1.4	tscaleHLTeleJet1=openFile(ana,"hltJetEle1FileName");
26			tscaleHLTeleJet2=openFile(ana,"hltJetEle2FileName");
27	arizzi	1.5	tscaleRecoEle=openFile(ana,"recoEleFileName");
28			// tscalePFMHTele=openFile(ana,"hltPFMHTEleFileName");
29			tscaleSingleEleMay=openFile(ana,"hltSingleEleMayFileName");
30			tscaleSingleEleV4=openFile(ana,"hltSingleEleV4FileName");
31	arizzi	1.6	tscaleHLTmuOr30=openFile(ana,"hltMuOr30FileName");
32	arizzi	1.4
33	arizzi	1.1	if(tscaleHLTmu == 0 \|\| tscaleIDmu == 0)
34			{
35			std::cout << "ERROR: cannot load Muon Trigger efficiencies" << std::endl;
36			}
37	arizzi	1.4
38
39	arizzi	1.1
40			}
41	arizzi	1.4
42			TTree * openFile(const edm::ParameterSet& ana, const char * name)
43			{
44	arizzi	1.5	TFile *hltMuFile = new TFile (ana.getParameter<std::string> (name).c_str(),"read");
45	arizzi	1.4	if(hltMuFile) return (TTree*) hltMuFile->Get("tree");
46			else return 0;
47			}
48
49			std::pair<float,float> efficiencyFromPtEta(float pt1, float eta1, TTree *t)
50			{
51			float s1 = 1.,err=1.;
52			std::pair<float,float> r(s1,err);
53			if(!t) return r;
54			float ptMin,ptMax,etaMin,etaMax,scale,error;
55			int count = 0;
56			t->SetBranchAddress("ptMin",&ptMin);
57			t->SetBranchAddress("ptMax",&ptMax);
58			t->SetBranchAddress("etaMin",&etaMin);
59			t->SetBranchAddress("etaMax",&etaMax);
60			t->SetBranchAddress("scale",&scale);
61			t->SetBranchAddress("error",&error);
62
63			for(int jentry = 0; jentry < t->GetEntries(); jentry++)
64			{
65			t->GetEntry(jentry);
66			if((pt1 > ptMin) && (pt1 < ptMax) && (eta1 > etaMin) && (eta1 < etaMax))
67			{
68			s1 = scale;
69			err=error;
70			count++;
71			}
72			}
73
74			if(count == 0 \|\| s1 == 0)
75			{
76			return r;
77			}
78
79			r.first=s1;
80			r.second = err;
81			return (r);
82			}
83	arizzi	1.1
84			float scaleMuIsoHLT(float pt1, float eta1)
85			{
86	arizzi	1.4	return efficiencyFromPtEta(pt1,eta1,tscaleHLTmu).first;
87	arizzi	1.1	}
88
89
90
91			float scaleMuID(float pt1, float eta1)
92			{
93	arizzi	1.4	return efficiencyFromPtEta(pt1,eta1,tscaleIDmu).first;
94	arizzi	1.1	}
95
96	arizzi	1.3	double scaleMetHLT( double met){
97
98			float s1 = 1;
99			TF1 * f = new TF1 ("f",TriggerZnunuCurve::trigMet, 0,99999, 0, "triggerZnunuCurve" );
100
101			s1 = f->Eval(met);
102
103			return (s1);
104
105			}
106
107
108	arizzi	1.5	double scaleDoubleEle17Ele8( std::vector<float> pt, std::vector<float> eta )
109			{
110			std::vector< std::vector<float> > allEleWithEffs;
111			for(unsigned int j=0; j< pt.size(); j++)
112			{
113			std::vector<float> thisEleEffs;
114			thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele1).first);
115			thisEleEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTele2).first);
116			allEleWithEffs.push_back(thisEleEffs);
117			}
118
119			return combiner2Thr.weight<Trigger1High2Loose>(allEleWithEffs);
120
121			}
122
123			double scaleSingleEleMay( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEleMay).first;}
124			double scaleSingleEleV4( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleSingleEleV4).first; }
125			double scaleID80Ele( float pt, float eta) { return efficiencyFromPtEta(pt,eta,tscaleID80Ele).first; }
126			double scaleID95Ele( float pt, float eta) { return efficiencyFromPtEta(pt,eta,tscaleID95Ele).first; }
127			double scaleRecoEle( float pt, float eta){ return efficiencyFromPtEta(pt,eta,tscaleRecoEle).first; }
128			double scalePFMHTEle( float MetPFPt){
129			double weightPFMHTrigger=0.;
130
131			//FIXME: read from file
132			if(MetPFPt>0. && MetPFPt<5.) weightPFMHTrigger=0.3834;
133			if(MetPFPt>5. && MetPFPt<10.) weightPFMHTrigger=0.4493;
134			if(MetPFPt>10. && MetPFPt<15.) weightPFMHTrigger=0.5676;
135			if(MetPFPt>15. && MetPFPt<20.) weightPFMHTrigger=0.6474;
136			if(MetPFPt>20. && MetPFPt<25.) weightPFMHTrigger=0.7695;
137			if(MetPFPt>25. && MetPFPt<30.) weightPFMHTrigger=0.8936;
138			if(MetPFPt>30. && MetPFPt<35.) weightPFMHTrigger=0.9304;
139			if(MetPFPt>35. && MetPFPt<40.) weightPFMHTrigger=0.9620;
140			if(MetPFPt>40. && MetPFPt<45.) weightPFMHTrigger=0.9894;
141			if(MetPFPt>45. && MetPFPt<50.) weightPFMHTrigger=0.9863;
142			if(MetPFPt>50. && MetPFPt<60.) weightPFMHTrigger=0.9978;
143			if(MetPFPt>60.) weightPFMHTrigger=1;
144			return weightPFMHTrigger;
145			}
146
147			double scaleJet30Jet25( std::vector<float> pt, std::vector<float> eta)
148	arizzi	1.4	{
149	arizzi	1.5
150			std::vector< std::vector<float> > allJetsWithEffs;
151			for(unsigned int j=0; j< pt.size(); j++)
152			{
153			std::vector<float> thisJetEffs;
154			thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet1).first);
155			thisJetEffs.push_back(efficiencyFromPtEta(pt[j],eta[j],tscaleHLTeleJet2).first);
156			allJetsWithEffs.push_back(thisJetEffs);
157			}
158
159			return combiner2Thr.weight<Trigger1High2Loose>(allJetsWithEffs);
160	arizzi	1.4	}
161	arizzi	1.3
162	arizzi	1.6	float scaleMuOr30IsoHLT(float pt1, float eta1)
163			{
164			return efficiencyFromPtEta(pt1,eta1,tscaleHLTmuOr30).first;
165			}
166	arizzi	1.3
167
168	arizzi	1.1	private:
169	arizzi	1.4	TTree * tscaleHLTele1;
170			TTree * tscaleHLTele2;
171			TTree * tscaleHLTeleJet1;
172			TTree * tscaleHLTeleJet2;
173	arizzi	1.5	TTree * tscaleID80Ele;
174			TTree * tscaleID95Ele;
175			TTree * tscaleRecoEle;
176	arizzi	1.6	TTree * tscaleHLTmuOr30;
177
178	arizzi	1.5	// TTree * tscalePFMHTele;
179			TTree * tscaleSingleEleMay;
180			TTree * tscaleSingleEleV4;
181	arizzi	1.4
182	arizzi	1.1	TTree * tscaleHLTmu;
183			TTree * tscaleIDmu;
184	arizzi	1.4	MultiThresholdEfficiency combiner2Thr;
185	arizzi	1.1	};
186
187			#endif