ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.13
Committed: Thu May 17 19:52:55 2012 UTC (12 years, 11 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.12: +59 -70 lines
Log Message: 
Significantly sped up shape analysis by using asymptotic upper limits as direct input; not recomputing systematics over again for full cls; two separate shape droplets for asymptotic and final results (return asymptotic ones if no full cls is needed); introduced upper and lower cutoff within which full cls is required; using timeout script to make sure nothing funny is going on in combination tool; overriding standard tmp directory; forcing hadd to replace any existing file; using different subdirectory for full cls to avoid clashes; writing out timing; catching combine tool errors directly; using timeout to detect upper limit cutoff and avoid further conflicts; added final grid point to ensure safe computation; removed additional comma that would freak out bash; doing full cls computation for sms at every point; many more changes but i think space here is limited ...

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