ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.6
Committed: Fri Apr 27 08:44:49 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.5: +53 -24 lines
Log Message:
Catching negative entries and eliminating them; using the final shape file directly

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