ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
Revision: 1.55
Committed: Thu Nov 24 08:15:39 2011 UTC (13 years, 5 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.54: +1 -1 lines
Log Message: 
Only display peak error result if requested (i.e. automatized mode is off)

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4 buchmann 1.31 #include <algorithm>
5     #include <cmath>
6 buchmann 1.1
7     #include <TMath.h>
8     #include <TColor.h>
9     #include <TPaveText.h>
10     #include <TRandom.h>
11     #include <TF1.h>
12    
13     #ifndef SampleClassLoaded
14     #include "ActiveSamples.C"
15     #endif
16    
17     #ifndef Verbosity
18     #define Verbosity 0
19     #endif
20    
21     #include <TFile.h>
22     #include <TTree.h>
23     #include <TH1.h>
24     #include <TCut.h>
25     #include <TMath.h>
26     #include <TLine.h>
27     #include <TCanvas.h>
28     #include <TProfile.h>
29     #include <TF1.h>
30    
31    
32    
33     Int_t nBins = 100;
34     Float_t jzbMin = -207;
35     Float_t jzbMax = 243;
36     Float_t jzbSel = 100;
37     int iplot=0;
38     int verbose=0;
39 buchmann 1.2 string geqleq;
40     string mcjzbexpression;
41 buchmann 1.4 bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
42 buchmann 1.2
43     TString geq_or_leq() {
44     if(geqleq=="geq") return TString(">=");
45     if(geqleq=="leq") return TString("<=");
46     return TString("GEQ_OR_LEQ_ERROR");
47     }
48 buchmann 1.1
49 buchmann 1.6 TString ngeq_or_leq() {
50     if(geqleq=="geq") return TString("<=");
51     if(geqleq=="leq") return TString(">=");
52     return TString("NGEQ_OR_LEQ_ERROR");
53     }
54    
55 buchmann 1.1 //______________________________________________________________________________
56     Double_t Interpolate(Double_t x, TH1 *histo)
57     {
58     // Given a point x, approximates the value via linear interpolation
59     // based on the two nearest bin centers
60     // Andy Mastbaum 10/21/08
61     // in newer ROOT versions but not in the one I have so I had to work around that ...
62    
63     Int_t xbin = histo->FindBin(x);
64     Double_t x0,x1,y0,y1;
65    
66     if(x<=histo->GetBinCenter(1)) {
67     return histo->GetBinContent(1);
68     } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
69     return histo->GetBinContent(histo->GetNbinsX());
70     } else {
71     if(x<=histo->GetBinCenter(xbin)) {
72     y0 = histo->GetBinContent(xbin-1);
73     x0 = histo->GetBinCenter(xbin-1);
74     y1 = histo->GetBinContent(xbin);
75     x1 = histo->GetBinCenter(xbin);
76     } else {
77     y0 = histo->GetBinContent(xbin);
78     x0 = histo->GetBinCenter(xbin);
79     y1 = histo->GetBinContent(xbin+1);
80     x1 = histo->GetBinCenter(xbin+1);
81     }
82     return y0 + (x-x0)*((y1-y0)/(x1-x0));
83     }
84     }
85    
86 buchmann 1.7 //____________________________________________________________________________________
87     // Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
88 buchmann 1.45 float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut, int flipped) {
89 buchmann 1.7 iplot++;
90     int count=iplot;
91 buchmann 1.45 string locmcjzbexpression=mcjzbexpression;
92 buchmann 1.7 // Define new histogram
93     string hname=GetNumericHistoName();
94     TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
95 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
96 buchmann 1.7 hname=GetNumericHistoName();
97     TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
98 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
99 buchmann 1.7
100     hname=GetNumericHistoName();
101     TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
103 buchmann 1.7 hname=GetNumericHistoName();
104     TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
105 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
106 buchmann 1.7
107 buchmann 1.42 float obs=0;
108     float pred=0;
109     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
110     if(PlottingSetup::RestrictToMassPeak) {
111     hname=GetNumericHistoName();
112     TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
114 buchmann 1.42 hname=GetNumericHistoName();
115     TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
117 buchmann 1.42
118     hname=GetNumericHistoName();
119     TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
120 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
121 buchmann 1.42 hname=GetNumericHistoName();
122     TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
123 buchmann 1.45 events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
124 buchmann 1.42
125     obs = hossfp->Integral();
126     pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
127 buchmann 1.45
128     if(flipped>0) {
129     obs = hossfn->Integral();
130     pred= hossfp->Integral() - (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
131     }
132 buchmann 1.42 delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133     } else {
134     obs = hossfp->Integral();
135     pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
136 buchmann 1.45 if(flipped>0) {
137     obs = hossfn->Integral();
138     pred= hossfp->Integral() - (hosofp->Integral() - hosofn->Integral());;
139     }
140 buchmann 1.42 }
141 buchmann 1.7
142     delete hossfp,hossfn,hosofp,hosofn;
143     return obs-pred;
144     }
145    
146 buchmann 1.1
147     //____________________________________________________________________________________
148     // Efficiency plot
149 buchmann 1.45 TH1F* plotEff(TTree* events, TCut kbase, TString informalname, int flipped) {
150 buchmann 1.1 iplot++;
151     int count=iplot;
152 buchmann 1.50 iplot++;
153     int count2=iplot;
154 buchmann 1.1 // Define new histogram
155     char hname[30]; sprintf(hname,"hJzbEff%d",count);
156 buchmann 1.50 char hname2[30]; sprintf(hname2,"hJzbEff%d",count2);
157 fronga 1.53 TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB [GeV]; Efficiency",nBins,jzbMin,jzbMax);
158     TH1F* hJzbEff2= new TH1F(hname2,"JZB selection efficiency ; JZB [GeV]; Efficiency",1,-14000,14000);
159 buchmann 1.1 Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
160 buchmann 1.50
161     if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
162 buchmann 1.51 else events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
163 buchmann 1.50 Float_t maxEff = hJzbEff2->Integral();
164 buchmann 1.5 if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff << std::endl;
165 buchmann 1.1
166 buchmann 1.5 if(verbose>0) dout << "JZB max = " << jzbMax << std::endl;
167 buchmann 1.1 // Loop over steps to get efficiency curve
168     char cut[256];
169     for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
170 buchmann 1.6 sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
171 buchmann 1.50 if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
172 buchmann 1.51 if(flipped>0) events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
173 buchmann 1.50 Float_t eff = static_cast<Float_t>(hJzbEff2->Integral())/maxEff;
174 buchmann 1.1 hJzbEff->SetBinContent(iBin+1,eff);
175     hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
176     }
177 buchmann 1.50 delete hJzbEff2;
178 buchmann 1.1 return hJzbEff;
179     }
180    
181    
182 fronga 1.53
183 buchmann 1.1 //________________________________________________________________________________________
184 pablom 1.23 // Master Formula
185     void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
186    
187     float x0 = eff[0];
188     float deltaPos = 0, deltaNeg = 0;
189     for(int k = 0; k < (eff.size()-1)/2; k++) {
190     float xneg = eff[2*k+2];
191     float xpos = eff[2*k+1];
192     if(xpos-x0>0 || xneg-x0>0) {
193     if(xpos-x0 > xneg-x0) {
194     deltaPos += (xpos-x0)*(xpos-x0);
195     } else {
196     deltaPos += (xneg-x0)*(xneg-x0);
197     }
198     }
199     if(x0-xpos>0 || x0-xneg>0) {
200     if(x0-xpos > x0-xneg) {
201     deltaNeg += (xpos-x0)*(xpos-x0);
202     } else {
203     deltaNeg += (xneg-x0)*(xneg-x0);
204     }
205     }
206     }
207     errHi = sqrt(deltaPos);
208     errLo = sqrt(deltaNeg);
209    
210     }
211    
212    
213     //________________________________________________________________________________________
214     // Get normalization factor for the PDFs
215 buchmann 1.42 float get_norm_pdf_factor(TTree *events, int k, string addcut) {
216 pablom 1.23
217 buchmann 1.25 TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
218 pablom 1.23 char nameVar[20];
219     sprintf(nameVar, "pdfW[%d]", k);
220 buchmann 1.42 events->Project("haux", nameVar, addcut.c_str());
221 buchmann 1.24 float thisW = haux->Integral();
222     events->Project("haux", "pdfW[0]");
223     float normW = haux->Integral();
224    
225     float factor=thisW/normW;
226 pablom 1.23
227     delete haux;
228    
229     return factor;
230    
231     }
232    
233    
234    
235     //________________________________________________________________________________________
236 buchmann 1.1 // Pile-up efficiency
237 buchmann 1.45 float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
238 buchmann 1.1 nBins = 16;
239     jzbMax = myJzbMax;
240    
241     // Acceptance cuts
242 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
243     TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
244 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
245    
246 buchmann 1.42 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
247 buchmann 1.45 TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
248 buchmann 1.1 hLM4->SetMinimum(0.);
249    
250     // Nominal function
251     TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
252     func->SetParameter(0,0.03);
253     func->SetParameter(1,0.);
254     hLM4->Fit(func,"Q");
255    
256     // Pimped-up function
257     TF1* funcUp = (TF1*)func->Clone();
258 buchmann 1.47 funcUp->SetParameter( 0, func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
259 buchmann 1.5 if(!automatized) dout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
260 buchmann 1.1 << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
261    
262 buchmann 1.18 return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
263 buchmann 1.1
264     }
265    
266     //____________________________________________________________________________________
267 buchmann 1.6 // Effect of peak shifting
268 buchmann 1.45 void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
269 buchmann 1.42 //Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
270 buchmann 1.6 TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
271     TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
272     TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
273     TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
274     TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
275     TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
276     nBins = 1;
277     string informalname="PeakErrorCalculation";
278     float resup,resdown,rescent;
279     for(int i=0;i<3;i++) {
280     string poscut,negcut;
281     if(i==0) {
282     poscut=peakcentral;
283     negcut=npeakcentral;
284     } else if(i==1) {
285     poscut=peakdown;
286     negcut=npeakdown;
287     } else if(i==2) {
288     poscut=peakup;
289     negcut=npeakup;
290     }
291 buchmann 1.10 float res;
292 buchmann 1.45 if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
293     else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
294 buchmann 1.7 if(i==0) rescent=res;
295     else if(i==1) resdown=res;
296     else if(i==2) resup=res;
297 buchmann 1.6 }
298 buchmann 1.41 if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
299     else result=(TMath::Abs(rescent-resdown)/(float)rescent);
300 buchmann 1.55 if(!automatized) cout << " " << result << endl;
301 buchmann 1.6 }
302    
303 pablom 1.43
304 buchmann 1.45 void MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
305 pablom 1.43 if(!events) {
306     write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
307     result=0;resulterr=0;
308     return;
309     }
310    
311     char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
312 pablom 1.54 char metSel[256]; sprintf(metSel, "met[4] > %f", jzbSel);
313     string metSelection(metSel);
314 pablom 1.43 // All acceptance cuts at gen. level
315     //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
316     TCut kbase("");
317     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
318     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
319     // Corresponding reco. cuts
320    
321     TCut acceptance("genPt2 != 0");
322     TCut massId(cutmass&&cutOSSF);
323     TCut njets(cutnJets);
324 buchmann 1.45 TCut jzbp;
325     TCut jzbn;
326 pablom 1.54 TCut met(("pfJetGoodNum > 1 && abs(mll-91.2) < 10.0 && id1 == id2 &&" + metSelection).c_str());
327 buchmann 1.45 if(flipped==0) {
328     jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
329     jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
330     } else {
331     jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
332     jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
333     }
334 pablom 1.54 float ntotal = events->Draw("pt1", addcut.c_str(), "goff");
335 pablom 1.43 TCut theCut;
336     switch(type) {
337     case 1:
338     theCut = kbase+acceptance;
339     break;
340     case 2:
341     theCut = kbase+massId;
342     break;
343     case 3:
344     theCut = kbase+massId+njets;
345     break;
346     case 4:
347     theCut = kbase+massId+njets+jzbn;
348     break;
349 pablom 1.54 case 5:
350     theCut = kbase + met;
351     break;
352 pablom 1.43 default:
353     theCut = kbase+massId+njets+jzbn;
354     break;
355     }
356 pablom 1.54
357 pablom 1.43 string stheCut(theCut);
358     char var[20];
359     sprintf(var, "pdfW[%d]", k);
360    
361     string svar(var);
362     string newtheCut;
363     if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
364     else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
365    
366     TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
367     if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
368     else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
369     Float_t sel = effh->Integral();
370     Float_t nsel=0;
371     //Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
372     float normFactor = 1;
373     if(k>=0) get_norm_pdf_factor(events, k, addcut);
374     sel = sel/normFactor;
375    
376     result=(sel)/ntotal;
377     resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
378    
379 buchmann 1.52 delete effh;
380 pablom 1.43 }
381    
382 buchmann 1.6 //____________________________________________________________________________________
383 buchmann 1.1 // Total selection efficiency (MC)
384 buchmann 1.33 //returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
385 buchmann 1.45 Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
386 buchmann 1.40 if(!events) {
387 buchmann 1.42 write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
388     result=0;
389     resulterr=0;
390     return Value(0,0);
391 buchmann 1.40 }
392 buchmann 1.42
393 buchmann 1.1 char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
394     // All acceptance cuts at gen. level
395 buchmann 1.14 //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
396     TCut kbase("");
397 buchmann 1.42
398     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
399     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
400 buchmann 1.9 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
401 buchmann 1.1 // Corresponding reco. cuts
402 buchmann 1.42
403     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
404     TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
405     TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
406     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
407     if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
408     ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
409     ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
410 buchmann 1.45 if(flipped>0) {
411     ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
412     ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
413     }
414 buchmann 1.42 } else {
415     //for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
416     ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
417     ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
418 buchmann 1.45 if(flipped>0) {
419     ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
420     ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
421     }
422 buchmann 1.42 }
423    
424 pablom 1.23 TCut posSide = kbase&&ksel;
425     TCut negSide = kbase&&ksel2;
426     string sposSide(posSide);
427     string snegSide(negSide);
428     char var[20];
429     sprintf(var, "pdfW[%d]", k);
430 buchmann 1.42 if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
431 pablom 1.23 string svar(var);
432 buchmann 1.42 string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
433     string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
434     string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
435     string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
436 pablom 1.23
437 buchmann 1.24 TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
438 buchmann 1.50 if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), TCut(newPosSide.c_str())*PlottingSetup::Weight,"goff");
439     else events->Draw((mcjzbexpression+">>effh").c_str(), TCut((sposSide+"&&(id1==id2)").c_str())*PlottingSetup::Weight,"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
440 buchmann 1.52
441 buchmann 1.24 Float_t sel = effh->Integral();
442 buchmann 1.21 Float_t nsel=0;
443 buchmann 1.42
444     ///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
445    
446 buchmann 1.21 if(ConsiderSignalContaminationForLimits) {
447 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
448     if(PlottingSetup::RestrictToMassPeak) {
449 buchmann 1.50 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
450 buchmann 1.42 nsel += effh->Integral();
451     } else {
452 buchmann 1.50 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
453 buchmann 1.42 nsel += effh->Integral();
454 buchmann 1.50 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewPosSide.c_str()))*PlottingSetup::Weight,"goff");
455 buchmann 1.42 nsel += effh->Integral();
456 buchmann 1.50 events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewNegSide.c_str()))*PlottingSetup::Weight,"goff");
457 buchmann 1.42 nsel -= effh->Integral();
458     }
459 buchmann 1.21 }
460 buchmann 1.42
461 pablom 1.23 //Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
462 buchmann 1.26 float normFactor = 1;
463 buchmann 1.42 if(k>=0) get_norm_pdf_factor(events, k, addcut);
464 pablom 1.23 sel = sel/normFactor;
465     nsel = nsel/normFactor;
466    
467 buchmann 1.13 Float_t tot = Neventsinfile;
468 buchmann 1.1
469 buchmann 1.33 Value result_wo_signalcont;
470    
471 buchmann 1.21 if(ConsiderSignalContaminationForLimits) {
472     result=(sel-nsel)/tot;
473     resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
474 buchmann 1.33 result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
475 buchmann 1.21 } else {//no signal contamination considered:
476     result=(sel)/tot;
477     resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
478 buchmann 1.33 result_wo_signalcont=Value(result,resulterr);
479 buchmann 1.21 }
480 pablom 1.43 if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
481 buchmann 1.42 if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
482 buchmann 1.24 delete effh;
483 buchmann 1.34 return result_wo_signalcont;
484 buchmann 1.1 }
485    
486 buchmann 1.29
487 buchmann 1.42
488 buchmann 1.29 //____________________________________________________________________________________
489 buchmann 1.30 // Selection efficiency for one process (MC)
490 buchmann 1.45 // not in use anymore.
491     /*
492     vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
493 buchmann 1.29 vector<float> process_efficiencies;
494     for(int iprocess=0;iprocess<=10;iprocess++) {
495     float this_process_efficiency,efferr;
496     stringstream addcutplus;
497     addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
498 buchmann 1.45 MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
499 buchmann 1.29 process_efficiencies.push_back(this_process_efficiency);
500     }
501     return process_efficiencies;
502     }
503 buchmann 1.45 */
504 buchmann 1.29
505 buchmann 1.45 void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
506 buchmann 1.42 TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
507 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
508 buchmann 1.42 if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
509 buchmann 1.45 TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
510 buchmann 1.1 Int_t bin = hLM4->FindBin(jzbSel); // To get the error
511 buchmann 1.8 jzbeff=Interpolate(jzbSel,hLM4);
512     jzbefferr=hLM4->GetBinError(bin);
513 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
514 buchmann 1.8 if(!automatized) dout << " " << jzbeff << "+-" << jzbefferr << std::endl;
515 buchmann 1.1 }
516    
517     //________________________________________________________________________
518     // Effect of energy scale on efficiency
519 buchmann 1.48 void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",Float_t jzbSelection=-1, TString plotName = "" ) {
520 buchmann 1.42 TCut kbase(genMassCut&&"genZPt>0");
521 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
522 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
523     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
524 buchmann 1.8
525 buchmann 1.42 TCut ksel(cutmass&&cutOSSF);
526 buchmann 1.1 TCut nJets("pfJetGoodNum>2");
527     stringstream down,up;
528 buchmann 1.48 down << "pfJetGoodNumn1sigma>=3";
529     up << "pfJetGoodNump1sigma>=3";
530 buchmann 1.1
531     TCut nJetsP(up.str().c_str());
532     TCut nJetsM(down.str().c_str());
533    
534     if ( !(plotName.Length()>1) ) plotName = informalname;
535    
536     nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
537 buchmann 1.51 TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname,flipped);
538 buchmann 1.45 TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname,flipped);
539     TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname,flipped);
540 buchmann 1.1
541     // Dump some information
542     Float_t eff = Interpolate(jzbSel,hist);
543     Float_t effp = Interpolate(jzbSel,histp);
544     Float_t effm = Interpolate(jzbSel,histm);
545 buchmann 1.5 if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
546     if(!automatized) dout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
547     if(!automatized) dout << " central: " << eff << std::endl;
548     if(!automatized) dout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
549 buchmann 1.18 jesup=(effp-eff)/eff;
550     jesdown=(effm-eff)/eff;
551 buchmann 1.1 }
552    
553     //________________________________________________________________________
554     // Effect of energy scale on JZB efficiency
555 buchmann 1.45 void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
556 buchmann 1.1
557 buchmann 1.42 TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
558 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
559 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
560     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
561     TCut ksel(cutmass&&cutOSSF);
562 buchmann 1.1
563     nBins = 50;
564     jzbMin = 0.5*jzbSel;
565     jzbMax = 2.0*jzbSel;
566    
567 buchmann 1.45 TH1F* hist = plotEff(events,kbase&&ksel,informalname,flipped);
568 buchmann 1.1
569     // Dump some information
570     Float_t eff = Interpolate(jzbSel,hist);
571     Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
572     Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
573 buchmann 1.33 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
574 buchmann 1.5 if(!automatized) dout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
575 buchmann 1.33 if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
576 buchmann 1.18 up=((effp-eff)/eff);
577     down=((effm-eff)/eff);
578 buchmann 1.1 }
579    
580     //________________________________________________________________________
581     // JZB response (true/reco. vs. true)
582 fronga 1.53 void JZBresponse(TTree *events, string name, bool requireZ, float &resp, float &resperr, int flipped, string addcut="", bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
583 buchmann 1.1
584     jzbMin = 20;
585 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
586     TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
587 buchmann 1.10 if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
588 buchmann 1.42 flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
589     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
590     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
591     TCut ksel(cutmass&&cutOSSF);
592 buchmann 1.1
593 buchmann 1.8 TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
594 buchmann 1.1
595 buchmann 1.45 string locmcjzbexpression=mcjzbexpression;
596     if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
597     string possibleminus="";
598     if(flipped>0) possibleminus="-";
599     if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
600 buchmann 1.1 else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
601    
602     hJzbResp->SetMaximum(1.2);
603     hJzbResp->SetMinimum(0.2);
604     hJzbResp->Fit("pol0","Q");
605     TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
606 buchmann 1.32 if(!fittedfunction) {
607     // in case there are not enough points passing our selection
608     cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
609     resp=1;
610     resperr=1;
611     } else {
612     resp=fittedfunction->GetParameter(0);
613     resperr=fittedfunction->GetParError(0);
614     if(!automatized) dout << " Response: " << resp << " +/- " << resperr << endl;
615     }
616 buchmann 1.5 delete hJzbResp;
617 buchmann 1.1 }
618    
619    
620 pablom 1.23 //________________________________________________________________________________________
621     // PDF uncertainty
622 buchmann 1.45 float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
623 pablom 1.23 std::vector<float> efficiency;
624 buchmann 1.24 for(int k = 1; k < NPdfs; k++) {
625 pablom 1.23 float result, resulterr;
626 buchmann 1.45 Value flipval;
627     MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, addcut, k);
628 pablom 1.23 efficiency.push_back(result);
629     }
630 buchmann 1.24 float errHi, errLow,err;
631 pablom 1.23 master_formula(efficiency, errHi, errLow);
632 buchmann 1.24 err=errLow;
633     if(errHi>errLow) err=errHi;
634     if(!automatized) dout << " Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
635     return err;
636 pablom 1.23
637     }
638    
639 buchmann 1.24 int get_npdfs(TTree *events) {
640     int NPDFs;
641     events->SetBranchAddress("NPdfs",&NPDFs);
642     events->GetEntry(1);
643     return NPDFs;
644     }
645    
646 pablom 1.23
647 buchmann 1.45 void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,int flipped, string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
648 buchmann 1.48 float JZBScaleUncert=0.05;
649 buchmann 1.2 mcjzbexpression=mcjzb;
650 buchmann 1.48 float triggereff=2.0/100;// in range [0,1]
651 buchmann 1.5 dout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
652 buchmann 1.18 float leptonseleff=2.0/100;// in range [0,1]
653 buchmann 1.32 leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
654 buchmann 1.5 dout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
655 buchmann 1.1
656 buchmann 1.27 int NPdfs=0;
657     if(ismSUGRA) NPdfs = get_npdfs(events);
658 buchmann 1.29
659 buchmann 1.8 float mceff,mcefferr,jzbeff,jzbefferr;
660 buchmann 1.5 if(!automatized) dout << "MC efficiencies:" << endl;
661 buchmann 1.45 Value flipefficiency;
662     Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,addcut,-1);
663 buchmann 1.35 if(!automatized) cout << " Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
664    
665 buchmann 1.45 if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
666 buchmann 1.8 if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
667 buchmann 1.1
668 buchmann 1.51 if(!automatized) dout << "Error from Peak position:";
669 buchmann 1.6 float sysfrompeak=0;
670 buchmann 1.45 PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
671 buchmann 1.6
672 buchmann 1.51 if(!automatized) dout << "Jet energy scale (JES): " << std::endl;
673 buchmann 1.42 float jesup,jesdown;
674 buchmann 1.48 JZBjetScale(events,jesdown,jesup,informalname,flipped,requireZ,addcut);
675 fronga 1.53
676 buchmann 1.5 if(!automatized) dout << "JZB scale: " << std::endl;
677 buchmann 1.42 float scaleup,scaledown,scalesyst;
678 buchmann 1.45 doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
679 buchmann 1.1
680 buchmann 1.5 if(!automatized) dout << "JZB response: " << std::endl;
681 buchmann 1.42 float resp,resperr;
682     if(PlottingSetup::computeJZBresponse) {
683 buchmann 1.32 if(!automatized) dout << "JZB response: " << std::endl;
684 fronga 1.53 if(!ismSUGRA) JZBresponse(events,informalname,requireZ,resp,resperr,flipped,addcut);
685 buchmann 1.32 }
686 buchmann 1.1
687 buchmann 1.5 if(!automatized) dout << "Pileup: " << std::endl;
688 buchmann 1.48 // float resolution;
689     //resolution=pileup(events,requireZ,informalname,flipped,addcut);
690 pablom 1.23
691 buchmann 1.22 float PDFuncert=0;
692 buchmann 1.42 if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
693 buchmann 1.45 if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
694 pablom 1.23
695 buchmann 1.5 dout << "_______________________________________________" << endl;
696 buchmann 1.10 dout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << " (all in %) ";
697     if(addcut!="") dout << "With additional cut: " << addcut;
698     dout << endl;
699 buchmann 1.18 dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
700     dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
701     dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
702     dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
703     dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
704 buchmann 1.48 // dout << "Resolution : " << resolution << endl; // in range [0,1]
705 buchmann 1.18 dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
706 buchmann 1.24 if(ismSUGRA) dout << "PDF uncertainty : " << PDFuncert << endl; // in range [0,1]
707 buchmann 1.32 if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
708     if(PlottingSetup::computeJZBresponse)dout << "JZB response : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
709 buchmann 1.1
710 buchmann 1.4 float toterr=0;
711 buchmann 1.18 toterr+=(triggereff)*(triggereff);
712     toterr+=(leptonseleff)*(leptonseleff);
713     if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
714     if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
715 buchmann 1.48 // toterr+=(resolution*resolution);
716 buchmann 1.18 toterr+=(sysfrompeak*sysfrompeak);
717 buchmann 1.22 if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
718 buchmann 1.18 dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
719 buchmann 1.19 float systerr=TMath::Sqrt(toterr)*mceff;
720     toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
721    
722     dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst) %" << endl;
723     dout << " we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
724 buchmann 1.36 dout << "_______________________________________________" << endl;
725 buchmann 1.19
726     //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
727 buchmann 1.4 vector<float> res;
728     res.push_back(jzbSel);
729     res.push_back(mceff);
730     res.push_back(mcefferr);
731     res.push_back(toterr);
732     res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
733 buchmann 1.18 if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
734     if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
735 buchmann 1.48 // res.push_back(fabs(resolution));
736 buchmann 1.49 res.push_back(0.0);
737 buchmann 1.35 res.push_back(mceff_nosigcont.getValue());
738     res.push_back(mceff_nosigcont.getError());
739     if(ismSUGRA) res.push_back(PDFuncert);
740 buchmann 1.34 results.push_back(res);
741 buchmann 1.1 }
742    
743 buchmann 1.45 vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
744 buchmann 1.4 automatized=true;
745     vector< vector<float> > systematics;
746 buchmann 1.1 for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
747 fronga 1.53 dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
748     for(int ibin=0;ibin<bins.size();ibin++) {
749     jzbSel=bins[ibin];
750     geqleq="geq";
751     do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,mcjzb,datajzb,mcpeakerror,requireZ);
752     }//end of bin loop
753 buchmann 1.1 }//end of signal loop
754 buchmann 1.4 return systematics;
755 buchmann 1.1 }