ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/Systematics.C
Revision: 1.7
Committed: Mon May 28 10:24:25 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.6: +7 -4 lines
Log Message: 
Correcting for (possible) reweighting bias

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3 buchmann 1.3 #include <assert.h>
4 buchmann 1.1 #include <sys/stat.h>
5     #include <algorithm>
6     #include <cmath>
7    
8     #include <TMath.h>
9     #include <TColor.h>
10     #include <TPaveText.h>
11     #include <TRandom.h>
12     #include <TF1.h>
13    
14     #ifndef SampleClassLoaded
15     #include "SampleClass.C"
16     #endif
17    
18 buchmann 1.3 /*#ifndef CrossSectionReaderLoaded
19     #include "CrossSectionReader.C"
20     #endif*/
21    
22 buchmann 1.1 #ifndef Verbosity
23     #define Verbosity 0
24     #endif
25    
26     #include <TFile.h>
27     #include <TTree.h>
28     #include <TH1.h>
29     #include <TCut.h>
30     #include <TMath.h>
31     #include <TLine.h>
32     #include <TCanvas.h>
33     #include <TProfile.h>
34     #include <TF1.h>
35    
36    
37    
38     Int_t nBins = 100;
39     Float_t jzbMin = -207;
40     Float_t jzbMax = 243;
41     Float_t jzbSel = 100;
42     int iplot=0;
43     int verbose=0;
44     string geqleq;
45     string mcjzbexpression;
46     bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
47    
48     TString geq_or_leq() {
49     if(geqleq=="geq") return TString(">=");
50     if(geqleq=="leq") return TString("<=");
51     return TString("GEQ_OR_LEQ_ERROR");
52     }
53    
54     TString ngeq_or_leq() {
55     if(geqleq=="geq") return TString("<=");
56     if(geqleq=="leq") return TString(">=");
57     return TString("NGEQ_OR_LEQ_ERROR");
58     }
59    
60     //______________________________________________________________________________
61     Double_t Interpolate(Double_t x, TH1 *histo)
62     {
63     // Given a point x, approximates the value via linear interpolation
64     // based on the two nearest bin centers
65     // Andy Mastbaum 10/21/08
66     // in newer ROOT versions but not in the one I have so I had to work around that ...
67    
68     Int_t xbin = histo->FindBin(x);
69     Double_t x0,x1,y0,y1;
70    
71     if(x<=histo->GetBinCenter(1)) {
72     return histo->GetBinContent(1);
73     } else if(x>=histo->GetBinCenter(histo->GetNbinsX())) {
74     return histo->GetBinContent(histo->GetNbinsX());
75     } else {
76     if(x<=histo->GetBinCenter(xbin)) {
77     y0 = histo->GetBinContent(xbin-1);
78     x0 = histo->GetBinCenter(xbin-1);
79     y1 = histo->GetBinContent(xbin);
80     x1 = histo->GetBinCenter(xbin);
81     } else {
82     y0 = histo->GetBinContent(xbin);
83     x0 = histo->GetBinCenter(xbin);
84     y1 = histo->GetBinContent(xbin+1);
85     x1 = histo->GetBinCenter(xbin+1);
86     }
87     return y0 + (x-x0)*((y1-y0)/(x1-x0));
88     }
89     }
90    
91     //____________________________________________________________________________________
92     // Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
93     float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut, int flipped) {
94     iplot++;
95     string locmcjzbexpression=mcjzbexpression;
96     // Define new histogram
97     string hname=GetNumericHistoName();
98     TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
99     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
100     hname=GetNumericHistoName();
101     TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
103    
104     hname=GetNumericHistoName();
105     TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
106     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
107     hname=GetNumericHistoName();
108     TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
109     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
110    
111     float obs=0;
112     float pred=0;
113     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
114     if(PlottingSetup::RestrictToMassPeak) {
115     hname=GetNumericHistoName();
116     TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
117     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
118     hname=GetNumericHistoName();
119     TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
120     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
121    
122     hname=GetNumericHistoName();
123     TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
124     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
125     hname=GetNumericHistoName();
126     TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
127     events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
128    
129     obs = hossfp->Integral();
130     pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
131    
132     if(flipped>0) {
133     obs = hossfn->Integral();
134     pred= hossfp->Integral() - (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
135     }
136 buchmann 1.6 delete sbhossfp;
137     delete sbhossfn;
138     delete sbhosofp;
139     delete sbhosofn;
140 buchmann 1.1 } else {
141     obs = hossfp->Integral();
142     pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
143     if(flipped>0) {
144     obs = hossfn->Integral();
145     pred= hossfp->Integral() - (hosofp->Integral() - hosofn->Integral());;
146     }
147     }
148    
149 buchmann 1.6 delete hossfp;
150     delete hossfn;
151     delete hosofp;
152     delete hosofn;
153    
154 buchmann 1.1 return obs-pred;
155     }
156    
157    
158     //____________________________________________________________________________________
159     // Efficiency plot
160     TH1F* plotEff(TTree* events, TCut kbase, TString informalname, int flipped) {
161     iplot++;
162     int count=iplot;
163     iplot++;
164     int count2=iplot;
165     // Define new histogram
166     char hname[30]; sprintf(hname,"hJzbEff%d",count);
167     char hname2[30]; sprintf(hname2,"hJzbEff%d",count2);
168     TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB [GeV]; Efficiency",nBins,jzbMin,jzbMax);
169     TH1F* hJzbEff2= new TH1F(hname2,"JZB selection efficiency ; JZB [GeV]; Efficiency",1,-14000,14000);
170     Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
171    
172     if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
173     else events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
174     Float_t maxEff = hJzbEff2->Integral();
175     if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff << std::endl;
176    
177     if(verbose>0) dout << "JZB max = " << jzbMax << std::endl;
178     // Loop over steps to get efficiency curve
179     char cut[256];
180     for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
181     sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
182     if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
183     if(flipped>0) events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
184     Float_t eff = static_cast<Float_t>(hJzbEff2->Integral())/maxEff;
185     hJzbEff->SetBinContent(iBin+1,eff);
186     hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
187     }
188     delete hJzbEff2;
189     return hJzbEff;
190     }
191    
192    
193    
194     //________________________________________________________________________________________
195     // Master Formula
196     void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
197    
198     float x0 = eff[0];
199     float deltaPos = 0, deltaNeg = 0;
200 buchmann 1.6 for(int k = 0; k < ((int)eff.size()-1)/2; k++) {
201 buchmann 1.1 float xneg = eff[2*k+2];
202     float xpos = eff[2*k+1];
203     if(xpos-x0>0 || xneg-x0>0) {
204     if(xpos-x0 > xneg-x0) {
205     deltaPos += (xpos-x0)*(xpos-x0);
206     } else {
207     deltaPos += (xneg-x0)*(xneg-x0);
208     }
209     }
210     if(x0-xpos>0 || x0-xneg>0) {
211     if(x0-xpos > x0-xneg) {
212     deltaNeg += (xpos-x0)*(xpos-x0);
213     } else {
214     deltaNeg += (xneg-x0)*(xneg-x0);
215     }
216     }
217     }
218     errHi = sqrt(deltaPos);
219     errLo = sqrt(deltaNeg);
220    
221     }
222    
223    
224     //________________________________________________________________________________________
225     // Get normalization factor for the PDFs
226     float get_norm_pdf_factor(TTree *events, int k, string addcut) {
227    
228     TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
229     char nameVar[20];
230     sprintf(nameVar, "pdfW[%d]", k);
231     events->Project("haux", nameVar, addcut.c_str());
232     float thisW = haux->Integral();
233     events->Project("haux", "pdfW[0]");
234     float normW = haux->Integral();
235    
236     float factor=thisW/normW;
237    
238     delete haux;
239    
240     return factor;
241    
242     }
243    
244    
245    
246     //________________________________________________________________________________________
247     // Pile-up efficiency
248     float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
249     nBins = 16;
250     jzbMax = myJzbMax;
251    
252     // Acceptance cuts
253     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
254     TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
255     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
256    
257     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
258     TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
259     hLM4->SetMinimum(0.);
260    
261     // Nominal function
262     TF1* func = new TF1("func","0.5*TMath::Erfc([0]*x-[1])",jzbMin,jzbMax);
263     func->SetParameter(0,0.03);
264     func->SetParameter(1,0.);
265     hLM4->Fit(func,"Q");
266    
267     // Pimped-up function
268     TF1* funcUp = (TF1*)func->Clone();
269     funcUp->SetParameter( 0, func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
270     if(!automatized) dout << " PU: " << funcUp->Eval(jzbSel) << " " << func->Eval(jzbSel)
271     << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
272    
273     return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
274    
275     }
276    
277     //____________________________________________________________________________________
278     // Effect of peak shifting
279     void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
280     //Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
281     TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
282     TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
283     TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
284     TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
285     TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
286     TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
287     nBins = 1;
288     string informalname="PeakErrorCalculation";
289     float resup,resdown,rescent;
290     for(int i=0;i<3;i++) {
291     string poscut,negcut;
292     if(i==0) {
293     poscut=peakcentral;
294     negcut=npeakcentral;
295     } else if(i==1) {
296     poscut=peakdown;
297     negcut=npeakdown;
298     } else if(i==2) {
299     poscut=peakup;
300     negcut=npeakup;
301     }
302     float res;
303     if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
304     else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
305     if(i==0) rescent=res;
306     else if(i==1) resdown=res;
307     else if(i==2) resup=res;
308     }
309     if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
310     else result=(TMath::Abs(rescent-resdown)/(float)rescent);
311     if(!automatized) cout << " " << result << endl;
312 buchmann 1.4 // if(rescent==0) result=0;
313 buchmann 1.1 }
314    
315    
316     void MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
317     if(!events) {
318     write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
319     result=0;resulterr=0;
320     return;
321     }
322    
323     char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
324     char metSel[256]; sprintf(metSel, "met[4] > %f", jzbSel);
325     string metSelection(metSel);
326     // All acceptance cuts at gen. level
327     //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
328     TCut kbase("");
329     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
330     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
331     // Corresponding reco. cuts
332    
333     TCut acceptance("genPt2 != 0");
334     TCut massId(cutmass&&cutOSSF);
335     TCut njets(cutnJets);
336     TCut jzbp;
337     TCut jzbn;
338     TCut met(("pfJetGoodNum > 1 && abs(mll-91.2) < 10.0 && id1 == id2 &&" + metSelection).c_str());
339     if(flipped==0) {
340     jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
341     jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
342     } else {
343     jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
344     jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
345     }
346     float ntotal = events->Draw("pt1", addcut.c_str(), "goff");
347     TCut theCut;
348     switch(type) {
349     case 1:
350     theCut = kbase+acceptance;
351     break;
352     case 2:
353     theCut = kbase+massId;
354     break;
355     case 3:
356     theCut = kbase+massId+njets;
357     break;
358     case 4:
359     theCut = kbase+massId+njets+jzbn;
360     break;
361     case 5:
362     theCut = kbase + met;
363     break;
364     default:
365     theCut = kbase+massId+njets+jzbn;
366     break;
367     }
368    
369     string stheCut(theCut);
370     char var[20];
371     sprintf(var, "pdfW[%d]", k);
372    
373     string svar(var);
374     string newtheCut;
375     if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
376     else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
377    
378     TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
379     if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
380     else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
381     Float_t sel = effh->Integral();
382     //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.
383     float normFactor = 1;
384     if(k>=0) get_norm_pdf_factor(events, k, addcut);
385     sel = sel/normFactor;
386    
387     result=(sel)/ntotal;
388     resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
389    
390     delete effh;
391     }
392    
393 buchmann 1.3
394    
395     float XSForProcessViaAddCutWrapper(string addcut, map < pair<float, float>, map<string, float> > xsec, int i) {
396     int position = addcut.find("Abs(M0-");
397     string M0string=addcut.substr(position+7,4);
398     position=M0string.find(")");
399     if(position>0&&position<5) M0string=M0string.substr(0,position);
400     position = addcut.find("Abs(M12-");
401     string M12string=addcut.substr(position+8,4);
402     position=M0string.find(")");
403     if(position>0&&position<5) M12string=M12string.substr(0,position);
404     float m0=atof(M0string.c_str());
405     float m12=atof(M12string.c_str());
406     return GetXSecForPointAndChannel(m0,m12,xsec,i);
407     }
408    
409     float sum(vector<float> v) {
410     float sum=0;
411 buchmann 1.6 for(int i=0;i<(int)v.size();i++) sum+=v[i];
412 buchmann 1.3 return sum;
413     }
414     Value mSUGRAefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, int scantype, map < pair<float, float>, map<string, float> > xsec, string addcut="", int kwrong = -2) {
415     if(kwrong>0) {
416     write_error(__FUNCTION__,"Watch out, evaluation of PDF uncerts is done differently now .... asserting this and exiting, so long!");
417     assert(kwrong<=0);
418     }
419    
420     if(!events) {
421     write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
422     result=0;
423     resulterr=0;
424     return Value(0,0);
425     }
426     char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
427     // All acceptance cuts at gen. level
428     //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
429     TCut kbase(basiccut&&leptoncut);
430    
431     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
432     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
433     TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
434    
435     TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
436     if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
437     ksel=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
438     ksel2=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
439     if(flipped>0) {
440     ksel=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
441     ksel2=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
442     }
443     }else {
444     //for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
445     ksel=TCut(cutnJets&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
446     ksel2=TCut(cutnJets&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
447     if(flipped>0) {
448     ksel=TCut(cutnJets&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
449     ksel2=TCut(cutnJets&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
450     }
451     }
452    
453     TCut posSide = kbase&&ksel;
454     TCut negSide = kbase&&ksel2;
455     string sposSide(posSide);
456     string snegSide(negSide);
457     string newPosSide = "((id1==id2)&&(" + sposSide + "))";
458     string newNegSide = "((id1==id2)&&(" + snegSide + "))";
459     string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))"; // only used for off peak analysis
460     string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))"; // only used for off peak analysis
461    
462     TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
463     vector<float> sel;
464     vector<float> nsel;
465     vector<float> Nproc;
466    
467     for(int i=0;i<11;i++) {
468     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(((string)sposSide+"&&(id1==id2)&&process=="+any2string(i)).c_str()))*cutWeight,"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
469     sel.push_back(effh->Integral());
470     events->Draw(("id1>>effh"), (addcut+"&&process=="+any2string(i)).c_str(),"goff");
471     Nproc.push_back(effh->Integral());
472     }
473    
474    
475     ///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
476    
477     if(ConsiderSignalContaminationForLimits) {
478     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
479     if(PlottingSetup::RestrictToMassPeak) {
480     for(int i=0;i<11;i++) {
481     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut((newNegSide+"&&process=="+any2string(i)).c_str()))*cutWeight,"goff");
482     nsel.push_back(effh->Integral());
483     }
484     } else {
485     for(int i=0;i<11;i++) {
486     float nselproc=0;
487     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut((newNegSide+"&&process=="+any2string(i)).c_str()))*cutWeight,"goff");
488     nselproc += effh->Integral();
489     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut((emnewPosSide+"&&process=="+any2string(i)).c_str()))*cutWeight,"goff");
490     nselproc += effh->Integral();
491     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut((emnewNegSide+"&&process=="+any2string(i)).c_str()))*cutWeight,"goff");
492     nselproc -= effh->Integral();
493     nsel.push_back(nselproc);
494     }
495     }
496     }
497    
498     Value result_wo_signalcont;
499    
500     float resultwosigcont;
501     float resultwosigconterr;
502     if(ConsiderSignalContaminationForLimits) {
503     result=0;
504     resulterr=0;
505     float totalXS=0;
506     resultwosigcont=0;
507     resultwosigconterr=0;
508     for(int i=0;i<11;i++) {
509     float xsi=XSForProcessViaAddCutWrapper(addcut,xsec,i);
510     if(Nproc[i]<1) continue;
511     result+=((sel[i]-nsel[i])/Nproc[i])*xsi;
512     totalXS+=xsi;
513     resulterr+=xsi*(sel[i]+nsel[i]+(sel[i]-nsel[i])*(sel[i]-nsel[i])/Nproc[i])/(Nproc[i]*Nproc[i]);
514     resultwosigcont+=(sel[i]/Nproc[i])*xsi;
515     resultwosigconterr+=xsi*(sel[i]+(sel[i]*sel[i])/Nproc[i])/(Nproc[i]*Nproc[i]);
516     }
517     result=result/totalXS;
518     resulterr=TMath::Sqrt((1/totalXS)*resulterr);
519     resultwosigcont=resultwosigcont/totalXS;
520     resultwosigconterr=TMath::Sqrt((1/totalXS)*resultwosigconterr);
521     result_wo_signalcont=Value(resultwosigcont,resultwosigconterr);
522     } else {//no signal contamination considered:
523     result=0;
524     resulterr=0;
525     float totalXS=0;
526     for(int i=0;i<11;i++) {
527     float xsi=XSForProcessViaAddCutWrapper(addcut,xsec,i);
528     result+=((sel[i])/Nproc[i])*xsi;
529     totalXS+=xsi;
530     resulterr+=xsi*(sel[i]+(sel[i]*sel[i])/Nproc[i])/(Nproc[i]*Nproc[i]);
531     }
532     result=result/totalXS;
533     resulterr=TMath::Sqrt((1/totalXS)*resulterr);
534     result_wo_signalcont=Value(result,resulterr);
535     }
536    
537    
538     if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sum(sel) << " , signal contamination : " << sum(nsel) << " and nevents=" << sum(Nproc) << ") " << std::endl;
539     delete effh;
540     return result_wo_signalcont;
541    
542     }
543    
544    
545 buchmann 1.1 //____________________________________________________________________________________
546     // Total selection efficiency (MC)
547     //returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
548 buchmann 1.3 Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, int scantype, map < pair<float, float>, map<string, float> > xsec, string addcut="", int k = 0) {
549 buchmann 1.4 if(scantype==mSUGRA) {
550     write_warning(__FUNCTION__,"Need to check the efficiency algorithm for mSUGRA!");
551     //return mSUGRAefficiency(events,result,resulterr,flipped,mcjzb,requireZ,Neventsinfile,scantype,xsec,addcut,k);
552     }
553 buchmann 1.1 if(!events) {
554     write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
555     result=0;
556     resulterr=0;
557     return Value(0,0);
558     }
559    
560     char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
561     // All acceptance cuts at gen. level
562     //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
563 buchmann 1.2 TCut kbase(basiccut&&leptoncut);
564 buchmann 1.1
565     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
566     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
567     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
568     // Corresponding reco. cuts
569     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
570     TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
571     TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
572     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
573     if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
574 buchmann 1.2 ksel=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
575     ksel2=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
576 buchmann 1.1 if(flipped>0) {
577 buchmann 1.2 ksel=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
578     ksel2=TCut(cutnJets&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
579 buchmann 1.1 }
580     } else {
581     //for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
582 buchmann 1.2 ksel=TCut(cutnJets&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
583     ksel2=TCut(cutnJets&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
584 buchmann 1.1 if(flipped>0) {
585 buchmann 1.2 ksel=TCut(cutnJets&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
586     ksel2=TCut(cutnJets&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
587 buchmann 1.1 }
588     }
589    
590     TCut posSide = kbase&&ksel;
591     TCut negSide = kbase&&ksel2;
592     string sposSide(posSide);
593     string snegSide(negSide);
594     char var[20];
595     sprintf(var, "pdfW[%d]", k);
596     if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
597     string svar(var);
598     string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
599     string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
600     string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
601     string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
602    
603     TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
604 buchmann 1.7
605     events->Draw(("1>>effh"), cutWeight,"goff");
606     float averageweight=effh->Integral()/effh->GetEntries();
607     if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), TCut(newPosSide.c_str())*cutWeight*((string)"(1.0/"+any2string(averageweight)+")").c_str(),"goff");
608     else events->Draw((mcjzbexpression+">>effh").c_str(), TCut((sposSide+"&&(id1==id2)").c_str())*cutWeight*((string)"(1.0/"+any2string(averageweight)+")").c_str(),"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
609 buchmann 1.1 Float_t sel = effh->Integral();
610     Float_t nsel=0;
611    
612     ///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
613    
614     if(ConsiderSignalContaminationForLimits) {
615     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
616     if(PlottingSetup::RestrictToMassPeak) {
617     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*cutWeight,"goff");
618     nsel += effh->Integral();
619 buchmann 1.7 cout << __FUNCTION__ << " :: Drawn with " << (TCut(newNegSide.c_str()))*cutWeight << endl;
620 buchmann 1.1 } else {
621     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*cutWeight,"goff");
622 buchmann 1.2
623 buchmann 1.1 nsel += effh->Integral();
624     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewPosSide.c_str()))*cutWeight,"goff");
625     nsel += effh->Integral();
626     events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewNegSide.c_str()))*cutWeight,"goff");
627     nsel -= effh->Integral();
628     }
629     }
630    
631     //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.
632     float normFactor = 1;
633     if(k>=0) get_norm_pdf_factor(events, k, addcut);
634     sel = sel/normFactor;
635     nsel = nsel/normFactor;
636    
637     Float_t tot = Neventsinfile;
638    
639     Value result_wo_signalcont;
640    
641     if(ConsiderSignalContaminationForLimits) {
642     result=(sel-nsel)/tot;
643     resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
644     result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
645     } else {//no signal contamination considered:
646     result=(sel)/tot;
647     resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
648     result_wo_signalcont=Value(result,resulterr);
649     }
650     if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
651 buchmann 1.7 if(!automatized) dout << " MC efficiency: " << result << "+-" << resulterr << " ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " raw MC events out of " << tot << ") raw MC events with PDF normalization factor " << normFactor << std::endl;
652 buchmann 1.1 delete effh;
653     return result_wo_signalcont;
654     }
655    
656    
657    
658     //____________________________________________________________________________________
659     // Selection efficiency for one process (MC)
660     // not in use anymore.
661     /*
662     vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
663     vector<float> process_efficiencies;
664     for(int iprocess=0;iprocess<=10;iprocess++) {
665     float this_process_efficiency,efferr;
666     stringstream addcutplus;
667     addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
668     MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
669     process_efficiencies.push_back(this_process_efficiency);
670     }
671     return process_efficiencies;
672     }
673     */
674    
675     void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
676     TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
677     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
678     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
679     TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
680     Int_t bin = hLM4->FindBin(jzbSel); // To get the error
681     jzbeff=Interpolate(jzbSel,hLM4);
682     jzbefferr=hLM4->GetBinError(bin);
683     if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
684     if(!automatized) dout << " " << jzbeff << "+-" << jzbefferr << std::endl;
685     }
686    
687     //________________________________________________________________________
688     // Effect of energy scale on efficiency
689 buchmann 1.3 void JZBjetScale(TTree *events, bool domSUGRA, map < pair<float, float>, map<string, float> > xsec, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",Float_t jzbSelection=-1, TString plotName = "" ) {
690 buchmann 1.1 TCut kbase(genMassCut&&"genZPt>0");
691     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
692     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
693     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
694    
695     TCut ksel(cutmass&&cutOSSF);
696    
697     if ( !(plotName.Length()>1) ) plotName = informalname;
698    
699 buchmann 1.3 TCut susyweight("1.0");
700     if(domSUGRA) {
701     stringstream susyweightS;
702     float sumweights=0;
703     susyweightS << "((";
704     for(int i=0;i<11;i++) {
705     if(i==0) susyweightS << "(";
706     if(i>0) susyweightS << " + ";
707     susyweightS << "(process==" << i << ")*";
708     float thisxsec=XSForProcessViaAddCutWrapper(addcut,xsec,i);
709     susyweightS << thisxsec;
710     sumweights+=thisxsec;
711     if(i==10) susyweightS << ")";
712     }
713     susyweightS << ")/" << sumweights << ")";
714     susyweight=TCut(susyweightS.str().c_str());
715     }
716    
717 buchmann 1.1 nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
718 buchmann 1.4 TH1F* hist = plotEff(events,((kbase&&ksel&&cutnJets)*susyweight),informalname,flipped);
719     TH1F* histp = plotEff(events,((kbase&&ksel&&cutnJetsJESup)*susyweight),informalname,flipped);
720     TH1F* histm = plotEff(events,((kbase&&ksel&&cutnJetsJESdown)*susyweight),informalname,flipped);
721 buchmann 1.1
722     // Dump some information
723     Float_t eff = Interpolate(jzbSel,hist);
724     Float_t effp = Interpolate(jzbSel,histp);
725     Float_t effm = Interpolate(jzbSel,histm);
726     if(!automatized) dout << " Efficiency at JZB==" << jzbSel << std::endl;
727     if(!automatized) dout << " JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
728     if(!automatized) dout << " central: " << eff << std::endl;
729     if(!automatized) dout << " JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
730     jesup=(effp-eff)/eff;
731     jesdown=(effm-eff)/eff;
732 buchmann 1.4 if(eff==0) {
733     jesup=1;
734     jesdown=1;
735     }
736 buchmann 1.1 }
737    
738     //________________________________________________________________________
739     // Effect of energy scale on JZB efficiency
740 buchmann 1.3 void doJZBscale(TTree *events, bool domSUGRA, map < pair<float, float>, map<string, float> > xsec, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
741 buchmann 1.1
742     TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
743     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
744     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
745     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
746     TCut ksel(cutmass&&cutOSSF);
747    
748     nBins = 50;
749     jzbMin = 0.5*jzbSel;
750     jzbMax = 2.0*jzbSel;
751    
752 buchmann 1.3 TCut susyweight("1.0");
753     if(domSUGRA) {
754     stringstream susyweightS;
755     float sumweights=0;
756     susyweightS << "((";
757     for(int i=0;i<11;i++) {
758     if(i==0) susyweightS << "(";
759     if(i>0) susyweightS << " + ";
760     susyweightS << "(process==" << i << ")*";
761     float thisxsec=XSForProcessViaAddCutWrapper(addcut,xsec,i);
762     susyweightS << thisxsec;
763     sumweights+=thisxsec;
764     if(i==10) susyweightS << ")";
765     }
766     susyweightS << ")/" << sumweights << ")";
767     susyweight=TCut(susyweightS.str().c_str());
768     }
769    
770    
771    
772     TH1F* hist = plotEff(events,((kbase&&ksel)*susyweight),informalname,flipped);
773 buchmann 1.1
774     // Dump some information
775     Float_t eff = Interpolate(jzbSel,hist);
776     Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
777     Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
778     if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.+systematic) << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)" << std::endl;
779     if(!automatized) dout << " efficiency at JZB==" << jzbSel << ": " << eff << std::endl;
780     if(!automatized) dout << " efficiency at JZB==" << jzbSel*(1.-systematic) << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)" << std::endl;
781     up=((effp-eff)/eff);
782     down=((effm-eff)/eff);
783 buchmann 1.4 if(eff==0) {
784     up=1;
785     down=1;
786     }
787 buchmann 1.1 }
788    
789     //________________________________________________________________________
790     // JZB response (true/reco. vs. true)
791     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 ) {
792    
793     jzbMin = 20;
794     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
795     TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
796     if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
797     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
798     if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
799     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
800     TCut ksel(cutmass&&cutOSSF);
801    
802     TProfile* hJzbResp = new TProfile("hJzbResp","JZB response ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
803    
804     string locmcjzbexpression=mcjzbexpression;
805     if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
806     string possibleminus="";
807     if(flipped>0) possibleminus="-";
808     if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
809     else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
810    
811     hJzbResp->SetMaximum(1.2);
812     hJzbResp->SetMinimum(0.2);
813     hJzbResp->Fit("pol0","Q");
814     TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
815     if(!fittedfunction) {
816     // in case there are not enough points passing our selection
817     cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
818     resp=1;
819     resperr=1;
820     } else {
821     resp=fittedfunction->GetParameter(0);
822     resperr=fittedfunction->GetParError(0);
823     if(!automatized) dout << " Response: " << resp << " +/- " << resperr << endl;
824     }
825     delete hJzbResp;
826     }
827    
828    
829     //________________________________________________________________________________________
830     // PDF uncertainty
831 buchmann 1.3 float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int scantype, map < pair<float, float>, map<string, float> > xsec, int NPdfs, string addcut="") {
832 buchmann 1.1 std::vector<float> efficiency;
833     for(int k = 1; k < NPdfs; k++) {
834     float result, resulterr;
835     Value flipval;
836 buchmann 1.3 MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, scantype, xsec, addcut, k);
837 buchmann 1.1 efficiency.push_back(result);
838     }
839     float errHi, errLow,err;
840     master_formula(efficiency, errHi, errLow);
841     err=errLow;
842     if(errHi>errLow) err=errHi;
843     if(!automatized) dout << " Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
844     return err;
845    
846     }
847    
848     int get_npdfs(TTree *events) {
849     int NPDFs;
850     events->SetBranchAddress("NPdfs",&NPDFs);
851     events->GetEntry(1);
852     return NPDFs;
853     }
854    
855    
856 buchmann 1.3 void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,int flipped, map < pair<float, float>, map<string, float> > xsec, string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool isscan=false,int scantype=-1) {
857 buchmann 1.5 float JZBScaleUncert=0.05; //if you modify this value please also adapt it in ShapeLimit.C not only here in Systematics.C
858 buchmann 1.1 mcjzbexpression=mcjzb;
859     float triggereff=2.0/100;// in range [0,1]
860     dout << "Trigger efficiency not implemented in this script yet, still using external one" << endl;
861     float leptonseleff=2.0/100;// in range [0,1]
862     leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
863     dout << "Lepton selection efficiency not implemented in this script yet, still using external one" << endl;
864    
865     int NPdfs=0;
866     if(isscan) NPdfs = get_npdfs(events);
867    
868     float mceff,mcefferr,jzbeff,jzbefferr;
869     if(!automatized) dout << "MC efficiencies:" << endl;
870     Value flipefficiency;
871 buchmann 1.3 Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,scantype,xsec,addcut,-1);
872 buchmann 1.1 if(!automatized) cout << " Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
873    
874     if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
875     if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
876    
877     if(!automatized) dout << "Error from Peak position:";
878     float sysfrompeak=0;
879     PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
880 buchmann 1.3
881 buchmann 1.5 bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in ShapeLimit.C not only here in Systematics.C
882 buchmann 1.3
883 buchmann 1.1 if(!automatized) dout << "Jet energy scale (JES): " << std::endl;
884     float jesup,jesdown;
885 buchmann 1.3 JZBjetScale(events,scantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,jesdown,jesup,informalname,flipped,requireZ,addcut);
886 buchmann 1.1
887     if(!automatized) dout << "JZB scale: " << std::endl;
888     float scaleup,scaledown,scalesyst;
889 buchmann 1.3 doJZBscale(events,scantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
890 buchmann 1.1
891     if(!automatized) dout << "JZB response: " << std::endl;
892     float resp,resperr;
893     if(PlottingSetup::computeJZBresponse) {
894     if(!automatized) dout << "JZB response: " << std::endl;
895     if(!isscan) JZBresponse(events,informalname,requireZ,resp,resperr,flipped,addcut);
896     }
897    
898     if(!automatized) dout << "Pileup: " << std::endl;
899     // float resolution;
900     //resolution=pileup(events,requireZ,informalname,flipped,addcut);
901    
902     float PDFuncert=0;
903     if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
904 buchmann 1.3 if(isscan&&scantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, scantype, xsec, NPdfs, addcut);//for mSUGRA this is done differently!
905 buchmann 1.1
906     dout << "_______________________________________________" << endl;
907     dout << " SUMMARY FOR " << informalname << " with JZB>" << jzbSel << " (all in %) ";
908     if(addcut!="") dout << "With additional cut: " << addcut;
909     dout << endl;
910     dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
911     dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
912     dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
913     dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
914     dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
915     // dout << "Resolution : " << resolution << endl; // in range [0,1]
916     dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
917     if(isscan) dout << "PDF uncertainty : " << PDFuncert << endl; // in range [0,1]
918     if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
919     if(PlottingSetup::computeJZBresponse)dout << "JZB response : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
920    
921     float toterr=0;
922     toterr+=(triggereff)*(triggereff);
923     toterr+=(leptonseleff)*(leptonseleff);
924     if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
925     if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
926     // toterr+=(resolution*resolution);
927     toterr+=(sysfrompeak*sysfrompeak);
928     if(isscan) toterr+=(PDFuncert*PDFuncert);
929 buchmann 1.4 toterr=TMath::Sqrt(toterr);
930     dout << "TOTAL SYSTEMATICS: " << toterr << " --> " << toterr*mceff << endl;
931     float systerr=toterr*mceff;
932     toterr=TMath::Sqrt(toterr*toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
933 buchmann 1.1
934     dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst) %" << endl;
935     dout << " we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
936     dout << "_______________________________________________" << endl;
937    
938     //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
939     vector<float> res;
940     res.push_back(jzbSel);
941     res.push_back(mceff);
942     res.push_back(mcefferr);
943 buchmann 1.4 res.push_back(systerr);
944 buchmann 1.1 res.push_back(toterr);
945     if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
946     if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
947     // res.push_back(fabs(resolution));
948     res.push_back(0.0);
949     res.push_back(mceff_nosigcont.getValue());
950     res.push_back(mceff_nosigcont.getError());
951     if(isscan) res.push_back(PDFuncert);
952     results.push_back(res);
953     }
954    
955     vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
956     automatized=true;
957     vector< vector<float> > systematics;
958 buchmann 1.3 map < pair<float, float>, map<string, float> > xsec; // only needed for mSUGRA, and in that case this funciton isn't called
959    
960 buchmann 1.6 for (int isignal=0; isignal<(int)signalsamples.collection.size();isignal++) {
961 buchmann 1.1 dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
962 buchmann 1.6 for(int ibin=0;ibin<(int)bins.size();ibin++) {
963 buchmann 1.1 jzbSel=bins[ibin];
964     geqleq="geq";
965 buchmann 1.3 do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,xsec,mcjzb,datajzb,mcpeakerror,requireZ);
966 buchmann 1.1 }//end of bin loop
967     }//end of signal loop
968     return systematics;
969     }