ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.11
Committed: Mon May 14 15:14:03 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.10: +168 -55 lines
Log Message: 
Updated shape implementation

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