Plotting/Modules/Systematics.C

#include <iostream>
#include <vector>
#include <sys/stat.h>

#include <TMath.h>
#include <TColor.h>
#include <TPaveText.h>
#include <TRandom.h>
#include <TF1.h>

#ifndef SampleClassLoaded
#include "ActiveSamples.C"
#endif

#ifndef Verbosity
#define Verbosity 0
#endif

#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TCut.h>
#include <TMath.h>
#include <TLine.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TF1.h>


Int_t nBins    =  100;
Float_t jzbMin = -207;
Float_t jzbMax =  243;
Float_t jzbSel =  100;
int iplot=0;
int verbose=0;
string geqleq;
string mcjzbexpression;

TString geq_or_leq() {
  if(geqleq=="geq") return TString(">=");
  if(geqleq=="leq") return TString("<=");
  return TString("GEQ_OR_LEQ_ERROR");
}

//______________________________________________________________________________
Double_t Interpolate(Double_t x, TH1 *histo)
{
   // Given a point x, approximates the value via linear interpolation
   // based on the two nearest bin centers
   // Andy Mastbaum 10/21/08
   // in newer ROOT versions but not in the one I have so I had to work around that ... 

   Int_t xbin = histo->FindBin(x);
   Double_t x0,x1,y0,y1;

   if(x<=histo->GetBinCenter(1)) {
      return histo->GetBinContent(1);
   } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
      return histo->GetBinContent(histo->GetNbinsX());
   } else {
      if(x<=histo->GetBinCenter(xbin)) {
         y0 = histo->GetBinContent(xbin-1);
         x0 = histo->GetBinCenter(xbin-1);
         y1 = histo->GetBinContent(xbin);
         x1 = histo->GetBinCenter(xbin);
      } else {
         y0 = histo->GetBinContent(xbin);
         x0 = histo->GetBinCenter(xbin);
         y1 = histo->GetBinContent(xbin+1);
         x1 = histo->GetBinCenter(xbin+1);
      }
      return y0 + (x-x0)*((y1-y0)/(x1-x0));
   }
}


//____________________________________________________________________________________
// Efficiency plot
TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
        iplot++;
        int count=iplot;
        // Define new histogram
        char hname[30]; sprintf(hname,"hJzbEff%d",count);
        TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
                                                         nBins,jzbMin,jzbMax);
        Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
        
        events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
        Float_t maxEff = events->GetSelectedRows();
        if(verbose>0) std::cout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
        
        if(verbose>0) std::cout <<  "JZB max = " << jzbMax << std::endl;
        // Loop over steps to get efficiency curve
        char cut[256];
        for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
                sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
                events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
                Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
                //     std::cout << "COUCOU " << __LINE__ << std::endl;
                hJzbEff->SetBinContent(iBin+1,eff);
                hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
        }
        return hJzbEff;
        
        
}


//________________________________________________________________________________________
// Pile-up efficiency
float pileup(TTree *events, string informalname, Float_t myJzbMax = 140. ) {
        nBins = 16;
        jzbMax = myJzbMax;
        
        // Acceptance cuts
        TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
        
        TH1F* hLM4 = plotEff(events,kbase,informalname);
        hLM4->SetMinimum(0.);
        
        // Nominal function
        TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
        func->SetParameter(0,0.03);
        func->SetParameter(1,0.);  
        hLM4->Fit(func,"Q");
        
        // Pimped-up function
        TF1* funcUp = (TF1*)func->Clone();
        funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
        std::cout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel) 
        << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
        
        return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
        
}

//____________________________________________________________________________________
// Total selection efficiency (MC)
void MCefficiency(TTree *events,float &res, float &reserr,string mcjzb) {
        
        char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
        // All acceptance cuts at gen. level
        TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
        // Corresponding reco. cuts
        TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
        
        events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
        Float_t sel = events->GetSelectedRows();
        events->Draw(mcjzbexpression.c_str(),kbase,"goff");
        Float_t tot = events->GetSelectedRows();
        
        res=sel/tot;
        reserr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
        std::cout << "  MC efficiency: " << res << "+-" << reserr << std::endl;
}

