ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
Revision: 1.5
Committed: Wed Jul 20 08:51:33 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.4: +35 -34 lines
Log Message:
Adapted the output; all output is now written on screen and to a file simultaneously

File Contents

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