ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.1
Committed: Wed Apr 11 15:37:24 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Log Message:
Added a very early version of the shape module; the data part works and has been verified for all scenarios (fp,np, for regular and flipped) so don't mess with it :-)

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     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) {
36     TH1F *histo = new TH1F(hname.c_str(),hname.c_str(),binning.size()-1,&binning[0]);
37     histo->Sumw2();
38     stringstream drawcommand;
39     drawcommand << variable << ">>" << hname;
40     if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
41     events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
42     } else {
43     float totalxs=0;
44     for(int i=0;i<12;i++) {
45     stringstream drawcommand2;
46     drawcommand2 << variable << ">>+" << hname;
47     stringstream mSUGRAweight;
48     float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
49     totalxs+=xs;
50     events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
51     }
52     histo->Scale(1.0/totalxs);
53     }
54     events->Draw("eventNum",addcut.c_str(),"goff");
55     float nevents = events->GetSelectedRows();
56     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
57    
58     return histo;
59     }
60    
61     /*void do_stat_up(TH1F *h, float sign=1.0) {
62     for(int i=1;i<=h->GetNbinsX();i++) {
63     h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
64     }
65     h->Write();
66     }
67    
68     void do_stat_dn(TH1F *h) {
69     do_stat_up(h,-1.0);
70     }*/
71    
72     void SQRT(TH1F *h) {
73     for (int i=1;i<=h->GetNbinsX();i++) {
74     h->SetBinContent(i,TMath::Sqrt(h->GetBinContent(i)));
75     }
76     }
77    
78     TH1F* Multiply(TH1F *h1, TH1F *h2) {
79     TH1F *h = (TH1F*)h1->Clone();
80     for(int i=1;i<=h1->GetNbinsX();i++) {
81     h->SetBinContent(i,h1->GetBinContent(i) * h2->GetBinContent(i));
82     }
83     return h;
84     }
85    
86     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) {
87     binning.push_back(8000);
88     limfile->cd();
89     dout << "Creatig shape templates ";
90     if(identifier!="") dout << "for "<<identifier;
91     dout << " ... " << endl;
92    
93     TCut limitnJetcut;
94     if(JES==noJES) limitnJetcut=cutnJets;
95     else {
96     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
97     if(JES==JESup) limitnJetcut=cutnJetsJESup;
98     }
99    
100     stringstream mSUGRAweight;
101     if(SUSYScanSpace::SUSYscantype==mSUGRA) {
102     //for(int i<
103     mSUGRAweight << "bla";
104     } else mSUGRAweight << "(1)";
105    
106     write_warning(__FUNCTION__,"Not done yet for mSUGRA");
107    
108     float zjetsestimateuncert=zjetsestimateuncertONPEAK;
109     float emuncert=emuncertONPEAK;
110     float emsidebanduncert=emsidebanduncertONPEAK;
111     float eemmsidebanduncert=eemmsidebanduncertONPEAK;
112    
113     if(!PlottingSetup::RestrictToMassPeak) {
114     zjetsestimateuncert=zjetsestimateuncertOFFPEAK;
115     emuncert=emuncertOFFPEAK;
116     emsidebanduncert=emsidebanduncertOFFPEAK;
117     eemmsidebanduncert=eemmsidebanduncertOFFPEAK;
118     }
119    
120     if(signalonly) {
121     cout << "Processing a signal with mcjzb: " << mcjzb << endl;
122     TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
123     TH1F *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
124     TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
125     TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
126    
127     TH1F *SBOSSFP;
128     TH1F *SBOSOFP;
129     TH1F *SBOSSFN;
130     TH1F *SBOSOFN;
131    
132     if(PlottingSetup::RestrictToMassPeak) {
133     SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
134     SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
135     SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
136     SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
137     }
138    
139     string signalname="signal";
140     if(identifier!="") {
141     signalname="signal_"+identifier;
142     }
143    
144     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
145     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
146    
147     TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
148     TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
149    
150     Lobs->Add(ZOSSFP);
151     Lpred->Add(ZOSSFN);
152    
153     flippedLobs->Add(ZOSSFN);
154     flippedLpred->Add(ZOSSFP);
155    
156     if(PlottingSetup::RestrictToMassPeak) {
157     Lpred->Add(ZOSOFP,1.0/3);
158     Lpred->Add(ZOSOFN,-1.0/3);
159     Lpred->Add(SBOSSFP,1.0/3);
160     Lpred->Add(SBOSSFN,-1.0/3);
161     Lpred->Add(SBOSOFP,1.0/3);
162     Lpred->Add(SBOSOFN,-1.0/3);
163    
164     //flipped prediction
165     flippedLpred->Add(ZOSOFP,-1.0/3);
166     flippedLpred->Add(ZOSOFN,1.0/3);
167     flippedLpred->Add(SBOSSFP,-1.0/3);
168     flippedLpred->Add(SBOSSFN,1.0/3);
169     flippedLpred->Add(SBOSOFP,-1.0/3);
170     flippedLpred->Add(SBOSOFN,1.0/3);
171    
172     } else {
173     Lpred->Add(ZOSOFP,1.0);
174     Lpred->Add(ZOSOFN,-1.0);
175    
176     //flipped prediction
177     flippedLpred->Add(ZOSOFP,-1.0);
178     flippedLpred->Add(ZOSOFN,1.0);
179     }
180    
181     TH1F *signal = (TH1F*)Lobs->Clone();
182     signal->Add(Lpred,-1);
183     signal->SetName(signalname.c_str());
184     signal->Write();
185    
186     TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
187     flippedsignal->Add(flippedLpred,-1);
188     flippedsignal->SetName(("flipped_"+signalname).c_str());
189     flippedsignal->Write();
190    
191     if(identifier=="") {
192     write_warning(__FUNCTION__,"Don't know right now how to define stat err on signal");
193     TH1F *signalStatUp = (TH1F*)signal->Clone();
194     signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
195     signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
196    
197     TH1F *signalStatDn = (TH1F*)signal->Clone();
198     signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
199     signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
200    
201     for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
202     signalStatDn->SetBinContent(i,signal->GetBinContent(i)-signal->GetBinError(i));
203     signalStatUp->SetBinContent(i,signal->GetBinContent(i)+signal->GetBinError(i));
204     }
205    
206     signalStatDn->Write();
207     signalStatUp->Write();
208    
209     delete signalStatDn;
210     delete signalStatUp;
211    
212     TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
213     flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
214     flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
215    
216     TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
217     flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
218     flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
219    
220     for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
221     flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-flippedsignal->GetBinError(i));
222     flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+flippedsignal->GetBinError(i));
223     }
224    
225     flippedsignalStatDn->Write();
226     flippedsignalStatUp->Write();
227    
228     delete flippedsignalStatDn;
229     delete flippedsignalStatUp;
230     }
231    
232     delete ZOSSFP;
233     delete ZOSOFP;
234     delete ZOSSFN;
235     delete ZOSOFN;
236    
237     if(PlottingSetup::RestrictToMassPeak) {
238     delete SBOSSFP;
239     delete SBOSOFP;
240     delete SBOSSFN;
241     delete SBOSOFN;
242     }
243    
244     delete Lobs;
245     delete Lpred;
246     delete flippedLobs;
247     delete flippedLpred;
248     }
249    
250    
251     if(dataonly) {
252     cout << "Processing data with datajzb: " << datajzb << endl;
253     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
254     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
255     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
256     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
257    
258     TH1F *SBOSSFP;
259     TH1F *SBOSOFP;
260     TH1F *SBOSSFN;
261     TH1F *SBOSOFN;
262    
263     if(PlottingSetup::RestrictToMassPeak) {
264     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
265     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
266     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
267     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
268     }
269    
270     string obsname="data_obs";
271     string predname="background";
272     string Zpredname="ZJetsBackground";
273     string Tpredname="TTbarBackground";
274     if(identifier!="") {
275     obsname=("data_"+identifier);
276     predname=("background_"+identifier);
277     Zpredname=("ZJetsBackground_"+identifier);
278     Tpredname=("TTbarBackground_"+identifier);
279     }
280    
281     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
282     obs->SetName(obsname.c_str());
283     obs->Write();
284    
285     TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
286     flippedobs->SetName(("flipped_"+obsname).c_str());
287     flippedobs->Write();
288    
289     TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
290     TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
291    
292     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
293    
294     TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
295     TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
296    
297     TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
298     if(PlottingSetup::RestrictToMassPeak) {
299     pred->Add(ZOSOFP,1.0/3);
300     pred->Add(ZOSOFN,-1.0/3);
301     pred->Add(SBOSSFP,1.0/3);
302     pred->Add(SBOSSFN,-1.0/3);
303     pred->Add(SBOSOFP,1.0/3);
304     pred->Add(SBOSOFN,-1.0/3);
305    
306     Tpred->Add(ZOSSFN,-1.0);
307     Tpred->Add(ZOSOFP,1.0/3);
308     Tpred->Add(SBOSSFP,1.0/3);
309     Tpred->Add(SBOSOFP,1.0/3);
310    
311     Zpred->Add(ZOSOFN,-1.0/3);
312     Zpred->Add(SBOSSFN,-1.0/3);
313     Zpred->Add(SBOSOFN,-1.0/3);
314    
315     //flipped
316     flippedpred->Add(ZOSOFP,-1.0/3);
317     flippedpred->Add(ZOSOFN,1.0/3);
318     flippedpred->Add(SBOSSFP,-1.0/3);
319     flippedpred->Add(SBOSSFN,1.0/3);
320     flippedpred->Add(SBOSOFP,-1.0/3);
321     flippedpred->Add(SBOSOFN,1.0/3);
322    
323     flippedTpred->Add(ZOSSFP,-1.0);
324     flippedTpred->Add(ZOSOFN,1.0/3);
325     flippedTpred->Add(SBOSSFN,1.0/3);
326     flippedTpred->Add(SBOSOFN,1.0/3);
327    
328     flippedZpred->Add(ZOSOFP,-1.0/3);
329     flippedZpred->Add(SBOSSFP,-1.0/3);
330     flippedZpred->Add(SBOSOFP,-1.0/3);
331     } else {
332     pred->Add(ZOSOFP,1.0);
333     pred->Add(ZOSOFN,-1.0);
334    
335     Tpred->Add(ZOSSFN,-1.0);
336     Tpred->Add(ZOSOFP,1.0);
337    
338     Zpred->Add(ZOSOFN,-1.0);
339    
340     //flipped
341     flippedpred->Add(ZOSOFN,1.0);
342     flippedpred->Add(ZOSOFP,-1.0);
343    
344     flippedTpred->Add(ZOSSFP,-1.0);
345     flippedTpred->Add(ZOSOFN,1.0);
346    
347     flippedZpred->Add(ZOSOFP,-1.0);
348     }
349    
350     pred->SetName(predname.c_str());
351     Tpred->SetName(Tpredname.c_str());
352     Zpred->SetName(Zpredname.c_str());
353    
354     flippedpred->SetName(("flipped_"+predname).c_str());
355     flippedTpred->SetName(("flipped_"+Tpredname).c_str());
356     flippedZpred->SetName(("flipped_"+Zpredname).c_str());
357    
358     pred->Write();
359     Tpred->Write();
360     Zpred->Write();
361    
362     flippedpred->Write();
363     flippedTpred->Write();
364     flippedZpred->Write();
365    
366     if(identifier=="") {
367     float stretchfactor=1.0;
368     bool isonpeak=false;
369     if(PlottingSetup::RestrictToMassPeak) {
370     stretchfactor=1.0/3.0;
371     isonpeak=true;
372     }
373    
374     //--------------------------------------------------------statistical uncertainty
375    
376     TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
377     predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
378     predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
379     if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
380     if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
381     if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
382     if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
383     SQRT(predstaterr);
384     TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
385     bgStatUp->Add(predstaterr);
386     bgStatUp->Write();
387     TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDn");
388     bgStatDn->Add(predstaterr,-1);
389     bgStatDn->Write();
390     // delete bgStatDn;
391     // delete bgStatUp;
392     delete predstaterr;
393    
394    
395     TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
396     flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
397     flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
398     if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
399     if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
400     if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
401     if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
402     SQRT(flippedpredstaterr);
403     TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
404     fbgStatUp->Add(predstaterr);
405     fbgStatUp->Write();
406     TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDn");
407     fbgStatDn->Add(predstaterr,-1);
408     fbgStatDn->Write();
409     // delete fbgStatDn;
410     // delete fbgStatUp;
411     delete predstaterr;
412    
413     TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
414     Tpredstaterr->Scale(stretchfactor*stretchfactor);
415     if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
416     if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
417     SQRT(Tpredstaterr);
418     TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
419     TpredStatUp->Add(Tpredstaterr);
420     TpredStatUp->Write();
421     TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDn");
422     TpredStatDn->Add(Tpredstaterr,-1);
423     TpredStatDn->Write();
424     // delete TpredStatDn;
425     // delete TpredStatUp;
426     delete Tpredstaterr;
427    
428     TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
429     fTpredstaterr->Scale(stretchfactor*stretchfactor);
430     if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
431     if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
432     SQRT(fTpredstaterr);
433     TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
434     fTpredStatUp->Add(fTpredstaterr);
435     fTpredStatUp->Write();
436     TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDn");
437     fTpredStatDn->Add(fTpredstaterr,-1);
438     fTpredStatDn->Write();
439     // delete fTpredStatDn;
440     // delete fTpredStatUp;
441     delete fTpredstaterr;
442    
443    
444     TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
445     Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
446     if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
447     if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
448     SQRT(Zpredstaterr);
449     TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
450     ZpredStatUp->Add(Zpredstaterr);
451     ZpredStatUp->Write();
452     TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDn");
453     ZpredStatDn->Add(Zpredstaterr,-1);
454     ZpredStatDn->Write();
455     // delete ZpredStatDn;
456     // delete ZpredStatUp;
457     delete Zpredstaterr;
458    
459     TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
460     Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
461     if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
462     if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
463     SQRT(fTpredstaterr);
464     TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
465     fZpredStatUp->Add(fZpredstaterr);
466     fZpredStatUp->Write();
467     TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDn");
468     fZpredStatDn->Add(fZpredstaterr,-1);
469     fZpredStatDn->Write();
470     // delete fZpredStatDn;
471     // delete fZpredStatUp;
472     delete fZpredstaterr;
473    
474     //--------------------------------------------------------systematic uncertainty
475    
476     TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
477     predsyserr->Add(pred,-1);//getting everything to zero.
478     predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
479     predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
480     predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
481     if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
482     if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
483     if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
484     if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
485     SQRT(predsyserr);
486     TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
487     bgSysUp->Add(predsyserr);
488     bgSysUp->Write();
489     TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDn");
490     bgSysDn->Add(predsyserr,-1);
491     bgSysDn->Write();
492     delete predsyserr;
493    
494     TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
495     fpredsyserr->Add(pred,-1);//getting everything to zero.
496     fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
497     fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
498     fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
499     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
500     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
501     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
502     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
503     SQRT(fpredsyserr);
504     TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
505     fbgSysUp->Add(fpredsyserr);
506     fbgSysUp->Write();
507     TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDn");
508     fbgSysDn->Add(fpredsyserr,-1);
509     fbgSysDn->Write();
510     delete fpredsyserr;
511    
512    
513     TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
514     Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
515     Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
516     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
517     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
518     SQRT(Tpredsyserr);
519     TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
520     TpredSysUp->Add(Tpredsyserr);
521     TpredSysUp->Write();
522     TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDn");
523     TpredSysDn->Add(Tpredsyserr,-1);
524     TpredSysDn->Write();
525     delete Tpredsyserr;
526    
527     TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
528     fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
529     fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
530     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
531     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
532     SQRT(fTpredsyserr);
533     TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
534     fTpredSysUp->Add(fTpredsyserr);
535     fTpredSysUp->Write();
536     TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDn");
537     fTpredSysDn->Add(fTpredsyserr,-1);
538     fTpredSysDn->Write();
539     delete fTpredsyserr;
540    
541    
542     //------------------------------------------------------------------------------------------------------------------------
543     TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
544     Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
545     Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
546     Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
547     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
548     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
549     SQRT(Zpredsyserr);
550     TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
551     ZpredSysUp->Add(Zpredsyserr);
552     ZpredSysUp->Write();
553     TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDn");
554     ZpredSysDn->Add(Zpredsyserr,-1);
555     ZpredSysDn->Write();
556     delete Zpredsyserr;
557    
558     TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
559     fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
560     fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
561     fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
562     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
563     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
564     SQRT(fZpredsyserr);
565     TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
566     fZpredSysUp->Add(fZpredsyserr);
567     fZpredSysUp->Write();
568     TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDn");
569     fZpredSysDn->Add(fZpredsyserr,-1);
570     fZpredSysDn->Write();
571     delete fZpredsyserr;
572     }
573    
574     if(identifier=="") {
575     for(int i=0;i<binning.size()-1;i++) {
576     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;
577     }
578     }
579     delete ZOSSFP;
580     delete ZOSOFP;
581     delete ZOSSFN;
582     delete ZOSOFN;
583    
584     if(PlottingSetup::RestrictToMassPeak) {
585     delete SBOSSFP;
586     delete SBOSOFP;
587     delete SBOSSFN;
588     delete SBOSOFN;
589     }
590     }
591    
592     }
593    
594     ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
595     map < pair<float, float>, map<string, float> > xsec;
596     if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
597    
598     write_info(__FUNCTION__,"Implementing the shape function");
599     time_t timestamp;
600     tm *now;
601     timestamp = time(0);
602     now = localtime(&timestamp);
603     stringstream RunDirectoryS;
604     RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
605     string RunDirectory = RunDirectoryS.str();
606    
607     ensure_directory_exists(RunDirectory);
608    
609     TFile *datafile = new TFile("../StoredShapes.root","READ");
610     if(datafile->IsZombie()) {
611     write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
612     assert(!datafile->IsZombie());
613     }
614     cout << "Run Directory: " << RunDirectory << endl;
615     TFile *limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
616    
617     TIter nextkey(datafile->GetListOfKeys());
618     TKey *key;
619     while ((key = (TKey*)nextkey()))
620     {
621     TH1F *obj =(TH1F*) key->ReadObj();
622     limfile->cd();
623     obj->Write();
624     }
625    
626     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
627    
628     bool signalonly=true;
629     bool dataonly=false;
630    
631     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
632     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
633     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
634     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
635     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
636    
637     // JSU & PDF: global factors
638     write_warning(__FUNCTION__,"NEED TO IMPLEMENT JSU & PDF ");
639    
640     float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
641     float scaledown=0,scaleup=0,scalesyst=0;
642     float mceff=0,mcefferr=0;
643     int flipped=0;
644     int Neventsinfile=10000;//doesn't matter.
645     string informalname="ShapeAnalysis";
646    
647     bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
648    
649     MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
650     if(mceff<0) flipped=1;
651     doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
652     float PDFuncert;
653     int NPdfs = get_npdfs(events);
654     if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
655    
656    
657     prepare_datacard(limfile);
658    
659     limfile->Close();
660     write_info(__FUNCTION__,"limitfile.root and datacard.txt have been generated");
661    
662     write_error(__FUNCTION__,"Implementation of limit calculation is still missing");
663    
664     write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up.");
665     dout << "Everything is saved in " << RunDirectory << endl;
666     // gSystem->Exec(("rm -r "+RunDirectory).c_str());
667     write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please store the information below.");
668     ShapeDroplet alpha;
669     alpha.expected=3;
670    
671     return alpha;
672     }
673    
674     void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
675     dout << "Generating data shape templates ... " << endl;
676     map < pair<float, float>, map<string, float> > xsec;
677     TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
678     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
679     TTree *faketree;
680     bool dataonly=true;
681     bool signalonly=false;
682     string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
683     float jzbpeakerrormc=0;
684     generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
685     // don't need these effects for obs & pred, only for signal!
686     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
687     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
688     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
689     // generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
690    
691    
692     datafile->Close();
693     }
694