ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.9
Committed: Fri May 4 11:47:57 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.8: +159 -30 lines
Log Message: 
adapted card

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 *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
296     TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
297     TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
298    
299     TH1F *SBOSSFP;
300     TH1F *SBOSOFP;
301     TH1F *SBOSSFN;
302     TH1F *SBOSOFN;
303    
304     if(PlottingSetup::RestrictToMassPeak) {
305     SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
306     SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
307     SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
308     SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
309     }
310    
311     string signalname="signal";
312     if(identifier!="") {
313     signalname="signal_"+identifier;
314     }
315    
316     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
317     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
318    
319     TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
320     TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
321    
322     Lobs->Add(ZOSSFP);
323     Lpred->Add(ZOSSFN);
324    
325 buchmann 1.7 dout << "SITUATION FOR SIGNAL: " << endl;
326     dout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << " JZB < 0 :" << ZOSSFN->Integral() << endl;
327     dout << " OSOF JZB> 0 : " << ZOSOFP->Integral() << " JZB < 0 :" << ZOSOFN->Integral() << endl;
328 buchmann 1.2 if(PlottingSetup::RestrictToMassPeak) {
329 buchmann 1.7 dout << " OSSF SB JZB> 0 : " << SBOSSFP->Integral() << " JZB < 0 :" << SBOSSFN->Integral() << endl;
330     dout << " OSOF SB JZB> 0 : " << SBOSOFP->Integral() << " JZB < 0 :" << SBOSOFN->Integral() << endl;
331 buchmann 1.2 }
332 buchmann 1.4
333 buchmann 1.2
334 buchmann 1.1 flippedLobs->Add(ZOSSFN);
335     flippedLpred->Add(ZOSSFP);
336    
337     if(PlottingSetup::RestrictToMassPeak) {
338     Lpred->Add(ZOSOFP,1.0/3);
339     Lpred->Add(ZOSOFN,-1.0/3);
340     Lpred->Add(SBOSSFP,1.0/3);
341     Lpred->Add(SBOSSFN,-1.0/3);
342     Lpred->Add(SBOSOFP,1.0/3);
343     Lpred->Add(SBOSOFN,-1.0/3);
344    
345     //flipped prediction
346     flippedLpred->Add(ZOSOFP,-1.0/3);
347     flippedLpred->Add(ZOSOFN,1.0/3);
348     flippedLpred->Add(SBOSSFP,-1.0/3);
349     flippedLpred->Add(SBOSSFN,1.0/3);
350     flippedLpred->Add(SBOSOFP,-1.0/3);
351     flippedLpred->Add(SBOSOFN,1.0/3);
352    
353     } else {
354     Lpred->Add(ZOSOFP,1.0);
355     Lpred->Add(ZOSOFN,-1.0);
356 buchmann 1.2
357 buchmann 1.1 //flipped prediction
358     flippedLpred->Add(ZOSOFP,-1.0);
359     flippedLpred->Add(ZOSOFN,1.0);
360     }
361    
362 buchmann 1.2 TH1F *signal = (TH1F*)Lobs->Clone("signal");
363 buchmann 1.1 signal->Add(Lpred,-1);
364     signal->SetName(signalname.c_str());
365 buchmann 1.2 signal->SetTitle(signalname.c_str());
366 buchmann 1.1 signal->Write();
367    
368     TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
369     flippedsignal->Add(flippedLpred,-1);
370     flippedsignal->SetName(("flipped_"+signalname).c_str());
371     flippedsignal->Write();
372    
373     if(identifier=="") {
374     TH1F *signalStatUp = (TH1F*)signal->Clone();
375     signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
376     signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
377    
378     TH1F *signalStatDn = (TH1F*)signal->Clone();
379     signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
380     signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
381    
382     for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
383 buchmann 1.2 float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
384 buchmann 1.7 dout << "Stat err in bin " << i << " : " << staterr << endl;
385     dout << " prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
386     dout << " we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
387 buchmann 1.5 if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
388     else signalStatDn->SetBinContent(i,0);
389 buchmann 1.2 signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
390     signal->SetBinError(i,staterr);
391 buchmann 1.1 }
392    
393     signalStatDn->Write();
394     signalStatUp->Write();
395    
396     delete signalStatDn;
397     delete signalStatUp;
398    
399     TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
400     flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
401     flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
402    
403     TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
404     flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
405     flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
406    
407     for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
408 buchmann 1.2 float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
409 buchmann 1.6 if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
410     else flippedsignalStatDn->SetBinContent(i,0);
411 buchmann 1.2 flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
412     flippedsignal->SetBinError(i,staterr);
413 buchmann 1.1 }
414    
415     flippedsignalStatDn->Write();
416     flippedsignalStatUp->Write();
417    
418     delete flippedsignalStatDn;
419     delete flippedsignalStatUp;
420     }
421    
422     delete ZOSSFP;
423     delete ZOSOFP;
424     delete ZOSSFN;
425     delete ZOSOFN;
426    
427     if(PlottingSetup::RestrictToMassPeak) {
428     delete SBOSSFP;
429     delete SBOSOFP;
430     delete SBOSSFN;
431     delete SBOSOFN;
432     }
433    
434     delete Lobs;
435     delete Lpred;
436     delete flippedLobs;
437     delete flippedLpred;
438     }
439    
440    
441     if(dataonly) {
442 buchmann 1.7 dout << "Processing data with datajzb: " << datajzb << endl;
443 buchmann 1.1 TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
444     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
445     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
446     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
447    
448     TH1F *SBOSSFP;
449     TH1F *SBOSOFP;
450     TH1F *SBOSSFN;
451     TH1F *SBOSOFN;
452    
453     if(PlottingSetup::RestrictToMassPeak) {
454     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
455     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
456     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
457     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
458     }
459    
460     string obsname="data_obs";
461     string predname="background";
462     string Zpredname="ZJetsBackground";
463     string Tpredname="TTbarBackground";
464     if(identifier!="") {
465     obsname=("data_"+identifier);
466     predname=("background_"+identifier);
467     Zpredname=("ZJetsBackground_"+identifier);
468     Tpredname=("TTbarBackground_"+identifier);
469     }
470    
471     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
472     obs->SetName(obsname.c_str());
473     obs->Write();
474    
475     TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
476     flippedobs->SetName(("flipped_"+obsname).c_str());
477     flippedobs->Write();
478    
479     TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
480     TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
481    
482     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
483    
484     TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
485     TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
486    
487     TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
488     if(PlottingSetup::RestrictToMassPeak) {
489     pred->Add(ZOSOFP,1.0/3);
490     pred->Add(ZOSOFN,-1.0/3);
491     pred->Add(SBOSSFP,1.0/3);
492     pred->Add(SBOSSFN,-1.0/3);
493     pred->Add(SBOSOFP,1.0/3);
494     pred->Add(SBOSOFN,-1.0/3);
495    
496     Tpred->Add(ZOSSFN,-1.0);
497     Tpred->Add(ZOSOFP,1.0/3);
498     Tpred->Add(SBOSSFP,1.0/3);
499     Tpred->Add(SBOSOFP,1.0/3);
500    
501     Zpred->Add(ZOSOFN,-1.0/3);
502     Zpred->Add(SBOSSFN,-1.0/3);
503     Zpred->Add(SBOSOFN,-1.0/3);
504    
505     //flipped
506     flippedpred->Add(ZOSOFP,-1.0/3);
507     flippedpred->Add(ZOSOFN,1.0/3);
508     flippedpred->Add(SBOSSFP,-1.0/3);
509     flippedpred->Add(SBOSSFN,1.0/3);
510     flippedpred->Add(SBOSOFP,-1.0/3);
511     flippedpred->Add(SBOSOFN,1.0/3);
512    
513     flippedTpred->Add(ZOSSFP,-1.0);
514     flippedTpred->Add(ZOSOFN,1.0/3);
515     flippedTpred->Add(SBOSSFN,1.0/3);
516     flippedTpred->Add(SBOSOFN,1.0/3);
517    
518     flippedZpred->Add(ZOSOFP,-1.0/3);
519     flippedZpred->Add(SBOSSFP,-1.0/3);
520     flippedZpred->Add(SBOSOFP,-1.0/3);
521     } else {
522     pred->Add(ZOSOFP,1.0);
523     pred->Add(ZOSOFN,-1.0);
524    
525     Tpred->Add(ZOSSFN,-1.0);
526     Tpred->Add(ZOSOFP,1.0);
527    
528     Zpred->Add(ZOSOFN,-1.0);
529    
530     //flipped
531     flippedpred->Add(ZOSOFN,1.0);
532     flippedpred->Add(ZOSOFP,-1.0);
533    
534     flippedTpred->Add(ZOSSFP,-1.0);
535     flippedTpred->Add(ZOSOFN,1.0);
536    
537     flippedZpred->Add(ZOSOFP,-1.0);
538     }
539    
540     pred->SetName(predname.c_str());
541     Tpred->SetName(Tpredname.c_str());
542     Zpred->SetName(Zpredname.c_str());
543    
544     flippedpred->SetName(("flipped_"+predname).c_str());
545     flippedTpred->SetName(("flipped_"+Tpredname).c_str());
546     flippedZpred->SetName(("flipped_"+Zpredname).c_str());
547    
548     pred->Write();
549     Tpred->Write();
550     Zpred->Write();
551    
552     flippedpred->Write();
553     flippedTpred->Write();
554     flippedZpred->Write();
555    
556     if(identifier=="") {
557     float stretchfactor=1.0;
558     bool isonpeak=false;
559     if(PlottingSetup::RestrictToMassPeak) {
560     stretchfactor=1.0/3.0;
561     isonpeak=true;
562     }
563    
564     //--------------------------------------------------------statistical uncertainty
565    
566     TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
567     predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
568     predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
569     if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
570     if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
571     if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
572     if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
573     SQRT(predstaterr);
574     TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
575     bgStatUp->Add(predstaterr);
576 buchmann 1.6 EliminateNegativeEntries(bgStatUp);
577 buchmann 1.1 bgStatUp->Write();
578 buchmann 1.2 TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
579 buchmann 1.1 bgStatDn->Add(predstaterr,-1);
580 buchmann 1.6 EliminateNegativeEntries(bgStatDn);
581 buchmann 1.1 bgStatDn->Write();
582     // delete bgStatDn;
583     // delete bgStatUp;
584     delete predstaterr;
585    
586    
587     TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
588     flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
589     flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
590     if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
591     if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
592     if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
593     if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
594     SQRT(flippedpredstaterr);
595     TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
596     fbgStatUp->Add(predstaterr);
597 buchmann 1.6 EliminateNegativeEntries(fbgStatUp);
598 buchmann 1.1 fbgStatUp->Write();
599 buchmann 1.2 TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
600 buchmann 1.1 fbgStatDn->Add(predstaterr,-1);
601 buchmann 1.6 EliminateNegativeEntries(fbgStatDn);
602 buchmann 1.1 fbgStatDn->Write();
603     // delete fbgStatDn;
604     // delete fbgStatUp;
605     delete predstaterr;
606    
607     TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
608     Tpredstaterr->Scale(stretchfactor*stretchfactor);
609     if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
610     if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
611     SQRT(Tpredstaterr);
612     TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
613     TpredStatUp->Add(Tpredstaterr);
614 buchmann 1.6 EliminateNegativeEntries(TpredStatUp);
615 buchmann 1.1 TpredStatUp->Write();
616 buchmann 1.2 TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
617 buchmann 1.1 TpredStatDn->Add(Tpredstaterr,-1);
618 buchmann 1.6 EliminateNegativeEntries(TpredStatDn);
619 buchmann 1.1 TpredStatDn->Write();
620     // delete TpredStatDn;
621     // delete TpredStatUp;
622     delete Tpredstaterr;
623    
624     TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
625     fTpredstaterr->Scale(stretchfactor*stretchfactor);
626     if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
627     if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
628     SQRT(fTpredstaterr);
629     TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
630     fTpredStatUp->Add(fTpredstaterr);
631 buchmann 1.6 EliminateNegativeEntries(fTpredStatUp);
632 buchmann 1.1 fTpredStatUp->Write();
633 buchmann 1.2 TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
634 buchmann 1.1 fTpredStatDn->Add(fTpredstaterr,-1);
635 buchmann 1.6 EliminateNegativeEntries(fTpredStatDn);
636 buchmann 1.1 fTpredStatDn->Write();
637     // delete fTpredStatDn;
638     // delete fTpredStatUp;
639     delete fTpredstaterr;
640    
641    
642     TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
643     Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
644     if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
645     if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
646     SQRT(Zpredstaterr);
647     TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
648     ZpredStatUp->Add(Zpredstaterr);
649 buchmann 1.6 EliminateNegativeEntries(ZpredStatUp);
650 buchmann 1.1 ZpredStatUp->Write();
651 buchmann 1.2 TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
652 buchmann 1.1 ZpredStatDn->Add(Zpredstaterr,-1);
653 buchmann 1.6 EliminateNegativeEntries(ZpredStatDn);
654 buchmann 1.1 ZpredStatDn->Write();
655     // delete ZpredStatDn;
656     // delete ZpredStatUp;
657     delete Zpredstaterr;
658    
659     TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
660     Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
661     if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
662     if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
663     SQRT(fTpredstaterr);
664     TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
665     fZpredStatUp->Add(fZpredstaterr);
666 buchmann 1.6 EliminateNegativeEntries(fZpredStatUp);
667 buchmann 1.1 fZpredStatUp->Write();
668 buchmann 1.2 TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
669 buchmann 1.1 fZpredStatDn->Add(fZpredstaterr,-1);
670 buchmann 1.6 EliminateNegativeEntries(fZpredStatDn);
671 buchmann 1.1 fZpredStatDn->Write();
672     // delete fZpredStatDn;
673     // delete fZpredStatUp;
674     delete fZpredstaterr;
675    
676     //--------------------------------------------------------systematic uncertainty
677    
678     TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
679     predsyserr->Add(pred,-1);//getting everything to zero.
680     predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
681     predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
682     predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
683     if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
684     if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
685     if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
686     if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
687     SQRT(predsyserr);
688     TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
689     bgSysUp->Add(predsyserr);
690 buchmann 1.6 EliminateNegativeEntries(bgSysUp);
691 buchmann 1.1 bgSysUp->Write();
692 buchmann 1.2 TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
693 buchmann 1.1 bgSysDn->Add(predsyserr,-1);
694 buchmann 1.6 EliminateNegativeEntries(bgSysDn);
695 buchmann 1.1 bgSysDn->Write();
696     delete predsyserr;
697    
698     TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
699     fpredsyserr->Add(pred,-1);//getting everything to zero.
700     fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
701     fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
702     fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
703     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
704     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
705     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
706     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
707     SQRT(fpredsyserr);
708     TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
709     fbgSysUp->Add(fpredsyserr);
710 buchmann 1.6 EliminateNegativeEntries(fbgSysUp);
711 buchmann 1.1 fbgSysUp->Write();
712 buchmann 1.2 TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
713 buchmann 1.1 fbgSysDn->Add(fpredsyserr,-1);
714 buchmann 1.6 EliminateNegativeEntries(fbgSysDn);
715 buchmann 1.1 fbgSysDn->Write();
716     delete fpredsyserr;
717    
718    
719     TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
720     Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
721     Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
722     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
723     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
724     SQRT(Tpredsyserr);
725     TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
726     TpredSysUp->Add(Tpredsyserr);
727 buchmann 1.6 EliminateNegativeEntries(TpredSysUp);
728 buchmann 1.1 TpredSysUp->Write();
729 buchmann 1.2 TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
730 buchmann 1.1 TpredSysDn->Add(Tpredsyserr,-1);
731 buchmann 1.6 EliminateNegativeEntries(TpredSysDn);
732 buchmann 1.1 TpredSysDn->Write();
733     delete Tpredsyserr;
734    
735     TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
736     fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
737     fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
738     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
739     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
740     SQRT(fTpredsyserr);
741     TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
742     fTpredSysUp->Add(fTpredsyserr);
743 buchmann 1.6 EliminateNegativeEntries(fTpredSysUp);
744 buchmann 1.1 fTpredSysUp->Write();
745 buchmann 1.2 TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
746 buchmann 1.1 fTpredSysDn->Add(fTpredsyserr,-1);
747 buchmann 1.6 EliminateNegativeEntries(fTpredSysDn);
748 buchmann 1.1 fTpredSysDn->Write();
749     delete fTpredsyserr;
750    
751    
752     //------------------------------------------------------------------------------------------------------------------------
753     TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
754     Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
755     Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
756     Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
757     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
758     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
759     SQRT(Zpredsyserr);
760     TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
761     ZpredSysUp->Add(Zpredsyserr);
762 buchmann 1.6 EliminateNegativeEntries(ZpredSysUp);
763 buchmann 1.1 ZpredSysUp->Write();
764 buchmann 1.2 TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
765 buchmann 1.1 ZpredSysDn->Add(Zpredsyserr,-1);
766 buchmann 1.6 EliminateNegativeEntries(ZpredSysDn);
767 buchmann 1.1 ZpredSysDn->Write();
768     delete Zpredsyserr;
769    
770     TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
771     fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
772     fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
773     fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
774     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
775     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
776     SQRT(fZpredsyserr);
777     TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
778     fZpredSysUp->Add(fZpredsyserr);
779 buchmann 1.6 EliminateNegativeEntries(fZpredSysUp);
780 buchmann 1.1 fZpredSysUp->Write();
781 buchmann 1.2 TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
782 buchmann 1.1 fZpredSysDn->Add(fZpredsyserr,-1);
783 buchmann 1.6 EliminateNegativeEntries(fZpredSysDn);
784 buchmann 1.1 fZpredSysDn->Write();
785     delete fZpredsyserr;
786     }
787    
788 buchmann 1.2 /*if(identifier=="") {
789 buchmann 1.1 for(int i=0;i<binning.size()-1;i++) {
790 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;
791 buchmann 1.1 }
792 buchmann 1.2 }*/
793 buchmann 1.1 delete ZOSSFP;
794     delete ZOSOFP;
795     delete ZOSSFN;
796     delete ZOSOFN;
797    
798     if(PlottingSetup::RestrictToMassPeak) {
799     delete SBOSSFP;
800     delete SBOSOFP;
801     delete SBOSSFN;
802     delete SBOSOFN;
803     }
804     }
805    
806     }
807    
808 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) {
809 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
810     if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
811    
812     time_t timestamp;
813     tm *now;
814     timestamp = time(0);
815     now = localtime(&timestamp);
816     stringstream RunDirectoryS;
817     RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
818     string RunDirectory = RunDirectoryS.str();
819    
820     ensure_directory_exists(RunDirectory);
821    
822 buchmann 1.7 TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
823 buchmann 1.1 if(datafile->IsZombie()) {
824     write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
825     assert(!datafile->IsZombie());
826     }
827 buchmann 1.7 dout << "Run Directory: " << RunDirectory << endl;
828 buchmann 1.2 TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
829 buchmann 1.1
830     TIter nextkey(datafile->GetListOfKeys());
831     TKey *key;
832     while ((key = (TKey*)nextkey()))
833     {
834     TH1F *obj =(TH1F*) key->ReadObj();
835     limfile->cd();
836     obj->Write();
837     }
838    
839     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
840    
841     bool signalonly=true;
842     bool dataonly=false;
843    
844     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
845 buchmann 1.2 // generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
846     // generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
847 buchmann 1.1 generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
848     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
849    
850     float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
851     float scaledown=0,scaleup=0,scalesyst=0;
852     float mceff=0,mcefferr=0;
853     int flipped=0;
854     int Neventsinfile=10000;//doesn't matter.
855     string informalname="ShapeAnalysis";
856    
857     bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
858    
859     MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
860 buchmann 1.6 if(mceff<0) {
861     flipped=1;
862     write_info(__FUNCTION__,"Doing flipping!");
863     }
864 buchmann 1.1 doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
865 buchmann 1.4 float PDFuncert=0;
866 buchmann 1.1 int NPdfs = get_npdfs(events);
867     if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
868    
869 buchmann 1.2 float JZBscale=scaledown;
870     if(scaleup>scaledown) JZBscale=scaleup;
871     dout << endl;
872     dout << endl;
873     dout << "Will use the following scalar (non-shape) systematics : " << endl;
874     dout << " JZB Scale (JSU) : " << JZBscale << endl;
875     dout << " PDF : " << PDFuncert << endl;
876    
877 buchmann 1.4 float SignalIntegral;
878     if(flipped) {
879     TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
880     SignalIntegral= flisi ->Integral();
881     delete flisi;
882     }
883     else {
884     TH1F *flisi = (TH1F*)(limfile->Get("signal"));
885     SignalIntegral= flisi ->Integral();
886     delete flisi;
887     }
888    
889 buchmann 1.2 TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
890    
891     TIter nnextkey(limfile->GetListOfKeys());
892     TKey *nkey;
893 buchmann 1.4
894     if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
895    
896 buchmann 1.2 while ((nkey = (TKey*)nnextkey()))
897     {
898     TH1F *obj =(TH1F*) nkey->ReadObj();
899     if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
900     if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
901     if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
902     final_limfile->cd();
903 buchmann 1.4 if(Contains(obj->GetName(),"signal")) {
904     obj->Scale(1.0/SignalIntegral);
905     }
906 buchmann 1.2 obj->Write();
907     }
908 buchmann 1.1
909 buchmann 1.6 prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
910 buchmann 1.4
911 buchmann 1.2 final_limfile->Close();
912 buchmann 1.1 limfile->Close();
913 buchmann 1.4 dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
914     stringstream command;
915     if(asymptotic) {
916     if(firstGuess>0) command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
917     else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
918     }
919 buchmann 1.9 else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.333 * firstGuess) << " " << int(firstGuess); // ASYMPTOTIC LIMITS
920 buchmann 1.4 dout <<"Going to run : " << command.str() << endl;
921     int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
922 buchmann 1.7 dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
923 buchmann 1.6 if(!(CreatedModelFileExitCode==0)) {
924     write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
925     ShapeDroplet alpha;
926     alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
927     alpha.SignalIntegral=1;
928     return alpha;
929     }
930 buchmann 1.2 ShapeDroplet alpha;
931     alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
932     alpha.PDF=PDFuncert;
933     alpha.JSU=JZBscale;
934 buchmann 1.4 alpha.SignalIntegral=SignalIntegral;
935 buchmann 1.2 dout << "Done reading limit information from droplet" << endl;
936     dout << alpha << endl;
937 buchmann 1.1
938     dout << "Everything is saved in " << RunDirectory << endl;
939 buchmann 1.7 dout << "Cleaning up ... " << std::flush;
940     /* dout << " 1) Make sure models directory exists ... " << std::flush;
941 buchmann 1.4 gSystem->Exec("mkdir -p models/");
942     dout << " ok!" << endl;
943     dout << " 2) Deleting any previous model files with the same name ... " << std::flush;
944     stringstream delcommand;
945     delcommand << "rm models/model_" << name << ".root";
946     gSystem->Exec(delcommand.str().c_str());
947     stringstream delcommand2;
948     delcommand2 << "rm models/model_" << name << ".txt";
949     gSystem->Exec(delcommand2.str().c_str());
950     dout << " ok!" << endl;
951     dout << " 3) Transfering the new model files ... " << std::flush;
952     stringstream copycommand;
953     copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
954     gSystem->Exec(copycommand.str().c_str());
955     stringstream copycommand2;
956     copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
957     gSystem->Exec(copycommand2.str().c_str());
958     stringstream copycommand3;
959     copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
960     gSystem->Exec(copycommand3.str().c_str());
961     dout << " ok!" << endl;
962 buchmann 1.7 dout << " 4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
963     gSystem->Exec(("rm -r "+RunDirectory).c_str());
964 buchmann 1.4 dout << " ok!" << endl;
965     delete limcan;
966     return alpha;
967 buchmann 1.1 }
968    
969     void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
970 buchmann 1.2 // dout << "Generating data shape templates ... " << endl;
971 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
972     TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
973     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
974     TTree *faketree;
975     bool dataonly=true;
976     bool signalonly=false;
977     string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
978 buchmann 1.8 // float jzbpeakerrormc=0;
979 buchmann 1.1 generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
980     // don't need these effects for obs & pred, only for signal!
981     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
982     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
983     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
984     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
985     datafile->Close();
986     }
987