float JZBefficiency(TTree *events, string informalname) {
        TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
        TH1F* hLM4 = plotEff(events,kbase,informalname);
        Int_t bin = hLM4->FindBin(jzbSel); // To get the error
        std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
        std::cout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
        return -1;
}

//________________________________________________________________________
// Effect of energy scale on efficiency
void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
        TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0");
        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
        TCut nJets("pfJetGoodNum>2");
        stringstream down,up;
        down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
        up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
        
        TCut nJetsP(up.str().c_str());
        TCut nJetsM(down.str().c_str());
        
        if ( jzbSelection>0 ) jzbSel = jzbSelection;
        
        if ( !(plotName.Length()>1) ) plotName = informalname;
        
        nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
        TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
        
        TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
        
        TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
        
        // Dump some information
        Float_t eff  = Interpolate(jzbSel,hist);
        Float_t effp = Interpolate(jzbSel,histp);
        Float_t effm = Interpolate(jzbSel,histm);
        std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
        std::cout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl; 
        std::cout << "    central:  " << eff << std::endl; 
        std::cout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl; 
        jesup=(effp-eff)/eff*100.;
        jesdown=(effm-eff)/eff*100.;
}

//________________________________________________________________________
// Effect of energy scale on JZB efficiency
void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
        
        TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&genNjets>2");
        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
        
        nBins =    50;
        jzbMin =   0.5*jzbSel;
        jzbMax =   2.0*jzbSel;
        
        TH1F* hist = plotEff(events,kbase&&ksel,informalname);
        
        // Dump some information
        Float_t eff  = Interpolate(jzbSel,hist);
        Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
        Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
        std::cout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl; 
        std::cout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl; 
        std::cout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
        up=((effp-eff)/eff)*100;
        down=((effm-eff)/eff)*100;
}

//________________________________________________________________________
// JZB response (true/reco. vs. true)
void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
        
        jzbMin = 20;
        TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&genNjets>2");
        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
        
        TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", 
                                                                          nPeriods, jzbMin, myJzbMax, "" );
        
        if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
        else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
        
        hJzbResp->SetMaximum(1.2);
        hJzbResp->SetMinimum(0.2);
        hJzbResp->Fit("pol0","Q");
        TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
        cout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
}


void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &uncertainties,string mcjzb,string datajzb) {
  
  float JetEnergyScaleUncert=0.1;
  float JZBScaleUncert=0.1;
  mcjzbexpression=mcjzb;
  
  float triggereff=4;//percent!
  cout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
  float leptonseleff=2;//percent!
  cout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
  
  float mceff,mcefferr;
  cout << "MC efficiencies:" << endl;
  MCefficiency(events,mceff,mcefferr,mcjzb);
  JZBefficiency(events,informalname);
  
  std::cout << "Jet energy scale: " << std::endl;
  float jesup,jesdown;
  JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
  
  std::cout << "JZB scale: " << std::endl;
  float scaleup,scaledown,scalesyst;
  doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
  
  std::cout << "JZB response: " << std::endl;
  JZBresponse(events);

  std::cout << "Pileup: " << std::endl;
  float resolution=pileup(events,informalname);
  
  cout << "_______________________________________________" << endl;
  cout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
  cout << "Trigger efficiency: " << triggereff << endl;
  cout << "Lepton Sel Eff: " << leptonseleff << endl;
  cout << "For JZB>" << jzbSel << endl;
  cout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
  cout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
  cout << "Resolution : " << resolution << endl;
  
  vector<float> uncert;
  uncert.push_back(jzbSel);
  uncert.push_back(triggereff);
  uncert.push_back(leptonseleff);
  uncert.push_back(0.5*(fabs(jesup)+fabs(jesdown)));
  uncert.push_back(0.5*(fabs(scaledown)+fabs(scaleup)));
  uncert.push_back(resolution);
  
  uncertainties.push_back(uncert);
}

vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
  vector< vector<float> > uncertainties;
  for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
      cout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
      for(int ibin=0;ibin<bins.size();ibin++) {
        jzbSel=bins[ibin];
        geqleq="geq";
        do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,uncertainties,mcjzb,datajzb);
      }//end of bin loop
  }//end of signal loop
  return uncertainties;
}
Revision:	1.3
Committed:	Tue Jul 19 09:28:30 2011 UTC (13 years, 9 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+5 -5 lines
Log Message:	Changed jet number
#	User	Rev	Content
1	buchmann	1.1	#include <iostream>
2			#include <vector>
3			#include <sys/stat.h>
4
5			#include <TMath.h>
6			#include <TColor.h>
7			#include <TPaveText.h>
8			#include <TRandom.h>
9			#include <TF1.h>
10
11			#ifndef SampleClassLoaded
12			#include "ActiveSamples.C"
13			#endif
14
15			#ifndef Verbosity
16			#define Verbosity 0
17			#endif
18
19			#include <TFile.h>
20			#include <TTree.h>
21			#include <TH1.h>
22			#include <TCut.h>
23			#include <TMath.h>
24			#include <TLine.h>
25			#include <TCanvas.h>
26			#include <TProfile.h>
27			#include <TF1.h>
28
29
30
31			Int_t nBins = 100;
32			Float_t jzbMin = -207;
33			Float_t jzbMax = 243;
34			Float_t jzbSel = 100;
35			int iplot=0;
36			int verbose=0;
37	buchmann	1.2	string geqleq;
38			string mcjzbexpression;
39
40			TString geq_or_leq() {
41			if(geqleq=="geq") return TString(">=");
42			if(geqleq=="leq") return TString("<=");
43			return TString("GEQ_OR_LEQ_ERROR");
44			}
45	buchmann	1.1
46			//______________________________________________________________________________
47			Double_t Interpolate(Double_t x, TH1 *histo)
48			{
49			// Given a point x, approximates the value via linear interpolation
50			// based on the two nearest bin centers
51			// Andy Mastbaum 10/21/08
52			// in newer ROOT versions but not in the one I have so I had to work around that ...
53
54			Int_t xbin = histo->FindBin(x);
55			Double_t x0,x1,y0,y1;
56
57			if(x<=histo->GetBinCenter(1)) {
58			return histo->GetBinContent(1);
59			} else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
60			return histo->GetBinContent(histo->GetNbinsX());
61			} else {
62			if(x<=histo->GetBinCenter(xbin)) {
63			y0 = histo->GetBinContent(xbin-1);
64			x0 = histo->GetBinCenter(xbin-1);
65			y1 = histo->GetBinContent(xbin);
66			x1 = histo->GetBinCenter(xbin);
67			} else {
68			y0 = histo->GetBinContent(xbin);
69			x0 = histo->GetBinCenter(xbin);
70			y1 = histo->GetBinContent(xbin+1);
71			x1 = histo->GetBinCenter(xbin+1);
72			}
73			return y0 + (x-x0)*((y1-y0)/(x1-x0));
74			}
75			}
76
77
78			//____________________________________________________________________________________
79			// Efficiency plot
80			TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
81			iplot++;
82			int count=iplot;
83			// Define new histogram
84			char hname[30]; sprintf(hname,"hJzbEff%d",count);
85			TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
86			nBins,jzbMin,jzbMax);
87			Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
88
89	buchmann	1.2	events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
90	buchmann	1.1	Float_t maxEff = events->GetSelectedRows();
91			if(verbose>0) std::cout << hname << " (" << informalname <<") " << maxEff << std::endl;
92
93			if(verbose>0) std::cout << "JZB max = " << jzbMax << std::endl;
94			// Loop over steps to get efficiency curve
95			char cut[256];
96			for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
97			sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
98	buchmann	1.2	events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
99	buchmann	1.1	Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
100			// std::cout << "COUCOU " << __LINE__ << std::endl;
101			hJzbEff->SetBinContent(iBin+1,eff);
102			hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
103			}
104			return hJzbEff;
105
106
107			}
108
109
110			//________________________________________________________________________________________
111			// Pile-up efficiency
112			float pileup(TTree *events, string informalname, Float_t myJzbMax = 140. ) {
113			nBins = 16;
114			jzbMax = myJzbMax;
115
116			// Acceptance cuts
117	buchmann	1.3	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)(id2+1)ch1*ch2)!=-2");
118	buchmann	1.1
119			TH1F* hLM4 = plotEff(events,kbase,informalname);
120			hLM4->SetMinimum(0.);
121
122			// Nominal function
123			TF1* func = new TF1("func","0.5TMath::Erfc([0]x-[1])",jzbMin,jzbMax);
124			func->SetParameter(0,0.03);
125			func->SetParameter(1,0.);
126			hLM4->Fit(func,"Q");
127
128			// Pimped-up function
129			TF1* funcUp = (TF1*)func->Clone();
130			funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
131			std::cout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
132			<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
133
134			return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
135
136			}
137
138			//____________________________________________________________________________________
139			// Total selection efficiency (MC)
140	buchmann	1.2	void MCefficiency(TTree *events,float &res, float &reserr,string mcjzb) {
141	buchmann	1.1
142			char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
143			// All acceptance cuts at gen. level
144	buchmann	1.3	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
145	buchmann	1.1	// Corresponding reco. cuts
146	buchmann	1.2	TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
147	buchmann	1.1
148	buchmann	1.2	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
149	buchmann	1.1	Float_t sel = events->GetSelectedRows();
150	buchmann	1.2	events->Draw(mcjzbexpression.c_str(),kbase,"goff");
151	buchmann	1.1	Float_t tot = events->GetSelectedRows();
152
153			res=sel/tot;
154			reserr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
155			std::cout << " MC efficiency: " << res << "+-" << reserr << std::endl;
156			}
157
158			float JZBefficiency(TTree *events, string informalname) {
159	buchmann	1.3	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)(id2+1)ch1*ch2)!=-2");
160	buchmann	1.1	TH1F* hLM4 = plotEff(events,kbase,informalname);
161			Int_t bin = hLM4->FindBin(jzbSel); // To get the error
162			std::cout << " Efficiency at JZB==" << jzbSel << std::endl;
163			std::cout << " " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin) << std::endl;
164			return -1;
165			}
166
167			//________________________________________________________________________
168			// Effect of energy scale on efficiency
169			void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
170			TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0");
171			TCut ksel("abs(mll-91.2)<20&&((id1+1)(id2+1)ch1*ch2)!=-2");
172			TCut nJets("pfJetGoodNum>2");
173			stringstream down,up;
174			down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
175			up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
176
177			TCut nJetsP(up.str().c_str());
178			TCut nJetsM(down.str().c_str());
179
180			if ( jzbSelection>0 ) jzbSel = jzbSelection;
181
182			if ( !(plotName.Length()>1) ) plotName = informalname;
183
184			nBins = 1; jzbMin = jzbSel0.95; jzbMax = jzbSel1.05;
185			TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
186
187			TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
188
189			TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
190
191			// Dump some information
192			Float_t eff = Interpolate(jzbSel,hist);
193			Float_t effp = Interpolate(jzbSel,histp);
194			Float_t effm = Interpolate(jzbSel,histm);
195			std::cout << " Efficiency at JZB==" << jzbSel << std::endl;
196			std::cout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
197			std::cout << " central: " << eff << std::endl;
198			std::cout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
199			jesup=(effp-eff)/eff*100.;
200			jesdown=(effm-eff)/eff*100.;
201			}
202
203			//________________________________________________________________________
204			// Effect of energy scale on JZB efficiency
205			void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
206
207	buchmann	1.3	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&genNjets>2");
208	buchmann	1.1	TCut ksel("abs(mll-91.2)<20&&((id1+1)(id2+1)ch1*ch2)!=-2");
209
210			nBins = 50;
211			jzbMin = 0.5*jzbSel;
212			jzbMax = 2.0*jzbSel;
213
214			TH1F* hist = plotEff(events,kbase&&ksel,informalname);
215
216			// Dump some information
217			Float_t eff = Interpolate(jzbSel,hist);
218			Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
219			Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
220			std::cout << " efficiency at JZB==" << jzbSel(1.+systematic) << "(-"<<syst100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
221			std::cout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
222			std::cout << " efficiency at JZB==" << jzbSel(1.-systematic) << "(-"<<syst100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
223			up=((effp-eff)/eff)*100;
224			down=((effm-eff)/eff)*100;
225			}
226
227			//________________________________________________________________________
228			// JZB response (true/reco. vs. true)
229			void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
230
231			jzbMin = 20;
232	buchmann	1.3	TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&genNjets>2");
233	buchmann	1.1	TCut ksel("abs(mll-91.2)<20&&((id1+1)(id2+1)ch1*ch2)!=-2");
234
235			TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true",
236			nPeriods, jzbMin, myJzbMax, "" );
237
238	buchmann	1.2	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
239	buchmann	1.1	else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
240
241			hJzbResp->SetMaximum(1.2);
242			hJzbResp->SetMinimum(0.2);
243			hJzbResp->Fit("pol0","Q");
244			TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
245			cout << " Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
246			}
247
248
249	buchmann	1.2	void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &uncertainties,string mcjzb,string datajzb) {
250	buchmann	1.1
251			float JetEnergyScaleUncert=0.1;
252			float JZBScaleUncert=0.1;
253	buchmann	1.2	mcjzbexpression=mcjzb;
254	buchmann	1.1
255			float triggereff=4;//percent!
256			cout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
257			float leptonseleff=2;//percent!
258			cout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
259
260			float mceff,mcefferr;
261			cout << "MC efficiencies:" << endl;
262	buchmann	1.2	MCefficiency(events,mceff,mcefferr,mcjzb);
263	buchmann	1.1	JZBefficiency(events,informalname);
264
265			std::cout << "Jet energy scale: " << std::endl;
266			float jesup,jesdown;
267			JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
268
269			std::cout << "JZB scale: " << std::endl;
270			float scaleup,scaledown,scalesyst;
271			doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
272
273			std::cout << "JZB response: " << std::endl;
274			JZBresponse(events);
275
276			std::cout << "Pileup: " << std::endl;
277			float resolution=pileup(events,informalname);
278
279			cout << "_______________________________________________" << endl;
280			cout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
281			cout << "Trigger efficiency: " << triggereff << endl;
282			cout << "Lepton Sel Eff: " << leptonseleff << endl;
283			cout << "For JZB>" << jzbSel << endl;
284			cout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
285			cout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
286			cout << "Resolution : " << resolution << endl;
287
288			vector<float> uncert;
289			uncert.push_back(jzbSel);
290			uncert.push_back(triggereff);
291			uncert.push_back(leptonseleff);
292			uncert.push_back(0.5*(fabs(jesup)+fabs(jesdown)));
293			uncert.push_back(0.5*(fabs(scaledown)+fabs(scaleup)));
294			uncert.push_back(resolution);
295
296			uncertainties.push_back(uncert);
297			}
298
299			vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
300			vector< vector<float> > uncertainties;
301			for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
302			cout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
303			for(int ibin=0;ibin<bins.size();ibin++) {
304			jzbSel=bins[ibin];
305	buchmann	1.2	geqleq="geq";
306			do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,uncertainties,mcjzb,datajzb);
307	buchmann	1.1	}//end of bin loop
308			}//end of signal loop
309			return uncertainties;
310			}