ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.3
Committed: Mon Apr 16 09:57:00 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.2: +8 -8 lines
Log Message:
Making stat uncertainties independently

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4     #include <fstream>
5    
6     #include <TCut.h>
7     #include <TROOT.h>
8     #include <TCanvas.h>
9     #include <TMath.h>
10     #include <TColor.h>
11     #include <TPaveText.h>
12     #include <TRandom.h>
13     #include <TH1.h>
14     #include <TH2.h>
15     #include <TF1.h>
16     #include <TSQLResult.h>
17     #include <TProfile.h>
18     #include <TSystem.h>
19     #include <TRandom3.h>
20    
21     //#include "TTbar_stuff.C"
22     using namespace std;
23    
24     using namespace PlottingSetup;
25    
26     namespace SUSYScanSpace {
27     vector<string> loaded_files;
28     int SUSYscantype;
29     float SavedMGlu;
30     float SavedMLSP;
31     string SavedMLSPname;
32     string SavedMGluname;
33     }
34    
35 buchmann 1.2 const char* strip_flip_away(string flipped_name) {
36     int pos = flipped_name.find("flipped_");
37     if(pos>=0&&pos<1000) flipped_name=flipped_name.substr(8,1000);
38     return flipped_name.c_str();
39     }
40    
41     void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF, int flipped) {
42     TH1F *dataob;
43     if(flipped) dataob = (TH1F*)f->Get("flipped_data_obs");
44     else dataob = (TH1F*)f->Get("data_obs");
45     TH1F *signal;
46    
47     if(flipped) signal = (TH1F*)f->Get("flipped_signal");
48     else signal = (TH1F*)f->Get("signal");
49    
50     TH1F *background1;
51     if(flipped) background1 = (TH1F*)f->Get("flipped_TTbarBackground");
52     else background1 = (TH1F*)f->Get("TTbarBackground");
53    
54     TH1F *background2;
55     if(flipped) background2 = (TH1F*)f->Get("flipped_ZJetsBackground");
56     else background2 = (TH1F*)f->Get("ZJetsBackground");
57    
58     assert(dataob);
59     assert(signal);
60     assert(background1);
61     assert(background2);
62    
63     ofstream datacard;
64 buchmann 1.3 write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
65 buchmann 1.2 datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
66     datacard << "Writing this to a file.\n";
67     datacard << "imax 1\n"; // number of channels
68     datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
69     datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
70     datacard << "---------------\n";
71     datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
72     datacard << "---------------\n";
73     datacard << "bin 1\n";
74     datacard << "observation "<<dataob->Integral()<<"\n";
75     datacard << "------------------------------\n";
76     datacard << "bin 1 1 1\n";
77     datacard << "process signal TTbarBackground ZJetsBackground\n";
78     datacard << "process 0 1 2\n";
79     datacard << "rate "<<signal->Integral()<<" "<<background1->Integral()<<" "<<background2->Integral()<<"\n";
80     datacard << "--------------------------------\n";
81 buchmann 1.3 datacard << "lumi lnN " << 1+ PlottingSetup::lumiuncert << " - - luminosity uncertainty\n"; // only affects MC -> signal!
82     datacard << "Stat shape 1 - - statistical uncertainty (signal)\n";
83     datacard << "Stat shape - 1 - statistical uncertainty (ttbar)\n";
84     datacard << "Stat shape - - 1 statistical uncertainty (zjets)\n";
85     datacard << "Sys shape - - 1 systematic uncertainty on ttbar\n";
86     datacard << "Sys shape - 1 - systematic uncertainty on zjets\n";
87 buchmann 1.2 datacard << "JES shape 1 - - uncertainty on Jet Energy Scale\n";
88     datacard << "JSU lnN " << 1+uJSU << " - - JZB Scale Uncertainty\n";
89     if(uPDF>0) datacard << "PDF lnN " << 1+uPDF << " - - uncertainty from PDFs\n";
90    
91 buchmann 1.3 // datacard << "peak shape 1 1 uncertainty on signal resolution.
92 buchmann 1.2 datacard.close();
93     }
94    
95     float QuickDrawNevents=0;
96    
97 buchmann 1.1 TH1F* QuickDraw(TTree *events, string hname, string variable, vector<float> binning, string xlabel, string ylabel, TCut cut, string addcut, bool dataormc, float luminosity, map < pair<float, float>, map<string, float> > &xsec) {
98     TH1F *histo = new TH1F(hname.c_str(),hname.c_str(),binning.size()-1,&binning[0]);
99     histo->Sumw2();
100     stringstream drawcommand;
101     drawcommand << variable << ">>" << hname;
102     if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
103     events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
104     } else {
105     float totalxs=0;
106     for(int i=0;i<12;i++) {
107     stringstream drawcommand2;
108     drawcommand2 << variable << ">>+" << hname;
109     stringstream mSUGRAweight;
110     float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
111     totalxs+=xs;
112     events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
113     }
114     histo->Scale(1.0/totalxs);
115     }
116     events->Draw("eventNum",addcut.c_str(),"goff");
117     float nevents = events->GetSelectedRows();
118 buchmann 1.2 QuickDrawNevents=nevents;
119 buchmann 1.1 histo->Scale(luminosity/nevents); // this will give a histogram which is normalized to 1 pb so that any limit will be with respect to 1 pb
120    
121     return histo;
122     }
123    
124     /*void do_stat_up(TH1F *h, float sign=1.0) {
125     for(int i=1;i<=h->GetNbinsX();i++) {
126     h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
127     }
128     h->Write();
129     }
130    
131     void do_stat_dn(TH1F *h) {
132     do_stat_up(h,-1.0);
133     }*/
134    
135     void SQRT(TH1F *h) {
136     for (int i=1;i<=h->GetNbinsX();i++) {
137     h->SetBinContent(i,TMath::Sqrt(h->GetBinContent(i)));
138     }
139     }
140    
141     TH1F* Multiply(TH1F *h1, TH1F *h2) {
142     TH1F *h = (TH1F*)h1->Clone();
143     for(int i=1;i<=h1->GetNbinsX();i++) {
144     h->SetBinContent(i,h1->GetBinContent(i) * h2->GetBinContent(i));
145     }
146     return h;
147     }
148    
149     void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map < pair<float, float>, map<string, float> > &xsec) {
150     binning.push_back(8000);
151     limfile->cd();
152     dout << "Creatig shape templates ";
153     if(identifier!="") dout << "for "<<identifier;
154     dout << " ... " << endl;
155    
156     TCut limitnJetcut;
157     if(JES==noJES) limitnJetcut=cutnJets;
158     else {
159     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
160     if(JES==JESup) limitnJetcut=cutnJetsJESup;
161     }
162    
163     stringstream mSUGRAweight;
164     if(SUSYScanSpace::SUSYscantype==mSUGRA) {
165     //for(int i<
166     mSUGRAweight << "bla";
167     } else mSUGRAweight << "(1)";
168    
169     write_warning(__FUNCTION__,"Not done yet for mSUGRA");
170    
171     float zjetsestimateuncert=zjetsestimateuncertONPEAK;
172     float emuncert=emuncertONPEAK;
173     float emsidebanduncert=emsidebanduncertONPEAK;
174     float eemmsidebanduncert=eemmsidebanduncertONPEAK;
175    
176     if(!PlottingSetup::RestrictToMassPeak) {
177     zjetsestimateuncert=zjetsestimateuncertOFFPEAK;
178     emuncert=emuncertOFFPEAK;
179     emsidebanduncert=emsidebanduncertOFFPEAK;
180     eemmsidebanduncert=eemmsidebanduncertOFFPEAK;
181     }
182    
183     if(signalonly) {
184 buchmann 1.2 cout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: " << identifier << ")" << endl;
185 buchmann 1.1 TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
186     TH1F *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
187     TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
188     TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
189    
190     TH1F *SBOSSFP;
191     TH1F *SBOSOFP;
192     TH1F *SBOSSFN;
193     TH1F *SBOSOFN;
194    
195     if(PlottingSetup::RestrictToMassPeak) {
196     SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
197     SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
198     SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
199     SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
200     }
201    
202     string signalname="signal";
203     if(identifier!="") {
204     signalname="signal_"+identifier;
205     }
206    
207     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
208     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
209    
210     TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
211     TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
212    
213     Lobs->Add(ZOSSFP);
214     Lpred->Add(ZOSSFN);
215    
216 buchmann 1.2 cout << "SITUATION FOR SIGNAL: " << endl;
217     cout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << " JZB < 0 :" << ZOSSFN->Integral() << endl;
218     cout << " OSOF JZB> 0 : " << ZOSOFP->Integral() << " JZB < 0 :" << ZOSOFN->Integral() << endl;
219     if(PlottingSetup::RestrictToMassPeak) {
220     cout << " OSSF SB JZB> 0 : " << SBOSSFP->Integral() << " JZB < 0 :" << SBOSSFN->Integral() << endl;
221     cout << " OSOF SB JZB> 0 : " << SBOSOFP->Integral() << " JZB < 0 :" << SBOSOFN->Integral() << endl;
222     }
223    
224    
225 buchmann 1.1 flippedLobs->Add(ZOSSFN);
226     flippedLpred->Add(ZOSSFP);
227    
228     if(PlottingSetup::RestrictToMassPeak) {
229     Lpred->Add(ZOSOFP,1.0/3);
230     Lpred->Add(ZOSOFN,-1.0/3);
231     Lpred->Add(SBOSSFP,1.0/3);
232     Lpred->Add(SBOSSFN,-1.0/3);
233     Lpred->Add(SBOSOFP,1.0/3);
234     Lpred->Add(SBOSOFN,-1.0/3);
235    
236     //flipped prediction
237     flippedLpred->Add(ZOSOFP,-1.0/3);
238     flippedLpred->Add(ZOSOFN,1.0/3);
239     flippedLpred->Add(SBOSSFP,-1.0/3);
240     flippedLpred->Add(SBOSSFN,1.0/3);
241     flippedLpred->Add(SBOSOFP,-1.0/3);
242     flippedLpred->Add(SBOSOFN,1.0/3);
243    
244     } else {
245     Lpred->Add(ZOSOFP,1.0);
246     Lpred->Add(ZOSOFN,-1.0);
247 buchmann 1.2
248 buchmann 1.1 //flipped prediction
249     flippedLpred->Add(ZOSOFP,-1.0);
250     flippedLpred->Add(ZOSOFN,1.0);
251     }
252    
253 buchmann 1.2 TH1F *signal = (TH1F*)Lobs->Clone("signal");
254 buchmann 1.1 signal->Add(Lpred,-1);
255     signal->SetName(signalname.c_str());
256 buchmann 1.2 signal->SetTitle(signalname.c_str());
257 buchmann 1.1 signal->Write();
258    
259     TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
260     flippedsignal->Add(flippedLpred,-1);
261     flippedsignal->SetName(("flipped_"+signalname).c_str());
262     flippedsignal->Write();
263    
264     if(identifier=="") {
265     TH1F *signalStatUp = (TH1F*)signal->Clone();
266     signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
267     signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
268    
269     TH1F *signalStatDn = (TH1F*)signal->Clone();
270     signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
271     signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
272    
273     for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
274 buchmann 1.2 float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
275     signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
276     signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
277     signal->SetBinError(i,staterr);
278 buchmann 1.1 }
279    
280     signalStatDn->Write();
281     signalStatUp->Write();
282    
283     delete signalStatDn;
284     delete signalStatUp;
285    
286     TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
287     flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
288     flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
289    
290     TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
291     flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
292     flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
293    
294     for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
295 buchmann 1.2 float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
296     flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
297     flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
298     flippedsignal->SetBinError(i,staterr);
299 buchmann 1.1 }
300    
301     flippedsignalStatDn->Write();
302     flippedsignalStatUp->Write();
303    
304     delete flippedsignalStatDn;
305     delete flippedsignalStatUp;
306     }
307    
308     delete ZOSSFP;
309     delete ZOSOFP;
310     delete ZOSSFN;
311     delete ZOSOFN;
312    
313     if(PlottingSetup::RestrictToMassPeak) {
314     delete SBOSSFP;
315     delete SBOSOFP;
316     delete SBOSSFN;
317     delete SBOSOFN;
318     }
319    
320     delete Lobs;
321     delete Lpred;
322     delete flippedLobs;
323     delete flippedLpred;
324     }
325    
326    
327     if(dataonly) {
328     cout << "Processing data with datajzb: " << datajzb << endl;
329     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
330     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
331     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
332     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
333    
334     TH1F *SBOSSFP;
335     TH1F *SBOSOFP;
336     TH1F *SBOSSFN;
337     TH1F *SBOSOFN;
338    
339     if(PlottingSetup::RestrictToMassPeak) {
340     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
341     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
342     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
343     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
344     }
345    
346     string obsname="data_obs";
347     string predname="background";
348     string Zpredname="ZJetsBackground";
349     string Tpredname="TTbarBackground";
350     if(identifier!="") {
351     obsname=("data_"+identifier);
352     predname=("background_"+identifier);
353     Zpredname=("ZJetsBackground_"+identifier);
354     Tpredname=("TTbarBackground_"+identifier);
355     }
356    
357     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
358     obs->SetName(obsname.c_str());
359     obs->Write();
360    
361     TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
362     flippedobs->SetName(("flipped_"+obsname).c_str());
363     flippedobs->Write();
364    
365     TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
366     TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
367    
368     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
369    
370     TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
371     TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
372    
373     TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
374     if(PlottingSetup::RestrictToMassPeak) {
375     pred->Add(ZOSOFP,1.0/3);
376     pred->Add(ZOSOFN,-1.0/3);
377     pred->Add(SBOSSFP,1.0/3);
378     pred->Add(SBOSSFN,-1.0/3);
379     pred->Add(SBOSOFP,1.0/3);
380     pred->Add(SBOSOFN,-1.0/3);
381    
382     Tpred->Add(ZOSSFN,-1.0);
383     Tpred->Add(ZOSOFP,1.0/3);
384     Tpred->Add(SBOSSFP,1.0/3);
385     Tpred->Add(SBOSOFP,1.0/3);
386    
387     Zpred->Add(ZOSOFN,-1.0/3);
388     Zpred->Add(SBOSSFN,-1.0/3);
389     Zpred->Add(SBOSOFN,-1.0/3);
390    
391     //flipped
392     flippedpred->Add(ZOSOFP,-1.0/3);
393     flippedpred->Add(ZOSOFN,1.0/3);
394     flippedpred->Add(SBOSSFP,-1.0/3);
395     flippedpred->Add(SBOSSFN,1.0/3);
396     flippedpred->Add(SBOSOFP,-1.0/3);
397     flippedpred->Add(SBOSOFN,1.0/3);
398    
399     flippedTpred->Add(ZOSSFP,-1.0);
400     flippedTpred->Add(ZOSOFN,1.0/3);
401     flippedTpred->Add(SBOSSFN,1.0/3);
402     flippedTpred->Add(SBOSOFN,1.0/3);
403    
404     flippedZpred->Add(ZOSOFP,-1.0/3);
405     flippedZpred->Add(SBOSSFP,-1.0/3);
406     flippedZpred->Add(SBOSOFP,-1.0/3);
407     } else {
408     pred->Add(ZOSOFP,1.0);
409     pred->Add(ZOSOFN,-1.0);
410    
411     Tpred->Add(ZOSSFN,-1.0);
412     Tpred->Add(ZOSOFP,1.0);
413    
414     Zpred->Add(ZOSOFN,-1.0);
415    
416     //flipped
417     flippedpred->Add(ZOSOFN,1.0);
418     flippedpred->Add(ZOSOFP,-1.0);
419    
420     flippedTpred->Add(ZOSSFP,-1.0);
421     flippedTpred->Add(ZOSOFN,1.0);
422    
423     flippedZpred->Add(ZOSOFP,-1.0);
424     }
425    
426     pred->SetName(predname.c_str());
427     Tpred->SetName(Tpredname.c_str());
428     Zpred->SetName(Zpredname.c_str());
429    
430     flippedpred->SetName(("flipped_"+predname).c_str());
431     flippedTpred->SetName(("flipped_"+Tpredname).c_str());
432     flippedZpred->SetName(("flipped_"+Zpredname).c_str());
433    
434     pred->Write();
435     Tpred->Write();
436     Zpred->Write();
437    
438     flippedpred->Write();
439     flippedTpred->Write();
440     flippedZpred->Write();
441    
442     if(identifier=="") {
443     float stretchfactor=1.0;
444     bool isonpeak=false;
445     if(PlottingSetup::RestrictToMassPeak) {
446     stretchfactor=1.0/3.0;
447     isonpeak=true;
448     }
449    
450     //--------------------------------------------------------statistical uncertainty
451    
452     TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
453     predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
454     predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
455     if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
456     if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
457     if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
458     if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
459     SQRT(predstaterr);
460     TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
461     bgStatUp->Add(predstaterr);
462     bgStatUp->Write();
463 buchmann 1.2 TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
464 buchmann 1.1 bgStatDn->Add(predstaterr,-1);
465     bgStatDn->Write();
466     // delete bgStatDn;
467     // delete bgStatUp;
468     delete predstaterr;
469    
470    
471     TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
472     flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
473     flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
474     if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
475     if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
476     if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
477     if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
478     SQRT(flippedpredstaterr);
479     TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
480     fbgStatUp->Add(predstaterr);
481     fbgStatUp->Write();
482 buchmann 1.2 TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
483 buchmann 1.1 fbgStatDn->Add(predstaterr,-1);
484     fbgStatDn->Write();
485     // delete fbgStatDn;
486     // delete fbgStatUp;
487     delete predstaterr;
488    
489     TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
490     Tpredstaterr->Scale(stretchfactor*stretchfactor);
491     if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
492     if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
493     SQRT(Tpredstaterr);
494     TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
495     TpredStatUp->Add(Tpredstaterr);
496     TpredStatUp->Write();
497 buchmann 1.2 TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
498 buchmann 1.1 TpredStatDn->Add(Tpredstaterr,-1);
499     TpredStatDn->Write();
500     // delete TpredStatDn;
501     // delete TpredStatUp;
502     delete Tpredstaterr;
503    
504     TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
505     fTpredstaterr->Scale(stretchfactor*stretchfactor);
506     if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
507     if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
508     SQRT(fTpredstaterr);
509     TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
510     fTpredStatUp->Add(fTpredstaterr);
511     fTpredStatUp->Write();
512 buchmann 1.2 TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
513 buchmann 1.1 fTpredStatDn->Add(fTpredstaterr,-1);
514     fTpredStatDn->Write();
515     // delete fTpredStatDn;
516     // delete fTpredStatUp;
517     delete fTpredstaterr;
518    
519    
520     TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
521     Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
522     if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
523     if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
524     SQRT(Zpredstaterr);
525     TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
526     ZpredStatUp->Add(Zpredstaterr);
527     ZpredStatUp->Write();
528 buchmann 1.2 TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
529 buchmann 1.1 ZpredStatDn->Add(Zpredstaterr,-1);
530     ZpredStatDn->Write();
531     // delete ZpredStatDn;
532     // delete ZpredStatUp;
533     delete Zpredstaterr;
534    
535     TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
536     Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
537     if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
538     if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
539     SQRT(fTpredstaterr);
540     TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
541     fZpredStatUp->Add(fZpredstaterr);
542     fZpredStatUp->Write();
543 buchmann 1.2 TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
544 buchmann 1.1 fZpredStatDn->Add(fZpredstaterr,-1);
545     fZpredStatDn->Write();
546     // delete fZpredStatDn;
547     // delete fZpredStatUp;
548     delete fZpredstaterr;
549    
550     //--------------------------------------------------------systematic uncertainty
551    
552     TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
553     predsyserr->Add(pred,-1);//getting everything to zero.
554     predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
555     predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
556     predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
557     if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
558     if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
559     if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
560     if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
561     SQRT(predsyserr);
562     TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
563     bgSysUp->Add(predsyserr);
564     bgSysUp->Write();
565 buchmann 1.2 TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
566 buchmann 1.1 bgSysDn->Add(predsyserr,-1);
567     bgSysDn->Write();
568     delete predsyserr;
569    
570     TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
571     fpredsyserr->Add(pred,-1);//getting everything to zero.
572     fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
573     fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
574     fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
575     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
576     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
577     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
578     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
579     SQRT(fpredsyserr);
580     TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
581     fbgSysUp->Add(fpredsyserr);
582     fbgSysUp->Write();
583 buchmann 1.2 TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
584 buchmann 1.1 fbgSysDn->Add(fpredsyserr,-1);
585     fbgSysDn->Write();
586     delete fpredsyserr;
587    
588    
589     TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
590     Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
591     Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
592     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
593     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
594     SQRT(Tpredsyserr);
595     TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
596     TpredSysUp->Add(Tpredsyserr);
597     TpredSysUp->Write();
598 buchmann 1.2 TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
599 buchmann 1.1 TpredSysDn->Add(Tpredsyserr,-1);
600     TpredSysDn->Write();
601     delete Tpredsyserr;
602    
603     TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
604     fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
605     fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
606     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
607     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
608     SQRT(fTpredsyserr);
609     TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
610     fTpredSysUp->Add(fTpredsyserr);
611     fTpredSysUp->Write();
612 buchmann 1.2 TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
613 buchmann 1.1 fTpredSysDn->Add(fTpredsyserr,-1);
614     fTpredSysDn->Write();
615     delete fTpredsyserr;
616    
617    
618     //------------------------------------------------------------------------------------------------------------------------
619     TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
620     Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
621     Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
622     Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
623     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
624     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
625     SQRT(Zpredsyserr);
626     TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
627     ZpredSysUp->Add(Zpredsyserr);
628     ZpredSysUp->Write();
629 buchmann 1.2 TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
630 buchmann 1.1 ZpredSysDn->Add(Zpredsyserr,-1);
631     ZpredSysDn->Write();
632     delete Zpredsyserr;
633    
634     TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
635     fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
636     fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
637     fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
638     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
639     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
640     SQRT(fZpredsyserr);
641     TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
642     fZpredSysUp->Add(fZpredsyserr);
643     fZpredSysUp->Write();
644 buchmann 1.2 TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
645 buchmann 1.1 fZpredSysDn->Add(fZpredsyserr,-1);
646     fZpredSysDn->Write();
647     delete fZpredsyserr;
648     }
649    
650 buchmann 1.2 /*if(identifier=="") {
651 buchmann 1.1 for(int i=0;i<binning.size()-1;i++) {
652     cout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
653     }
654 buchmann 1.2 }*/
655 buchmann 1.1 delete ZOSSFP;
656     delete ZOSOFP;
657     delete ZOSSFN;
658     delete ZOSOFN;
659    
660     if(PlottingSetup::RestrictToMassPeak) {
661     delete SBOSSFP;
662     delete SBOSOFP;
663     delete SBOSSFN;
664     delete SBOSOFN;
665     }
666     }
667    
668     }
669    
670     ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
671     map < pair<float, float>, map<string, float> > xsec;
672     if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
673    
674     write_info(__FUNCTION__,"Implementing the shape function");
675     time_t timestamp;
676     tm *now;
677     timestamp = time(0);
678     now = localtime(&timestamp);
679     stringstream RunDirectoryS;
680     RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
681     string RunDirectory = RunDirectoryS.str();
682    
683     ensure_directory_exists(RunDirectory);
684    
685     TFile *datafile = new TFile("../StoredShapes.root","READ");
686     if(datafile->IsZombie()) {
687     write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
688     assert(!datafile->IsZombie());
689     }
690     cout << "Run Directory: " << RunDirectory << endl;
691 buchmann 1.2 TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
692 buchmann 1.1
693     TIter nextkey(datafile->GetListOfKeys());
694     TKey *key;
695     while ((key = (TKey*)nextkey()))
696     {
697     TH1F *obj =(TH1F*) key->ReadObj();
698     limfile->cd();
699     obj->Write();
700     }
701    
702     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
703    
704     bool signalonly=true;
705     bool dataonly=false;
706    
707     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
708 buchmann 1.2 // generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
709     // generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
710 buchmann 1.1 generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
711     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
712    
713     float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
714     float scaledown=0,scaleup=0,scalesyst=0;
715     float mceff=0,mcefferr=0;
716     int flipped=0;
717     int Neventsinfile=10000;//doesn't matter.
718     string informalname="ShapeAnalysis";
719    
720     bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
721    
722     MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
723     if(mceff<0) flipped=1;
724     doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
725     float PDFuncert;
726     int NPdfs = get_npdfs(events);
727     if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
728    
729    
730 buchmann 1.2 float JZBscale=scaledown;
731     if(scaleup>scaledown) JZBscale=scaleup;
732     dout << endl;
733     dout << endl;
734     dout << "Will use the following scalar (non-shape) systematics : " << endl;
735     dout << " JZB Scale (JSU) : " << JZBscale << endl;
736     dout << " PDF : " << PDFuncert << endl;
737    
738     prepare_limit_datacard(RunDirectory,limfile,JZBscale,PDFuncert,flipped);
739     TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
740    
741     TIter nnextkey(limfile->GetListOfKeys());
742     TKey *nkey;
743     while ((nkey = (TKey*)nnextkey()))
744     {
745     TH1F *obj =(TH1F*) nkey->ReadObj();
746     if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
747     if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
748     if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
749     final_limfile->cd();
750     obj->Write();
751     }
752 buchmann 1.1
753 buchmann 1.2 final_limfile->Close();
754 buchmann 1.1 limfile->Close();
755 buchmann 1.2 write_info(__FUNCTION__,"Shape root file and datacard have been generated in "+RunDirectory);
756    
757     int CreatedModelFileExitCode = gSystem->Exec(("bash CreateModel.sh "+RunDirectory+" susydatacard.txt").c_str());
758     cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
759     assert(CreatedModelFileExitCode==0);
760     ShapeDroplet alpha;
761     alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
762     alpha.PDF=PDFuncert;
763     alpha.JSU=JZBscale;
764     dout << "Done reading limit information from droplet" << endl;
765     dout << alpha << endl;
766 buchmann 1.1
767     write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up.");
768     dout << "Everything is saved in " << RunDirectory << endl;
769     // gSystem->Exec(("rm -r "+RunDirectory).c_str());
770    
771     return alpha;
772     }
773    
774     void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
775 buchmann 1.2 // dout << "Generating data shape templates ... " << endl;
776 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
777     TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
778     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
779     TTree *faketree;
780     bool dataonly=true;
781     bool signalonly=false;
782     string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
783     float jzbpeakerrormc=0;
784     generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
785     // don't need these effects for obs & pred, only for signal!
786     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
787     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
788     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
789     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
790     datafile->Close();
791     }
792