ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.10
Committed: Mon May 7 16:44:27 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.9: +62 -45 lines
Log Message:
made CreateModel paths absolute (necessary to run on the batch system); added some checks when filling histos

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 buchmann 1.6 void EliminateNegativeEntries(TH1F *histo) {
42     for(int i=1;i<=histo->GetNbinsX();i++) {
43     if(histo->GetBinContent(i)<0) histo->SetBinContent(i,0);
44     }
45     }
46    
47     void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
48     TH1F *dataob = (TH1F*)f->Get("data_obs");
49     TH1F *signal = (TH1F*)f->Get("signal");
50     TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
51     TH1F *Zbackground = (TH1F*)f->Get("ZJetsBackground");
52 buchmann 1.2
53     assert(dataob);
54     assert(signal);
55 buchmann 1.6 assert(Tbackground);
56     assert(Zbackground);
57 buchmann 1.2
58    
59 buchmann 1.9 if(FullMCAnalysis) {
60     //dealing with MC analysis - need to use shapes.
61     write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
62     ofstream datacard;
63     write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
64     datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
65     datacard << "Writing this to a file.\n";
66     datacard << "imax 1\n"; // number of channels
67     datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
68     datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
69     datacard << "---------------\n";
70     datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
71     datacard << "---------------\n";
72     datacard << "bin 1\n";
73     datacard << "observation "<<dataob->Integral()<<"\n";
74     datacard << "------------------------------\n";
75     datacard << "bin 1 1 1\n";
76     datacard << "process signal TTbarBackground ZJetsBackground\n";
77     datacard << "process 0 1 2\n";
78     datacard << "rate "<<signal->Integral()<<" "<<Tbackground->Integral()<<" "<<Zbackground->Integral()<<"\n";
79     datacard << "--------------------------------\n";
80     datacard << "lumi lnN " << 1+ PlottingSetup::lumiuncert << " - - luminosity uncertainty\n"; // only affects MC -> signal!
81     datacard << "Stat shape - 1 - statistical uncertainty (ttbar)\n";
82     datacard << "Stat shape - - 1 statistical uncertainty (zjets)\n";
83     datacard << "Sys shape - 1 - systematic uncertainty on ttbar\n";
84     datacard << "Sys shape - - 1 systematic uncertainty on zjets\n";
85     datacard << "Stat shape 1 - - statistical uncertainty (signal)\n";
86     datacard << "JES shape 1 - - uncertainty on Jet Energy Scale\n";
87     datacard << "JSU lnN " << 1+uJSU << " - - JZB Scale Uncertainty\n";
88     if(uPDF>0) datacard << "PDF lnN " << 1+uPDF << " - - uncertainty from PDFs\n";
89     datacard.close();
90     } else {
91     //doing mutibin analysis
92     ofstream datacard;
93     write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
94     datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
95     datacard << "Writing this to a file.\n";
96     datacard << "imax " << dataob->GetNbinsX() << "\n"; // number of channels
97     datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
98     datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
99     datacard << "---------------\n";
100     datacard << "bin\t";
101     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << " \t";
102     datacard << "\n";
103     datacard << "observation\t";
104     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinContent(i) << " \t";
105     datacard<<"\n";
106     datacard << "------------------------------\n";
107     datacard << "bin\t";
108     for(int i=1;i<=dataob->GetNbinsX();i++) {
109     datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
110     " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
111     " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t";
112     }
113     datacard << "\n";
114     datacard << "process\t";
115     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t signal \t TTbarBackground \t ZJetsBackground";
116     datacard << "\n";
117     datacard << "process\t";
118     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t 0 \t 1 \t 2";
119     datacard << "\n";
120    
121    
122     datacard << "rate\t";
123     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t " << signal->GetBinContent(i) << " \t " << Tbackground->GetBinContent(i) << " \t " << Zbackground->GetBinContent(i) << "\t";
124     datacard<<"\n";
125    
126     datacard << "--------------------------------\n";
127    
128    
129     datacard << "lumi lnN \t";
130     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+ PlottingSetup::lumiuncert << "\t - \t -";
131     datacard << " luminosity uncertainty\n"; // only affects MC -> signal!
132    
133     // Statistical uncertainty
134     datacard << "Stat lnN \t";
135     for(int i=1;i<=dataob->GetNbinsX();i++) {
136     //Signal
137     float central = signal->GetBinContent(i);
138     float up = ((TH1F*)f->Get("signal_StatDown"))->GetBinContent(i);
139     float down = ((TH1F*)f->Get("signal_StatUp"))->GetBinContent(i);
140     float suncert=0;
141     if(central>0) {
142     if(abs(up-central)>abs(down-central)) suncert=abs(up-central)/central;
143     else suncert=abs(central-down)/central;
144     }
145    
146     //TTbar
147     central = Tbackground->GetBinContent(i);
148     up = ((TH1F*)f->Get("TTbarBackground_StatUp"))->GetBinContent(i);
149     down = ((TH1F*)f->Get("TTbarBackground_StatDown"))->GetBinContent(i);
150     float tuncert=0;
151     if(central>0) {
152     if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
153     else tuncert=abs(central-down)/central;
154     }
155     //ZJets
156     central = Zbackground->GetBinContent(i);
157     up = ((TH1F*)f->Get("ZJetsBackground_StatUp"))->GetBinContent(i);
158     down = ((TH1F*)f->Get("ZJetsBackground_StatDown"))->GetBinContent(i);
159     float zuncert=0;
160     if(central>0) {
161     if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
162     else zuncert=abs(central-down)/central;
163     }
164     datacard << " " << 1+suncert << " \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
165     }
166    
167     datacard << " statistical uncertainty\n";
168    
169     // Statistical uncertainty
170     datacard << "Sys lnN \t";
171     for(int i=1;i<=dataob->GetNbinsX();i++) {
172     float central = Tbackground->GetBinContent(i);
173     float up = ((TH1F*)f->Get("TTbarBackground_SysUp"))->GetBinContent(i);
174     float down = ((TH1F*)f->Get("TTbarBackground_SysDown"))->GetBinContent(i);
175     float tuncert=0;
176     if(central>0) {
177     if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
178     else tuncert=abs(central-down)/central;
179     }
180     central = Zbackground->GetBinContent(i);
181     up = ((TH1F*)f->Get("ZJetsBackground_SysUp"))->GetBinContent(i);
182     down = ((TH1F*)f->Get("ZJetsBackground_SysDown"))->GetBinContent(i);
183     float zuncert=0;
184     if(central>0) {
185     if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
186     else zuncert=abs(central-down)/central;
187     }
188     datacard << " - \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
189     }
190    
191     datacard << " systematic uncertainty\n";
192    
193    
194     float JESup = ((TH1F*)f->Get("signal_JESUp"))->Integral();
195     float JESdn = ((TH1F*)f->Get("signal_JESDown"))->Integral();
196     float central = signal->Integral();
197     float uJES=0;
198     if(abs(JESup-central)>abs(JESdn-central)) uJES=abs(JESup-central)/central;
199     else uJES=abs(central-JESdn)/central;
200     datacard << "JES lnN";
201     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJES << "\t - \t - \t";
202     datacard << "uncertainty on Jet Energy Scale\n";
203    
204     datacard << "JSU lnN ";
205     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJSU << "\t - \t - \t";
206     datacard << "JZB Scale Uncertainty\n";
207    
208     if(uPDF>0) {
209     datacard << "PDF lnN ";
210     for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uPDF << "\t - \t - \t";
211     datacard << "uncertainty from PDFs\n";
212     }
213    
214     datacard.close();
215     }
216    
217 buchmann 1.2 }
218    
219     float QuickDrawNevents=0;
220    
221 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) {
222     TH1F *histo = new TH1F(hname.c_str(),hname.c_str(),binning.size()-1,&binning[0]);
223     histo->Sumw2();
224     stringstream drawcommand;
225     drawcommand << variable << ">>" << hname;
226     if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
227     events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
228     } else {
229     float totalxs=0;
230     for(int i=0;i<12;i++) {
231     stringstream drawcommand2;
232     drawcommand2 << variable << ">>+" << hname;
233     stringstream mSUGRAweight;
234     float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
235     totalxs+=xs;
236 buchmann 1.4 mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
237 buchmann 1.1 events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
238     }
239 buchmann 1.4 histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
240 buchmann 1.1 }
241     events->Draw("eventNum",addcut.c_str(),"goff");
242     float nevents = events->GetSelectedRows();
243 buchmann 1.2 QuickDrawNevents=nevents;
244 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
245    
246     return histo;
247     }
248    
249    
250     void SQRT(TH1F *h) {
251     for (int i=1;i<=h->GetNbinsX();i++) {
252     h->SetBinContent(i,TMath::Sqrt(h->GetBinContent(i)));
253     }
254     }
255    
256     TH1F* Multiply(TH1F *h1, TH1F *h2) {
257     TH1F *h = (TH1F*)h1->Clone();
258     for(int i=1;i<=h1->GetNbinsX();i++) {
259     h->SetBinContent(i,h1->GetBinContent(i) * h2->GetBinContent(i));
260     }
261     return h;
262     }
263    
264     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) {
265 buchmann 1.4
266 buchmann 1.1 binning.push_back(8000);
267     limfile->cd();
268 buchmann 1.4 dout << "Creating shape templates ";
269 buchmann 1.1 if(identifier!="") dout << "for "<<identifier;
270     dout << " ... " << endl;
271    
272     TCut limitnJetcut;
273 buchmann 1.4
274 buchmann 1.1 if(JES==noJES) limitnJetcut=cutnJets;
275     else {
276     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
277     if(JES==JESup) limitnJetcut=cutnJetsJESup;
278     }
279    
280     float zjetsestimateuncert=zjetsestimateuncertONPEAK;
281     float emuncert=emuncertONPEAK;
282     float emsidebanduncert=emsidebanduncertONPEAK;
283     float eemmsidebanduncert=eemmsidebanduncertONPEAK;
284    
285     if(!PlottingSetup::RestrictToMassPeak) {
286     zjetsestimateuncert=zjetsestimateuncertOFFPEAK;
287     emuncert=emuncertOFFPEAK;
288     emsidebanduncert=emsidebanduncertOFFPEAK;
289     eemmsidebanduncert=eemmsidebanduncertOFFPEAK;
290     }
291    
292     if(signalonly) {
293 buchmann 1.7 dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
294 buchmann 1.1 TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
295     TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
296 buchmann 1.10 TH1F *ZOSOFP;
297     TH1F *ZOSOFN;
298    
299     if(!PlottingSetup::FullMCAnalysis) {
300     ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
301     ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
302     }
303 buchmann 1.1
304     TH1F *SBOSSFP;
305     TH1F *SBOSOFP;
306     TH1F *SBOSSFN;
307     TH1F *SBOSOFN;
308    
309 buchmann 1.10 if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
310 buchmann 1.1 SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
311     SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
312     SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
313     SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
314     }
315    
316     string signalname="signal";
317     if(identifier!="") {
318     signalname="signal_"+identifier;
319     }
320    
321     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
322     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
323    
324     TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
325     TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
326    
327     Lobs->Add(ZOSSFP);
328 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
329 buchmann 1.1
330 buchmann 1.7 dout << "SITUATION FOR SIGNAL: " << endl;
331 buchmann 1.10
332    
333     if(PlottingSetup::FullMCAnalysis) {
334     dout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << endl;
335     } else {
336     dout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << " JZB < 0 :" << ZOSSFN->Integral() << endl;
337     dout << " OSOF JZB> 0 : " << ZOSOFP->Integral() << " JZB < 0 :" << ZOSOFN->Integral() << endl;
338     if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
339     dout << " OSSF SB JZB> 0 : " << SBOSSFP->Integral() << " JZB < 0 :" << SBOSSFN->Integral() << endl;
340     dout << " OSOF SB JZB> 0 : " << SBOSOFP->Integral() << " JZB < 0 :" << SBOSOFN->Integral() << endl;
341     }
342 buchmann 1.2 }
343 buchmann 1.4
344 buchmann 1.2
345 buchmann 1.1 flippedLobs->Add(ZOSSFN);
346 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
347 buchmann 1.1
348 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) {
349     if(PlottingSetup::RestrictToMassPeak) {
350     Lpred->Add(ZOSOFP,1.0/3);
351     Lpred->Add(ZOSOFN,-1.0/3);
352     Lpred->Add(SBOSSFP,1.0/3);
353     Lpred->Add(SBOSSFN,-1.0/3);
354     Lpred->Add(SBOSOFP,1.0/3);
355     Lpred->Add(SBOSOFN,-1.0/3);
356    
357     //flipped prediction
358     flippedLpred->Add(ZOSOFP,-1.0/3);
359     flippedLpred->Add(ZOSOFN,1.0/3);
360     flippedLpred->Add(SBOSSFP,-1.0/3);
361     flippedLpred->Add(SBOSSFN,1.0/3);
362     flippedLpred->Add(SBOSOFP,-1.0/3);
363     flippedLpred->Add(SBOSOFN,1.0/3);
364     } else {
365     Lpred->Add(ZOSOFP,1.0);
366     Lpred->Add(ZOSOFN,-1.0);
367    
368     //flipped prediction
369     flippedLpred->Add(ZOSOFP,-1.0);
370     flippedLpred->Add(ZOSOFN,1.0);
371     }
372 buchmann 1.1 }
373    
374 buchmann 1.2 TH1F *signal = (TH1F*)Lobs->Clone("signal");
375 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
376 buchmann 1.1 signal->SetName(signalname.c_str());
377 buchmann 1.2 signal->SetTitle(signalname.c_str());
378 buchmann 1.1 signal->Write();
379    
380     TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
381 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
382 buchmann 1.1 flippedsignal->SetName(("flipped_"+signalname).c_str());
383     flippedsignal->Write();
384    
385     if(identifier=="") {
386     TH1F *signalStatUp = (TH1F*)signal->Clone();
387     signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
388     signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
389    
390     TH1F *signalStatDn = (TH1F*)signal->Clone();
391     signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
392     signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
393    
394     for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
395 buchmann 1.10 float staterr;
396     if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
397     else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
398    
399     if(!PlottingSetup::FullMCAnalysis) {
400     dout << "Stat err in bin " << i << " : " << staterr << endl;
401     dout << " prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
402     dout << " we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
403     }
404 buchmann 1.5 if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
405     else signalStatDn->SetBinContent(i,0);
406 buchmann 1.2 signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
407     signal->SetBinError(i,staterr);
408 buchmann 1.1 }
409    
410     signalStatDn->Write();
411     signalStatUp->Write();
412    
413     delete signalStatDn;
414     delete signalStatUp;
415    
416     TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
417     flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
418     flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
419    
420     TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
421     flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
422     flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
423    
424     for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
425 buchmann 1.10 float staterr;
426     if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
427     else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
428 buchmann 1.6 if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
429     else flippedsignalStatDn->SetBinContent(i,0);
430 buchmann 1.2 flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
431     flippedsignal->SetBinError(i,staterr);
432 buchmann 1.1 }
433    
434     flippedsignalStatDn->Write();
435     flippedsignalStatUp->Write();
436    
437     delete flippedsignalStatDn;
438     delete flippedsignalStatUp;
439     }
440    
441     delete ZOSSFP;
442     delete ZOSOFP;
443     delete ZOSSFN;
444     delete ZOSOFN;
445    
446     if(PlottingSetup::RestrictToMassPeak) {
447     delete SBOSSFP;
448     delete SBOSOFP;
449     delete SBOSSFN;
450     delete SBOSOFN;
451     }
452    
453     delete Lobs;
454     delete Lpred;
455     delete flippedLobs;
456     delete flippedLpred;
457     }
458    
459    
460     if(dataonly) {
461 buchmann 1.7 dout << "Processing data with datajzb: " << datajzb << endl;
462 buchmann 1.1 TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
463     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
464     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
465     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
466    
467     TH1F *SBOSSFP;
468     TH1F *SBOSOFP;
469     TH1F *SBOSSFN;
470     TH1F *SBOSOFN;
471    
472     if(PlottingSetup::RestrictToMassPeak) {
473     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
474     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
475     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
476     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
477     }
478    
479     string obsname="data_obs";
480     string predname="background";
481     string Zpredname="ZJetsBackground";
482     string Tpredname="TTbarBackground";
483     if(identifier!="") {
484     obsname=("data_"+identifier);
485     predname=("background_"+identifier);
486     Zpredname=("ZJetsBackground_"+identifier);
487     Tpredname=("TTbarBackground_"+identifier);
488     }
489    
490     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
491     obs->SetName(obsname.c_str());
492     obs->Write();
493    
494     TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
495     flippedobs->SetName(("flipped_"+obsname).c_str());
496     flippedobs->Write();
497    
498     TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
499     TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
500    
501     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
502    
503     TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
504     TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
505    
506     TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
507     if(PlottingSetup::RestrictToMassPeak) {
508     pred->Add(ZOSOFP,1.0/3);
509     pred->Add(ZOSOFN,-1.0/3);
510     pred->Add(SBOSSFP,1.0/3);
511     pred->Add(SBOSSFN,-1.0/3);
512     pred->Add(SBOSOFP,1.0/3);
513     pred->Add(SBOSOFN,-1.0/3);
514    
515     Tpred->Add(ZOSSFN,-1.0);
516     Tpred->Add(ZOSOFP,1.0/3);
517     Tpred->Add(SBOSSFP,1.0/3);
518     Tpred->Add(SBOSOFP,1.0/3);
519    
520     Zpred->Add(ZOSOFN,-1.0/3);
521     Zpred->Add(SBOSSFN,-1.0/3);
522     Zpred->Add(SBOSOFN,-1.0/3);
523    
524     //flipped
525     flippedpred->Add(ZOSOFP,-1.0/3);
526     flippedpred->Add(ZOSOFN,1.0/3);
527     flippedpred->Add(SBOSSFP,-1.0/3);
528     flippedpred->Add(SBOSSFN,1.0/3);
529     flippedpred->Add(SBOSOFP,-1.0/3);
530     flippedpred->Add(SBOSOFN,1.0/3);
531    
532     flippedTpred->Add(ZOSSFP,-1.0);
533     flippedTpred->Add(ZOSOFN,1.0/3);
534     flippedTpred->Add(SBOSSFN,1.0/3);
535     flippedTpred->Add(SBOSOFN,1.0/3);
536    
537     flippedZpred->Add(ZOSOFP,-1.0/3);
538     flippedZpred->Add(SBOSSFP,-1.0/3);
539     flippedZpred->Add(SBOSOFP,-1.0/3);
540     } else {
541     pred->Add(ZOSOFP,1.0);
542     pred->Add(ZOSOFN,-1.0);
543    
544     Tpred->Add(ZOSSFN,-1.0);
545     Tpred->Add(ZOSOFP,1.0);
546    
547     Zpred->Add(ZOSOFN,-1.0);
548    
549     //flipped
550     flippedpred->Add(ZOSOFN,1.0);
551     flippedpred->Add(ZOSOFP,-1.0);
552    
553     flippedTpred->Add(ZOSSFP,-1.0);
554     flippedTpred->Add(ZOSOFN,1.0);
555    
556     flippedZpred->Add(ZOSOFP,-1.0);
557     }
558    
559     pred->SetName(predname.c_str());
560     Tpred->SetName(Tpredname.c_str());
561     Zpred->SetName(Zpredname.c_str());
562    
563     flippedpred->SetName(("flipped_"+predname).c_str());
564     flippedTpred->SetName(("flipped_"+Tpredname).c_str());
565     flippedZpred->SetName(("flipped_"+Zpredname).c_str());
566    
567     pred->Write();
568     Tpred->Write();
569     Zpred->Write();
570    
571     flippedpred->Write();
572     flippedTpred->Write();
573     flippedZpred->Write();
574    
575     if(identifier=="") {
576     float stretchfactor=1.0;
577     bool isonpeak=false;
578     if(PlottingSetup::RestrictToMassPeak) {
579     stretchfactor=1.0/3.0;
580     isonpeak=true;
581     }
582    
583     //--------------------------------------------------------statistical uncertainty
584    
585     TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
586     predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
587     predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
588     if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
589     if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
590     if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
591     if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
592     SQRT(predstaterr);
593     TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
594     bgStatUp->Add(predstaterr);
595 buchmann 1.6 EliminateNegativeEntries(bgStatUp);
596 buchmann 1.1 bgStatUp->Write();
597 buchmann 1.2 TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
598 buchmann 1.1 bgStatDn->Add(predstaterr,-1);
599 buchmann 1.6 EliminateNegativeEntries(bgStatDn);
600 buchmann 1.1 bgStatDn->Write();
601     // delete bgStatDn;
602     // delete bgStatUp;
603     delete predstaterr;
604    
605    
606     TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
607     flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
608     flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
609     if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
610     if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
611     if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
612     if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
613     SQRT(flippedpredstaterr);
614     TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
615     fbgStatUp->Add(predstaterr);
616 buchmann 1.6 EliminateNegativeEntries(fbgStatUp);
617 buchmann 1.1 fbgStatUp->Write();
618 buchmann 1.2 TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
619 buchmann 1.1 fbgStatDn->Add(predstaterr,-1);
620 buchmann 1.6 EliminateNegativeEntries(fbgStatDn);
621 buchmann 1.1 fbgStatDn->Write();
622     // delete fbgStatDn;
623     // delete fbgStatUp;
624     delete predstaterr;
625    
626     TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
627     Tpredstaterr->Scale(stretchfactor*stretchfactor);
628     if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
629     if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
630     SQRT(Tpredstaterr);
631     TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
632     TpredStatUp->Add(Tpredstaterr);
633 buchmann 1.6 EliminateNegativeEntries(TpredStatUp);
634 buchmann 1.1 TpredStatUp->Write();
635 buchmann 1.2 TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
636 buchmann 1.1 TpredStatDn->Add(Tpredstaterr,-1);
637 buchmann 1.6 EliminateNegativeEntries(TpredStatDn);
638 buchmann 1.1 TpredStatDn->Write();
639     // delete TpredStatDn;
640     // delete TpredStatUp;
641     delete Tpredstaterr;
642    
643     TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
644     fTpredstaterr->Scale(stretchfactor*stretchfactor);
645     if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
646     if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
647     SQRT(fTpredstaterr);
648     TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
649     fTpredStatUp->Add(fTpredstaterr);
650 buchmann 1.6 EliminateNegativeEntries(fTpredStatUp);
651 buchmann 1.1 fTpredStatUp->Write();
652 buchmann 1.2 TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
653 buchmann 1.1 fTpredStatDn->Add(fTpredstaterr,-1);
654 buchmann 1.6 EliminateNegativeEntries(fTpredStatDn);
655 buchmann 1.1 fTpredStatDn->Write();
656     // delete fTpredStatDn;
657     // delete fTpredStatUp;
658     delete fTpredstaterr;
659    
660    
661     TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
662     Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
663     if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
664     if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
665     SQRT(Zpredstaterr);
666     TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
667     ZpredStatUp->Add(Zpredstaterr);
668 buchmann 1.6 EliminateNegativeEntries(ZpredStatUp);
669 buchmann 1.1 ZpredStatUp->Write();
670 buchmann 1.2 TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
671 buchmann 1.1 ZpredStatDn->Add(Zpredstaterr,-1);
672 buchmann 1.6 EliminateNegativeEntries(ZpredStatDn);
673 buchmann 1.1 ZpredStatDn->Write();
674     // delete ZpredStatDn;
675     // delete ZpredStatUp;
676     delete Zpredstaterr;
677    
678     TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
679     Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
680     if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
681     if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
682     SQRT(fTpredstaterr);
683     TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
684     fZpredStatUp->Add(fZpredstaterr);
685 buchmann 1.6 EliminateNegativeEntries(fZpredStatUp);
686 buchmann 1.1 fZpredStatUp->Write();
687 buchmann 1.2 TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
688 buchmann 1.1 fZpredStatDn->Add(fZpredstaterr,-1);
689 buchmann 1.6 EliminateNegativeEntries(fZpredStatDn);
690 buchmann 1.1 fZpredStatDn->Write();
691     // delete fZpredStatDn;
692     // delete fZpredStatUp;
693     delete fZpredstaterr;
694    
695     //--------------------------------------------------------systematic uncertainty
696    
697     TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
698     predsyserr->Add(pred,-1);//getting everything to zero.
699     predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
700     predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
701     predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
702     if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
703     if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
704     if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
705     if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
706     SQRT(predsyserr);
707     TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
708     bgSysUp->Add(predsyserr);
709 buchmann 1.6 EliminateNegativeEntries(bgSysUp);
710 buchmann 1.1 bgSysUp->Write();
711 buchmann 1.2 TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
712 buchmann 1.1 bgSysDn->Add(predsyserr,-1);
713 buchmann 1.6 EliminateNegativeEntries(bgSysDn);
714 buchmann 1.1 bgSysDn->Write();
715     delete predsyserr;
716    
717     TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
718     fpredsyserr->Add(pred,-1);//getting everything to zero.
719     fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
720     fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
721     fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
722     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
723     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
724     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
725     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
726     SQRT(fpredsyserr);
727     TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
728     fbgSysUp->Add(fpredsyserr);
729 buchmann 1.6 EliminateNegativeEntries(fbgSysUp);
730 buchmann 1.1 fbgSysUp->Write();
731 buchmann 1.2 TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
732 buchmann 1.1 fbgSysDn->Add(fpredsyserr,-1);
733 buchmann 1.6 EliminateNegativeEntries(fbgSysDn);
734 buchmann 1.1 fbgSysDn->Write();
735     delete fpredsyserr;
736    
737    
738     TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
739     Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
740     Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
741     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
742     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
743     SQRT(Tpredsyserr);
744     TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
745     TpredSysUp->Add(Tpredsyserr);
746 buchmann 1.6 EliminateNegativeEntries(TpredSysUp);
747 buchmann 1.1 TpredSysUp->Write();
748 buchmann 1.2 TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
749 buchmann 1.1 TpredSysDn->Add(Tpredsyserr,-1);
750 buchmann 1.6 EliminateNegativeEntries(TpredSysDn);
751 buchmann 1.1 TpredSysDn->Write();
752     delete Tpredsyserr;
753    
754     TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
755     fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
756     fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
757     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
758     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
759     SQRT(fTpredsyserr);
760     TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
761     fTpredSysUp->Add(fTpredsyserr);
762 buchmann 1.6 EliminateNegativeEntries(fTpredSysUp);
763 buchmann 1.1 fTpredSysUp->Write();
764 buchmann 1.2 TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
765 buchmann 1.1 fTpredSysDn->Add(fTpredsyserr,-1);
766 buchmann 1.6 EliminateNegativeEntries(fTpredSysDn);
767 buchmann 1.1 fTpredSysDn->Write();
768     delete fTpredsyserr;
769    
770    
771     //------------------------------------------------------------------------------------------------------------------------
772     TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
773     Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
774     Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
775     Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
776     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
777     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
778     SQRT(Zpredsyserr);
779     TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
780     ZpredSysUp->Add(Zpredsyserr);
781 buchmann 1.6 EliminateNegativeEntries(ZpredSysUp);
782 buchmann 1.1 ZpredSysUp->Write();
783 buchmann 1.2 TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
784 buchmann 1.1 ZpredSysDn->Add(Zpredsyserr,-1);
785 buchmann 1.6 EliminateNegativeEntries(ZpredSysDn);
786 buchmann 1.1 ZpredSysDn->Write();
787     delete Zpredsyserr;
788    
789     TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
790     fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
791     fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
792     fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
793     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
794     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
795     SQRT(fZpredsyserr);
796     TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
797     fZpredSysUp->Add(fZpredsyserr);
798 buchmann 1.6 EliminateNegativeEntries(fZpredSysUp);
799 buchmann 1.1 fZpredSysUp->Write();
800 buchmann 1.2 TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
801 buchmann 1.1 fZpredSysDn->Add(fZpredsyserr,-1);
802 buchmann 1.6 EliminateNegativeEntries(fZpredSysDn);
803 buchmann 1.1 fZpredSysDn->Write();
804     delete fZpredsyserr;
805     }
806    
807 buchmann 1.2 /*if(identifier=="") {
808 buchmann 1.1 for(int i=0;i<binning.size()-1;i++) {
809 buchmann 1.7 dout << "[ " << 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;
810 buchmann 1.1 }
811 buchmann 1.2 }*/
812 buchmann 1.1 delete ZOSSFP;
813     delete ZOSOFP;
814     delete ZOSSFN;
815     delete ZOSOFN;
816    
817     if(PlottingSetup::RestrictToMassPeak) {
818     delete SBOSSFP;
819     delete SBOSOFP;
820     delete SBOSSFN;
821     delete SBOSOFN;
822     }
823     }
824    
825     }
826    
827 buchmann 1.4 ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
828 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
829     if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
830    
831     time_t timestamp;
832     tm *now;
833     timestamp = time(0);
834     now = localtime(&timestamp);
835     stringstream RunDirectoryS;
836     RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
837     string RunDirectory = RunDirectoryS.str();
838    
839     ensure_directory_exists(RunDirectory);
840    
841 buchmann 1.7 TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
842 buchmann 1.1 if(datafile->IsZombie()) {
843     write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
844     assert(!datafile->IsZombie());
845     }
846 buchmann 1.7 dout << "Run Directory: " << RunDirectory << endl;
847 buchmann 1.2 TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
848 buchmann 1.1
849     TIter nextkey(datafile->GetListOfKeys());
850     TKey *key;
851     while ((key = (TKey*)nextkey()))
852     {
853     TH1F *obj =(TH1F*) key->ReadObj();
854     limfile->cd();
855     obj->Write();
856     }
857    
858     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
859    
860     bool signalonly=true;
861     bool dataonly=false;
862    
863     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
864     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
865     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
866    
867     float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
868     float scaledown=0,scaleup=0,scalesyst=0;
869     float mceff=0,mcefferr=0;
870     int flipped=0;
871     int Neventsinfile=10000;//doesn't matter.
872     string informalname="ShapeAnalysis";
873    
874     bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
875    
876     MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
877 buchmann 1.6 if(mceff<0) {
878     flipped=1;
879     write_info(__FUNCTION__,"Doing flipping!");
880     }
881 buchmann 1.1 doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
882 buchmann 1.4 float PDFuncert=0;
883 buchmann 1.1 int NPdfs = get_npdfs(events);
884     if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
885    
886 buchmann 1.2 float JZBscale=scaledown;
887     if(scaleup>scaledown) JZBscale=scaleup;
888     dout << endl;
889     dout << endl;
890     dout << "Will use the following scalar (non-shape) systematics : " << endl;
891     dout << " JZB Scale (JSU) : " << JZBscale << endl;
892     dout << " PDF : " << PDFuncert << endl;
893    
894 buchmann 1.4 float SignalIntegral;
895     if(flipped) {
896     TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
897     SignalIntegral= flisi ->Integral();
898     delete flisi;
899     }
900     else {
901     TH1F *flisi = (TH1F*)(limfile->Get("signal"));
902     SignalIntegral= flisi ->Integral();
903     delete flisi;
904     }
905    
906 buchmann 1.2 TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
907    
908     TIter nnextkey(limfile->GetListOfKeys());
909     TKey *nkey;
910 buchmann 1.4
911     if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
912    
913 buchmann 1.2 while ((nkey = (TKey*)nnextkey()))
914     {
915     TH1F *obj =(TH1F*) nkey->ReadObj();
916     if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
917     if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
918     if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
919     final_limfile->cd();
920 buchmann 1.4 if(Contains(obj->GetName(),"signal")) {
921     obj->Scale(1.0/SignalIntegral);
922     }
923 buchmann 1.2 obj->Write();
924     }
925 buchmann 1.1
926 buchmann 1.6 prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
927 buchmann 1.4
928 buchmann 1.2 final_limfile->Close();
929 buchmann 1.1 limfile->Close();
930 buchmann 1.4 dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
931     stringstream command;
932     if(asymptotic) {
933 buchmann 1.10 if(firstGuess>0) command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
934     else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
935 buchmann 1.4 }
936 buchmann 1.10 else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.333 * firstGuess) << " " << int(firstGuess); // ASYMPTOTIC LIMITS
937 buchmann 1.4 dout <<"Going to run : " << command.str() << endl;
938     int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
939 buchmann 1.7 dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
940 buchmann 1.6 if(!(CreatedModelFileExitCode==0)) {
941     write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
942     ShapeDroplet alpha;
943     alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
944     alpha.SignalIntegral=1;
945     return alpha;
946     }
947 buchmann 1.2 ShapeDroplet alpha;
948     alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
949     alpha.PDF=PDFuncert;
950     alpha.JSU=JZBscale;
951 buchmann 1.4 alpha.SignalIntegral=SignalIntegral;
952 buchmann 1.2 dout << "Done reading limit information from droplet" << endl;
953     dout << alpha << endl;
954 buchmann 1.1
955     dout << "Everything is saved in " << RunDirectory << endl;
956 buchmann 1.7 dout << "Cleaning up ... " << std::flush;
957     /* dout << " 1) Make sure models directory exists ... " << std::flush;
958 buchmann 1.4 gSystem->Exec("mkdir -p models/");
959     dout << " ok!" << endl;
960     dout << " 2) Deleting any previous model files with the same name ... " << std::flush;
961     stringstream delcommand;
962     delcommand << "rm models/model_" << name << ".root";
963     gSystem->Exec(delcommand.str().c_str());
964     stringstream delcommand2;
965     delcommand2 << "rm models/model_" << name << ".txt";
966     gSystem->Exec(delcommand2.str().c_str());
967     dout << " ok!" << endl;
968     dout << " 3) Transfering the new model files ... " << std::flush;
969     stringstream copycommand;
970     copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
971     gSystem->Exec(copycommand.str().c_str());
972     stringstream copycommand2;
973     copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
974     gSystem->Exec(copycommand2.str().c_str());
975     stringstream copycommand3;
976     copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
977     gSystem->Exec(copycommand3.str().c_str());
978     dout << " ok!" << endl;
979 buchmann 1.7 dout << " 4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
980     gSystem->Exec(("rm -r "+RunDirectory).c_str());
981 buchmann 1.4 dout << " ok!" << endl;
982     delete limcan;
983     return alpha;
984 buchmann 1.1 }
985    
986     void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
987 buchmann 1.2 // dout << "Generating data shape templates ... " << endl;
988 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
989     TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
990     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
991     TTree *faketree;
992     bool dataonly=true;
993     bool signalonly=false;
994     string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
995 buchmann 1.8 // float jzbpeakerrormc=0;
996 buchmann 1.1 generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
997     // don't need these effects for obs & pred, only for signal!
998     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
999     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1000     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
1001     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
1002     datafile->Close();
1003     }
1004