ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.4
Committed: Wed Apr 25 12:59:53 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.3: +72 -25 lines
Log Message: 
Various changes in shape analysis: not normalizing empty histos, more verbose, implemented mSUGRA histos, added first guessing for asymptotic limits, added second full run if necessary with arguments, and made the datacard more legible

File Contents

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