ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C
Revision: 1.17
Committed: Fri Sep 7 13:40:44 2012 UTC (12 years, 8 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.16: +1 -1 lines
Log Message: 
Minor fixes for comparisons between signed and unsigned integer expressions

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 buchmann 1.17 while(currentpoint<25*observed && (int)points.size()<npoints) {
54 buchmann 1.12
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 buchmann 1.15 if((int)points.size()==npoints) break;
59 buchmann 1.12 }
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 buchmann 1.14 TH2F* empty2d(string hname) {
300     TH2F *h = new TH2F(hname.c_str(),hname.c_str(),1000,0,1000,1000,0,1000);
301     return h;
302     }
303 buchmann 1.1
304 buchmann 1.14 TH2F* QuickDraw2d(TTree *events, string hname, string variable, string xlabel, string ylabel, TCut cut, string addcut, bool dataormc, float luminosity, map < pair<float, float>, map<string, float> > &xsec) {
305     TH2F *histo = empty2d(hname);
306     histo->Sumw2();
307     stringstream drawcommand;
308     drawcommand << variable << ":mll>>" << hname;
309     if(SUSYScanSpace::SUSYscantype!=mSUGRA) {
310     events->Draw(drawcommand.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight);
311     } else {
312     float totalxs=0;
313     for(int i=0;i<12;i++) {
314     stringstream drawcommand2;
315     drawcommand2 << variable << ">>+" << hname;
316     stringstream mSUGRAweight;
317     float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
318     totalxs+=xs;
319     mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
320     events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
321     }
322     histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
323     }
324     events->Draw("eventNum",addcut.c_str(),"goff");
325     float nevents = events->GetSelectedRows();
326     QuickDrawNevents=nevents;
327     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
328    
329     return histo;
330 buchmann 1.1 }
331    
332     TH1F* Multiply(TH1F *h1, TH1F *h2) {
333     TH1F *h = (TH1F*)h1->Clone();
334     for(int i=1;i<=h1->GetNbinsX();i++) {
335     h->SetBinContent(i,h1->GetBinContent(i) * h2->GetBinContent(i));
336     }
337     return h;
338     }
339    
340 buchmann 1.11
341     string ReduceNumericHistoName(string origname) {
342     int pos = origname.find("__h_");
343     if(pos == -1) return origname;
344     return origname.substr(0,pos);
345     }
346    
347    
348    
349     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) {
350     //weight can be varied to take care of efficiency errors (weightEffUp, weightEffDown)
351     vector<float> mbinning;
352     mbinning.push_back(-8000);
353     for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
354     mbinning.push_back(0);
355 buchmann 1.15 for(int i=0;i<(int)binning.size();i++) mbinning.push_back(binning[i]);
356 buchmann 1.11 mbinning.push_back(8000);
357    
358    
359     TCanvas *shapecan = new TCanvas("shapecan","shapecan");
360     TH1F* Signal;
361    
362     stringstream sInverseCutWeight;
363     sInverseCutWeight << "(1.0/" << (const char*)cutWeight<<")";
364     TCut InverseCutWeight(sInverseCutWeight.str().c_str());
365     if(weight==cutWeight) InverseCutWeight = TCut("1.0");
366     if(!dotree) {
367     //dealing with ALL MC samples (when we save the standard file)
368     THStack McPredicted = allsamples.DrawStack("McPred",mcjzb,mbinning, "JZB [GeV]", "events", ((cutmass&&cutOSSF&&JetCut)*(weight*InverseCutWeight)),mc,luminosity);
369     int nhists = McPredicted.GetHists()->GetEntries();
370     for (int i = 0; i < nhists; i++) {
371     TH1* hist = static_cast<TH1*>(McPredicted.GetHists()->At(i));
372     // cout << hist->GetName() << " : " << hist->Integral() << endl;
373     // cout << " " << ReduceNumericHistoName(hist->GetName()) << endl;
374     if(identifier=="StatUp") {
375     float appliedweight = hist->Integral() / hist->GetEntries();
376     for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)+appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
377     }
378     if(identifier=="StatDown") {
379     float appliedweight = hist->Integral() / hist->GetEntries();
380     for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)-appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
381     }
382     hist->SetName(("mc_"+ReduceNumericHistoName(hist->GetName())+"_"+identifier).c_str());
383     hist->Write();
384     }
385     } else {
386     string histoname="mc_signal";
387     if(identifier!="") histoname=histoname+"_"+identifier;
388     Signal = QuickDraw(signalevents,histoname,mcjzb,binning,"JZB", "events",((cutmass&&cutOSSF&&JetCut&&basiccut))*(weight*InverseCutWeight),addcut,mc,luminosity,xsec);
389     if(identifier=="StatUp") {
390     float appliedweight = Signal->Integral() / Signal->GetEntries();
391     for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)+appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
392     }
393     if(identifier=="StatDown") {
394     float appliedweight = Signal->Integral() / Signal->GetEntries();
395     for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)-appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
396     }
397     Signal->Write();
398     }
399    
400     /*
401     *
402     * Jet energy scale (easy, use different JetCut)
403     * JZB Scale Uncertainty (will be a number - signal only)
404     * Efficiency (will be weightEffUp, weightEffDown)
405     * PDFs (signal only) -> Will be a number yet again, computed in the traditional way
406     */
407     delete shapecan;
408    
409     }
410    
411    
412 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) {
413 buchmann 1.4
414 buchmann 1.1 binning.push_back(8000);
415     limfile->cd();
416 buchmann 1.4 dout << "Creating shape templates ";
417 buchmann 1.1 if(identifier!="") dout << "for "<<identifier;
418     dout << " ... " << endl;
419    
420     TCut limitnJetcut;
421 buchmann 1.4
422 buchmann 1.1 if(JES==noJES) limitnJetcut=cutnJets;
423     else {
424     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
425     if(JES==JESup) limitnJetcut=cutnJetsJESup;
426     }
427    
428     float zjetsestimateuncert=zjetsestimateuncertONPEAK;
429     float emuncert=emuncertONPEAK;
430     float emsidebanduncert=emsidebanduncertONPEAK;
431     float eemmsidebanduncert=eemmsidebanduncertONPEAK;
432    
433     if(!PlottingSetup::RestrictToMassPeak) {
434     zjetsestimateuncert=zjetsestimateuncertOFFPEAK;
435     emuncert=emuncertOFFPEAK;
436     emsidebanduncert=emsidebanduncertOFFPEAK;
437     eemmsidebanduncert=eemmsidebanduncertOFFPEAK;
438     }
439    
440     if(signalonly) {
441 buchmann 1.7 dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
442 buchmann 1.1 TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
443     TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
444 buchmann 1.10 TH1F *ZOSOFP;
445     TH1F *ZOSOFN;
446 buchmann 1.14
447     TH2F *ZOSSFP2d = QuickDraw2d(signalevents,"ZOSSFP2d",mcjzb, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
448     TH2F *ZOSSFN2d = QuickDraw2d(signalevents,"ZOSSFN2d","-"+mcjzb,"JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
449     TH2F *ZOSOFP2d;
450     TH2F *ZOSOFN2d;
451    
452 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) {
453     ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
454     ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
455 buchmann 1.14
456     ZOSOFP2d = QuickDraw2d(signalevents,"ZOSOFP2d",mcjzb, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
457     ZOSOFN2d = QuickDraw2d(signalevents,"ZOSOFN2d","-"+mcjzb, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
458 buchmann 1.10 }
459 buchmann 1.1
460     TH1F *SBOSSFP;
461     TH1F *SBOSOFP;
462     TH1F *SBOSSFN;
463     TH1F *SBOSOFN;
464    
465 buchmann 1.14 TH2F *SBOSSFP2d;
466     TH2F *SBOSOFP2d;
467     TH2F *SBOSSFN2d;
468     TH2F *SBOSOFN2d;
469    
470 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis&&PlottingSetup::UseSidebandsForcJZB) {
471 buchmann 1.1 SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
472     SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
473     SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
474     SBOSOFN = QuickDraw(signalevents,"SBOSOFN","-"+mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
475 buchmann 1.14
476     SBOSSFP2d = QuickDraw2d(signalevents,"SBOSSFP2d",mcjzb,"JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
477     SBOSOFP2d = QuickDraw2d(signalevents,"SBOSOFP2d",mcjzb,"JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
478     SBOSSFN2d = QuickDraw2d(signalevents,"SBOSSFN2d","-"+mcjzb,"JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
479     SBOSOFN2d = QuickDraw2d(signalevents,"SBOSOFN2d","-"+mcjzb,"JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
480 buchmann 1.1 }
481    
482     string signalname="signal";
483 buchmann 1.14 if(identifier!="") signalname="signal_"+identifier;
484 buchmann 1.1
485     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
486     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
487    
488 buchmann 1.14 TH2F *Lobs2d = empty2d("Lobs2d");
489     TH2F *Lpred2d = empty2d("Lobs2d");
490    
491 buchmann 1.1 TH1F *flippedLobs = new TH1F("flippedLobs","flippedLobs",binning.size()-1,&binning[0]);
492     TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
493    
494 buchmann 1.14 TH2F *flippedLobs2d = empty2d("flippedLobs2d");
495     TH2F *flippedLpred2d = empty2d("flippedLpred2d");
496    
497 buchmann 1.1 Lobs->Add(ZOSSFP);
498 buchmann 1.14 Lobs2d->Add(ZOSSFP2d);
499     if(!PlottingSetup::FullMCAnalysis) {
500     Lpred->Add(ZOSSFN);
501     Lpred2d->Add(ZOSSFN2d);
502     }
503 buchmann 1.1
504 buchmann 1.7 dout << "SITUATION FOR SIGNAL: " << endl;
505 buchmann 1.10 if(PlottingSetup::FullMCAnalysis) {
506 buchmann 1.14 dout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << endl;
507 buchmann 1.10 } else {
508 buchmann 1.14 dout << " OSSF JZB> 0 : " << ZOSSFP->Integral() << " JZB < 0 :" << ZOSSFN->Integral() << endl;
509     dout << " OSOF JZB> 0 : " << ZOSOFP->Integral() << " JZB < 0 :" << ZOSOFN->Integral() << endl;
510    
511 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis&&PlottingSetup::UseSidebandsForcJZB) {
512 buchmann 1.14 dout << " OSSF SB JZB> 0 : " << SBOSSFP->Integral() << " JZB < 0 :" << SBOSSFN->Integral() << endl;
513     dout << " OSOF SB JZB> 0 : " << SBOSOFP->Integral() << " JZB < 0 :" << SBOSOFN->Integral() << endl;
514     }
515 buchmann 1.2 }
516    
517 buchmann 1.1 flippedLobs->Add(ZOSSFN);
518 buchmann 1.14 flippedLobs2d->Add(ZOSSFN2d);
519    
520     if(!PlottingSetup::FullMCAnalysis) {
521     flippedLpred->Add(ZOSSFP);
522     flippedLpred2d->Add(ZOSSFP2d);
523     }
524 buchmann 1.1
525 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) {
526 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
527 buchmann 1.14 Lpred->Add(ZOSOFP,1.0/3);
528     Lpred->Add(ZOSOFN,-1.0/3);
529     Lpred->Add(SBOSSFP,1.0/3);
530     Lpred->Add(SBOSSFN,-1.0/3);
531     Lpred->Add(SBOSOFP,1.0/3);
532     Lpred->Add(SBOSOFN,-1.0/3);
533    
534     //flipped prediction
535     flippedLpred->Add(ZOSOFP,-1.0/3);
536     flippedLpred->Add(ZOSOFN,1.0/3);
537     flippedLpred->Add(SBOSSFP,-1.0/3);
538     flippedLpred->Add(SBOSSFN,1.0/3);
539     flippedLpred->Add(SBOSOFP,-1.0/3);
540     flippedLpred->Add(SBOSOFN,1.0/3);
541    
542     //2d stuff: regular
543     Lpred2d->Add(ZOSOFP2d,1.0/3);
544     Lpred2d->Add(ZOSOFN2d,-1.0/3);
545     Lpred2d->Add(SBOSSFP2d,1.0/3);
546     Lpred2d->Add(SBOSSFN2d,-1.0/3);
547     Lpred2d->Add(SBOSOFP2d,1.0/3);
548     Lpred2d->Add(SBOSOFN2d,-1.0/3);
549    
550     //3d stuff: flipped prediction
551     flippedLpred2d->Add(ZOSOFP2d,-1.0/3);
552     flippedLpred2d->Add(ZOSOFN2d,1.0/3);
553     flippedLpred2d->Add(SBOSSFP2d,-1.0/3);
554     flippedLpred2d->Add(SBOSSFN2d,1.0/3);
555     flippedLpred2d->Add(SBOSOFP2d,-1.0/3);
556     flippedLpred2d->Add(SBOSOFN2d,1.0/3);
557 buchmann 1.10 } else {
558 buchmann 1.14 Lpred->Add(ZOSOFP,1.0);
559     Lpred->Add(ZOSOFN,-1.0);
560    
561     //flipped prediction
562     flippedLpred->Add(ZOSOFP,-1.0);
563     flippedLpred->Add(ZOSOFN,1.0);
564    
565     //2d stuff
566     Lpred2d->Add(ZOSOFP2d,1.0);
567     Lpred2d->Add(ZOSOFN2d,-1.0);
568    
569     //flipped prediction
570     flippedLpred2d->Add(ZOSOFP2d,-1.0);
571     flippedLpred2d->Add(ZOSOFN2d,1.0);
572 buchmann 1.10 }
573 buchmann 1.1 }
574    
575 buchmann 1.2 TH1F *signal = (TH1F*)Lobs->Clone("signal");
576 buchmann 1.14 TH1F *signal2d = (TH1F*)Lobs2d->Clone("signal2d");
577    
578 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
579 buchmann 1.1 signal->SetName(signalname.c_str());
580 buchmann 1.2 signal->SetTitle(signalname.c_str());
581 buchmann 1.1 signal->Write();
582    
583 buchmann 1.14 if(!PlottingSetup::FullMCAnalysis) signal2d->Add(Lpred2d,-1);
584     signal2d->SetName((signalname+"2d").c_str());
585     signal2d->SetTitle((signalname+"2d").c_str());
586     signal2d->Write();
587    
588 buchmann 1.1 TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
589 buchmann 1.10 if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
590 buchmann 1.1 flippedsignal->SetName(("flipped_"+signalname).c_str());
591     flippedsignal->Write();
592    
593 buchmann 1.14 TH2F *flippedsignal2d = (TH2F*)flippedLobs2d->Clone();
594     if(!PlottingSetup::FullMCAnalysis) flippedsignal2d->Add(flippedLpred,-1);
595     flippedsignal2d->SetName(("flipped2d_"+signalname).c_str());
596     flippedsignal2d->Write();
597    
598 buchmann 1.1 if(identifier=="") {
599     TH1F *signalStatUp = (TH1F*)signal->Clone();
600     signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
601     signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
602    
603     TH1F *signalStatDn = (TH1F*)signal->Clone();
604     signalStatDn->SetName(((string)signal->GetName()+"_StatDown").c_str());
605     signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
606    
607 buchmann 1.14 TH2F *signalStatUp2d = (TH2F*)signal->Clone();
608     signalStatUp2d->SetName(((string)signal2d->GetName()+"_StatUp").c_str());
609     signalStatUp2d->SetTitle(((string)signal2d->GetTitle()+"_StatUp").c_str());
610    
611     TH2F *signalStatDn2d = (TH2F*)signal->Clone();
612     signalStatDn2d->SetName(((string)signal2d->GetName()+"_StatDown").c_str());
613     signalStatDn2d->SetTitle(((string)signal2d->GetTitle()+"_StatDown").c_str());
614    
615    
616 buchmann 1.1 for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
617 buchmann 1.10 float staterr;
618     if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
619     else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
620    
621     if(!PlottingSetup::FullMCAnalysis) {
622     dout << "Stat err in bin " << i << " : " << staterr << endl;
623     dout << " prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
624     dout << " we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
625     }
626 buchmann 1.5 if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
627     else signalStatDn->SetBinContent(i,0);
628 buchmann 1.2 signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
629     signal->SetBinError(i,staterr);
630 buchmann 1.1 }
631    
632 buchmann 1.14 for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
633     for(int j=1;j<=signalStatDn->GetNbinsY();j++) {
634     float staterr;
635     if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred2d->GetBinContent(i,j) + Lobs2d->GetBinContent(i,j));
636     else staterr = TMath::Sqrt(Lobs2d->GetBinContent(i,j));
637    
638     if(!PlottingSetup::FullMCAnalysis) {
639     dout << "Stat err in bin " << i << ", " << j << " : " << staterr << endl;
640     dout << " prediction: " << Lpred2d->GetBinContent(i,j) << " , observation: " << Lobs2d->GetBinContent(i,j) << " --> signal: " << signal2d->GetBinContent(i,j) << endl;
641     dout << " we obtain : " << signal2d->GetBinContent(i,j)-staterr << " , " << signal2d->GetBinContent(i,j)+staterr << endl;
642     }
643     if(signal2d->GetBinContent(i,j)-staterr>0) signalStatDn2d->SetBinContent(i,j,signal2d->GetBinContent(i,j)-staterr);
644     else signalStatDn2d->SetBinContent(i,j,0);
645     signalStatUp2d->SetBinContent(i,signal2d->GetBinContent(i,j)+staterr);
646     signal2d->SetBinError(i,j,staterr);
647     }
648     }
649    
650     //----------------------------------- ******** GOTTEN HERE *******---------------------------------------
651    
652 buchmann 1.1 signalStatDn->Write();
653     signalStatUp->Write();
654    
655     delete signalStatDn;
656     delete signalStatUp;
657    
658     TH1F *flippedsignalStatUp = (TH1F*)flippedsignal->Clone();
659     flippedsignalStatUp->SetName(((string)flippedsignal->GetName()+"_StatUp").c_str());
660     flippedsignalStatUp->SetTitle(((string)flippedsignal->GetTitle()+"_StatUp").c_str());
661    
662     TH1F *flippedsignalStatDn = (TH1F*)flippedsignal->Clone();
663     flippedsignalStatDn->SetName(((string)flippedsignal->GetName()+"_StatDown").c_str());
664     flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
665    
666     for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
667 buchmann 1.10 float staterr;
668     if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
669     else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
670 buchmann 1.6 if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
671     else flippedsignalStatDn->SetBinContent(i,0);
672 buchmann 1.2 flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
673     flippedsignal->SetBinError(i,staterr);
674 buchmann 1.1 }
675    
676     flippedsignalStatDn->Write();
677     flippedsignalStatUp->Write();
678    
679     delete flippedsignalStatDn;
680     delete flippedsignalStatUp;
681     }
682    
683     delete ZOSSFP;
684     delete ZOSOFP;
685     delete ZOSSFN;
686     delete ZOSOFN;
687    
688 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
689 buchmann 1.1 delete SBOSSFP;
690     delete SBOSOFP;
691     delete SBOSSFN;
692     delete SBOSOFN;
693     }
694    
695     delete Lobs;
696     delete Lpred;
697     delete flippedLobs;
698     delete flippedLpred;
699     }
700    
701    
702     if(dataonly) {
703 buchmann 1.7 dout << "Processing data with datajzb: " << datajzb << endl;
704 buchmann 1.1 TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
705     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
706     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
707     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
708    
709     TH1F *SBOSSFP;
710     TH1F *SBOSOFP;
711     TH1F *SBOSSFN;
712     TH1F *SBOSOFN;
713    
714 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
715 buchmann 1.1 SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
716     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
717     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
718     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,data,luminosity);
719     }
720    
721     string obsname="data_obs";
722     string predname="background";
723     string Zpredname="ZJetsBackground";
724     string Tpredname="TTbarBackground";
725     if(identifier!="") {
726     obsname=("data_"+identifier);
727     predname=("background_"+identifier);
728     Zpredname=("ZJetsBackground_"+identifier);
729     Tpredname=("TTbarBackground_"+identifier);
730     }
731    
732     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
733     obs->SetName(obsname.c_str());
734     obs->Write();
735    
736     TH1F *flippedobs = (TH1F*)ZOSSFN->Clone("flipped_observation");
737     flippedobs->SetName(("flipped_"+obsname).c_str());
738     flippedobs->Write();
739    
740     TH1F *Tpred = (TH1F*)ZOSSFN->Clone("TTbarprediction");
741     TH1F *Zpred = (TH1F*)ZOSSFN->Clone("ZJetsprediction");
742    
743     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
744    
745     TH1F *flippedTpred = (TH1F*)ZOSSFP->Clone("flipped_TTbarprediction");
746     TH1F *flippedZpred = (TH1F*)ZOSSFP->Clone("flipped_ZJetsprediction");
747    
748     TH1F *flippedpred = (TH1F*)ZOSSFP->Clone("flipped_prediction");
749 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
750 buchmann 1.1 pred->Add(ZOSOFP,1.0/3);
751     pred->Add(ZOSOFN,-1.0/3);
752     pred->Add(SBOSSFP,1.0/3);
753     pred->Add(SBOSSFN,-1.0/3);
754     pred->Add(SBOSOFP,1.0/3);
755     pred->Add(SBOSOFN,-1.0/3);
756    
757     Tpred->Add(ZOSSFN,-1.0);
758     Tpred->Add(ZOSOFP,1.0/3);
759     Tpred->Add(SBOSSFP,1.0/3);
760     Tpred->Add(SBOSOFP,1.0/3);
761    
762     Zpred->Add(ZOSOFN,-1.0/3);
763     Zpred->Add(SBOSSFN,-1.0/3);
764     Zpred->Add(SBOSOFN,-1.0/3);
765    
766     //flipped
767     flippedpred->Add(ZOSOFP,-1.0/3);
768     flippedpred->Add(ZOSOFN,1.0/3);
769     flippedpred->Add(SBOSSFP,-1.0/3);
770     flippedpred->Add(SBOSSFN,1.0/3);
771     flippedpred->Add(SBOSOFP,-1.0/3);
772     flippedpred->Add(SBOSOFN,1.0/3);
773    
774     flippedTpred->Add(ZOSSFP,-1.0);
775     flippedTpred->Add(ZOSOFN,1.0/3);
776     flippedTpred->Add(SBOSSFN,1.0/3);
777     flippedTpred->Add(SBOSOFN,1.0/3);
778    
779     flippedZpred->Add(ZOSOFP,-1.0/3);
780     flippedZpred->Add(SBOSSFP,-1.0/3);
781     flippedZpred->Add(SBOSOFP,-1.0/3);
782     } else {
783     pred->Add(ZOSOFP,1.0);
784     pred->Add(ZOSOFN,-1.0);
785    
786     Tpred->Add(ZOSSFN,-1.0);
787     Tpred->Add(ZOSOFP,1.0);
788    
789     Zpred->Add(ZOSOFN,-1.0);
790    
791     //flipped
792     flippedpred->Add(ZOSOFN,1.0);
793     flippedpred->Add(ZOSOFP,-1.0);
794    
795     flippedTpred->Add(ZOSSFP,-1.0);
796     flippedTpred->Add(ZOSOFN,1.0);
797    
798     flippedZpred->Add(ZOSOFP,-1.0);
799     }
800    
801     pred->SetName(predname.c_str());
802     Tpred->SetName(Tpredname.c_str());
803     Zpred->SetName(Zpredname.c_str());
804    
805     flippedpred->SetName(("flipped_"+predname).c_str());
806     flippedTpred->SetName(("flipped_"+Tpredname).c_str());
807     flippedZpred->SetName(("flipped_"+Zpredname).c_str());
808    
809     pred->Write();
810     Tpred->Write();
811     Zpred->Write();
812    
813     flippedpred->Write();
814     flippedTpred->Write();
815     flippedZpred->Write();
816    
817     if(identifier=="") {
818     float stretchfactor=1.0;
819     bool isonpeak=false;
820 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
821 buchmann 1.1 stretchfactor=1.0/3.0;
822     isonpeak=true;
823     }
824    
825     //--------------------------------------------------------statistical uncertainty
826    
827     TH1F *predstaterr = (TH1F*)ZOSSFN->Clone("predstaterr");
828     predstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
829     predstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
830     if(isonpeak) predstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
831     if(isonpeak) predstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
832     if(isonpeak) predstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
833     if(isonpeak) predstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
834     SQRT(predstaterr);
835     TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
836     bgStatUp->Add(predstaterr);
837 buchmann 1.6 EliminateNegativeEntries(bgStatUp);
838 buchmann 1.1 bgStatUp->Write();
839 buchmann 1.2 TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
840 buchmann 1.1 bgStatDn->Add(predstaterr,-1);
841 buchmann 1.6 EliminateNegativeEntries(bgStatDn);
842 buchmann 1.1 bgStatDn->Write();
843     // delete bgStatDn;
844     // delete bgStatUp;
845     delete predstaterr;
846    
847    
848     TH1F *flippedpredstaterr = (TH1F*)ZOSSFP->Clone("predstaterr");
849     flippedpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
850     flippedpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
851     if(isonpeak) flippedpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
852     if(isonpeak) flippedpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
853     if(isonpeak) flippedpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
854     if(isonpeak) flippedpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
855     SQRT(flippedpredstaterr);
856     TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
857     fbgStatUp->Add(predstaterr);
858 buchmann 1.6 EliminateNegativeEntries(fbgStatUp);
859 buchmann 1.1 fbgStatUp->Write();
860 buchmann 1.2 TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
861 buchmann 1.1 fbgStatDn->Add(predstaterr,-1);
862 buchmann 1.6 EliminateNegativeEntries(fbgStatDn);
863 buchmann 1.1 fbgStatDn->Write();
864     // delete fbgStatDn;
865     // delete fbgStatUp;
866     delete predstaterr;
867    
868     TH1F *Tpredstaterr = (TH1F*)ZOSOFP->Clone("Tpredstaterr");
869     Tpredstaterr->Scale(stretchfactor*stretchfactor);
870     if(isonpeak) Tpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
871     if(isonpeak) Tpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
872     SQRT(Tpredstaterr);
873     TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
874     TpredStatUp->Add(Tpredstaterr);
875 buchmann 1.6 EliminateNegativeEntries(TpredStatUp);
876 buchmann 1.1 TpredStatUp->Write();
877 buchmann 1.2 TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
878 buchmann 1.1 TpredStatDn->Add(Tpredstaterr,-1);
879 buchmann 1.6 EliminateNegativeEntries(TpredStatDn);
880 buchmann 1.1 TpredStatDn->Write();
881     // delete TpredStatDn;
882     // delete TpredStatUp;
883     delete Tpredstaterr;
884    
885     TH1F *fTpredstaterr = (TH1F*)ZOSOFN->Clone("flipped_Tpredstaterr");
886     fTpredstaterr->Scale(stretchfactor*stretchfactor);
887     if(isonpeak) fTpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
888     if(isonpeak) fTpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
889     SQRT(fTpredstaterr);
890     TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
891     fTpredStatUp->Add(fTpredstaterr);
892 buchmann 1.6 EliminateNegativeEntries(fTpredStatUp);
893 buchmann 1.1 fTpredStatUp->Write();
894 buchmann 1.2 TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
895 buchmann 1.1 fTpredStatDn->Add(fTpredstaterr,-1);
896 buchmann 1.6 EliminateNegativeEntries(fTpredStatDn);
897 buchmann 1.1 fTpredStatDn->Write();
898     // delete fTpredStatDn;
899     // delete fTpredStatUp;
900     delete fTpredstaterr;
901    
902    
903     TH1F *Zpredstaterr = (TH1F*)ZOSSFN->Clone("ZJetsstaterr");
904     Zpredstaterr->Add(ZOSOFN,stretchfactor*stretchfactor);
905     if(isonpeak) Zpredstaterr->Add(SBOSSFN,stretchfactor*stretchfactor);
906     if(isonpeak) Zpredstaterr->Add(SBOSOFN,stretchfactor*stretchfactor);
907     SQRT(Zpredstaterr);
908     TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
909     ZpredStatUp->Add(Zpredstaterr);
910 buchmann 1.6 EliminateNegativeEntries(ZpredStatUp);
911 buchmann 1.1 ZpredStatUp->Write();
912 buchmann 1.2 TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
913 buchmann 1.1 ZpredStatDn->Add(Zpredstaterr,-1);
914 buchmann 1.6 EliminateNegativeEntries(ZpredStatDn);
915 buchmann 1.1 ZpredStatDn->Write();
916     // delete ZpredStatDn;
917     // delete ZpredStatUp;
918     delete Zpredstaterr;
919    
920     TH1F *fZpredstaterr = (TH1F*)ZOSSFP->Clone("flipped_ZJetsstaterr");
921     Zpredstaterr->Add(ZOSOFP,stretchfactor*stretchfactor);
922     if(isonpeak) fZpredstaterr->Add(SBOSSFP,stretchfactor*stretchfactor);
923     if(isonpeak) fZpredstaterr->Add(SBOSOFP,stretchfactor*stretchfactor);
924     SQRT(fTpredstaterr);
925     TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
926     fZpredStatUp->Add(fZpredstaterr);
927 buchmann 1.6 EliminateNegativeEntries(fZpredStatUp);
928 buchmann 1.1 fZpredStatUp->Write();
929 buchmann 1.2 TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
930 buchmann 1.1 fZpredStatDn->Add(fZpredstaterr,-1);
931 buchmann 1.6 EliminateNegativeEntries(fZpredStatDn);
932 buchmann 1.1 fZpredStatDn->Write();
933     // delete fZpredStatDn;
934     // delete fZpredStatUp;
935     delete fZpredstaterr;
936    
937     //--------------------------------------------------------systematic uncertainty
938    
939     TH1F *predsyserr = (TH1F*)pred->Clone("SysError");
940     predsyserr->Add(pred,-1);//getting everything to zero.
941     predsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
942     predsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
943     predsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
944     if(isonpeak) predsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
945     if(isonpeak) predsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
946     if(isonpeak) predsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
947     if(isonpeak) predsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
948     SQRT(predsyserr);
949     TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
950     bgSysUp->Add(predsyserr);
951 buchmann 1.6 EliminateNegativeEntries(bgSysUp);
952 buchmann 1.1 bgSysUp->Write();
953 buchmann 1.2 TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
954 buchmann 1.1 bgSysDn->Add(predsyserr,-1);
955 buchmann 1.6 EliminateNegativeEntries(bgSysDn);
956 buchmann 1.1 bgSysDn->Write();
957     delete predsyserr;
958    
959     TH1F *fpredsyserr = (TH1F*)pred->Clone("SysError");
960     fpredsyserr->Add(pred,-1);//getting everything to zero.
961     fpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
962     fpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
963     fpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
964     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
965     if(isonpeak) fpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
966     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
967     if(isonpeak) fpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
968     SQRT(fpredsyserr);
969     TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
970     fbgSysUp->Add(fpredsyserr);
971 buchmann 1.6 EliminateNegativeEntries(fbgSysUp);
972 buchmann 1.1 fbgSysUp->Write();
973 buchmann 1.2 TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
974 buchmann 1.1 fbgSysDn->Add(fpredsyserr,-1);
975 buchmann 1.6 EliminateNegativeEntries(fbgSysDn);
976 buchmann 1.1 fbgSysDn->Write();
977     delete fpredsyserr;
978    
979    
980     TH1F *Tpredsyserr = (TH1F*)Tpred->Clone("SysError");
981     Tpredsyserr->Add(Tpred,-1);//getting everything to zero.
982     Tpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
983     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
984     if(isonpeak) Tpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
985     SQRT(Tpredsyserr);
986     TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
987     TpredSysUp->Add(Tpredsyserr);
988 buchmann 1.6 EliminateNegativeEntries(TpredSysUp);
989 buchmann 1.1 TpredSysUp->Write();
990 buchmann 1.2 TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
991 buchmann 1.1 TpredSysDn->Add(Tpredsyserr,-1);
992 buchmann 1.6 EliminateNegativeEntries(TpredSysDn);
993 buchmann 1.1 TpredSysDn->Write();
994     delete Tpredsyserr;
995    
996     TH1F *fTpredsyserr = (TH1F*)Tpred->Clone("SysError");
997     fTpredsyserr->Add(Tpred,-1);//getting everything to zero.
998     fTpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
999     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
1000     if(isonpeak) fTpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
1001     SQRT(fTpredsyserr);
1002     TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
1003     fTpredSysUp->Add(fTpredsyserr);
1004 buchmann 1.6 EliminateNegativeEntries(fTpredSysUp);
1005 buchmann 1.1 fTpredSysUp->Write();
1006 buchmann 1.2 TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
1007 buchmann 1.1 fTpredSysDn->Add(fTpredsyserr,-1);
1008 buchmann 1.6 EliminateNegativeEntries(fTpredSysDn);
1009 buchmann 1.1 fTpredSysDn->Write();
1010     delete fTpredsyserr;
1011    
1012    
1013     //------------------------------------------------------------------------------------------------------------------------
1014     TH1F *Zpredsyserr = (TH1F*)Zpred->Clone("SysError");
1015     Zpredsyserr->Add(Zpred,-1);//getting everything to zero.
1016     Zpredsyserr->Add(Multiply(ZOSSFN,ZOSSFN),zjetsestimateuncert*zjetsestimateuncert);
1017     Zpredsyserr->Add(Multiply(ZOSOFN,ZOSOFN),stretchfactor*stretchfactor*emuncert*emuncert);
1018     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSSFN,SBOSSFN),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
1019     if(isonpeak) Zpredsyserr->Add(Multiply(SBOSOFN,SBOSOFN),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
1020     SQRT(Zpredsyserr);
1021     TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
1022     ZpredSysUp->Add(Zpredsyserr);
1023 buchmann 1.6 EliminateNegativeEntries(ZpredSysUp);
1024 buchmann 1.1 ZpredSysUp->Write();
1025 buchmann 1.2 TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
1026 buchmann 1.1 ZpredSysDn->Add(Zpredsyserr,-1);
1027 buchmann 1.6 EliminateNegativeEntries(ZpredSysDn);
1028 buchmann 1.1 ZpredSysDn->Write();
1029     delete Zpredsyserr;
1030    
1031     TH1F *fZpredsyserr = (TH1F*)Zpred->Clone("SysError");
1032     fZpredsyserr->Add(Zpred,-1);//getting everything to zero.
1033     fZpredsyserr->Add(Multiply(ZOSSFP,ZOSSFP),zjetsestimateuncert*zjetsestimateuncert);
1034     fZpredsyserr->Add(Multiply(ZOSOFP,ZOSOFP),stretchfactor*stretchfactor*emuncert*emuncert);
1035     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSSFP,SBOSSFP),stretchfactor*stretchfactor*eemmsidebanduncert*eemmsidebanduncert);
1036     if(isonpeak) fZpredsyserr->Add(Multiply(SBOSOFP,SBOSOFP),stretchfactor*stretchfactor*emsidebanduncert*emsidebanduncert);
1037     SQRT(fZpredsyserr);
1038     TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
1039     fZpredSysUp->Add(fZpredsyserr);
1040 buchmann 1.6 EliminateNegativeEntries(fZpredSysUp);
1041 buchmann 1.1 fZpredSysUp->Write();
1042 buchmann 1.2 TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
1043 buchmann 1.1 fZpredSysDn->Add(fZpredsyserr,-1);
1044 buchmann 1.6 EliminateNegativeEntries(fZpredSysDn);
1045 buchmann 1.1 fZpredSysDn->Write();
1046     delete fZpredsyserr;
1047     }
1048    
1049 buchmann 1.2 /*if(identifier=="") {
1050 buchmann 1.1 for(int i=0;i<binning.size()-1;i++) {
1051 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;
1052 buchmann 1.1 }
1053 buchmann 1.2 }*/
1054 buchmann 1.1 delete ZOSSFP;
1055     delete ZOSOFP;
1056     delete ZOSSFN;
1057     delete ZOSOFN;
1058    
1059 buchmann 1.16 if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
1060 buchmann 1.1 delete SBOSSFP;
1061     delete SBOSOFP;
1062     delete SBOSSFN;
1063     delete SBOSOFN;
1064     }
1065     }
1066 buchmann 1.14
1067    
1068     //----------------------------------- ******** GOTTEN HERE *******---------------------------------------
1069    
1070    
1071    
1072    
1073    
1074 buchmann 1.1 }
1075    
1076 buchmann 1.12 void DoFullCls(string RunDirectory,float firstGuess) {
1077     vector<float> epoints = get_expected_points(firstGuess,25,RunDirectory,true);//get 20 points;
1078     ofstream RealScript;
1079     dout << "Going to write the full CLs script to " << RunDirectory << "/LimitScript.sh" << endl;
1080     RealScript.open ((RunDirectory+"/LimitScript.sh").c_str());
1081     RealScript << "#!/bin/bash \n";
1082     RealScript << "\n";
1083     RealScript << "RunDirectory=" << RunDirectory << "\n";
1084     RealScript << "CMSSWDirectory=" << PlottingSetup::CMSSW_dir << "\n";
1085 buchmann 1.13 RealScript << "ORIGTMP=$TMP\n";
1086     RealScript << "ORIGTMPDIR=$TMPDIR\n";
1087     RealScript << "export TMP=" << RunDirectory << "\n";
1088     RealScript << "export TMPDIR=" << RunDirectory << "\n";
1089 buchmann 1.12 RealScript << "\n";
1090     RealScript << "echo \"Going to compute limit in $RunDirectory\"\n";
1091     RealScript << "ORIGDIR=`pwd`\n";
1092     RealScript << "\n";
1093     RealScript << "pushd $CMSSWDirectory\n";
1094     RealScript << "cd ~/final_production_2011/CMSSW_4_2_8/src/HiggsAnalysis/\n";
1095     RealScript << "origscramarch=$SCRAM_ARCH\n";
1096     RealScript << "origbase=$CMSSW_BASE\n";
1097     RealScript << "export SCRAM_ARCH=slc5_amd64_gcc434\n";
1098     RealScript << "eval `scram ru -sh`\n";
1099     RealScript << "\n";
1100     RealScript << "mkdir -p $RunDirectory\n";
1101     RealScript << "cd $RunDirectory\n";
1102     RealScript << "echo `pwd`\n";
1103     RealScript << "mkdir -p FullCLs\n";
1104     RealScript << "cd FullCLs\n";
1105     RealScript << "\n";
1106     RealScript << "combineEXE=\"${CMSSW_BASE}/bin/${SCRAM_ARCH}/combine\"\n";
1107     RealScript << "\n";
1108 buchmann 1.13 RealScript << "dobreak=0\n";
1109     RealScript << "npoints=0\n";
1110 buchmann 1.12 RealScript << "time for i in {";
1111     RealScript << epoints[0]/4 << ","; // adding a VERY VERY low point just to be totally safe
1112     RealScript << epoints[0]/2 << ","; // adding a VERY low point just to be safe
1113 buchmann 1.15 for(int i=0;i<(int)epoints.size();i++) RealScript << epoints[i] << ",";
1114 buchmann 1.12 RealScript << 2*epoints[epoints.size()-1] << ","; // adding a VERY high point just to be safe
1115 buchmann 1.13 RealScript << 4*epoints[epoints.size()-1]; // adding a VERY VERY high point just to be totally safe
1116     RealScript << "}; do \n";
1117     RealScript << "let \"npoints = npoints +1\"\n";
1118     RealScript << "if [[ $dobreak -gt 0 ]]; then \n";
1119     RealScript << " continue \n";
1120     RealScript << "fi \n";
1121     RealScript << "echo -e \" Going to compute CLs for $i \\n\"\n";
1122     RealScript << "echo $(date +%d%b%Y,%H:%M:%S)\n";
1123     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";
1124     // 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";
1125     RealScript << "errorsize=$?\n";
1126     RealScript << "rm " << RunDirectory << "/rstat* 2>/dev/null \n";
1127     RealScript << "if [[ $errorsize -gt 0 ]]; then \n";
1128     RealScript << " echo \"Apparently something went wrong (had to be aborted)\"\n";
1129     RealScript << " if [[ $npoints -gt 10 ]]; then \n";
1130     RealScript << " dobreak=1 \n";
1131     RealScript << " fi \n";
1132     RealScript << "fi \n";
1133     RealScript << "done\n";
1134    
1135 buchmann 1.12 RealScript << "\n";
1136 buchmann 1.13 RealScript << "hadd -f FullCLs.root higgsCombine*.root\n";
1137 buchmann 1.12 RealScript << "mkdir originalFiles\n";
1138     RealScript << "mv higgsCombine* originalFiles\n";
1139     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root\n";
1140     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.5\n";
1141     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.84\n";
1142     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.16\n";
1143     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.975\n";
1144     RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.025\n";
1145     RealScript << "\n";
1146     RealScript << "hadd FinalResult.root higgsCombineTest.HybridNew*\n";
1147     RealScript << "\n";
1148 buchmann 1.13 RealScript << "g++ " << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/ShapeLimits/ReadAndSave.C -o ReadAndSave.exec `root-config --glibs --cflags` \n";
1149 buchmann 1.12 RealScript << "./ReadAndSave.exec FinalResult.root " << RunDirectory << "/FullCLsShapeDropletResult.txt \n";
1150     RealScript << "\n";
1151 buchmann 1.13 RealScript << "rm " << RunDirectory << "/rstat* \n";
1152     RealScript << "\n";
1153 buchmann 1.12 RealScript << "#####\n";
1154     RealScript << "echo \"Finalizing ...\"\n";
1155     RealScript << "cd $ORIGDIR\n";
1156     RealScript << "echo $ORIGDIR\n";
1157     RealScript << "ls -ltrh ../SandBox/ | grep combine\n";
1158     RealScript << "export SCRAM_ARCH=$origscramarch\n";
1159 buchmann 1.13 RealScript << "export TMP=$ORIGTMP\n";
1160     RealScript << "export TMPDIR=$ORIGTMPDIR\n";
1161 buchmann 1.12 RealScript << "exit 0\n";
1162     RealScript << "\n";
1163     RealScript.close();
1164     gSystem->Exec(("bash "+RunDirectory+"/LimitScript.sh").c_str());
1165     }
1166    
1167    
1168 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) {
1169 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
1170     if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
1171    
1172     time_t timestamp;
1173     tm *now;
1174     timestamp = time(0);
1175     now = localtime(&timestamp);
1176     stringstream RunDirectoryS;
1177     RunDirectoryS << PlottingSetup::cbafbasedir << "/exchange/ShapeComputation___" << name << "__" << now->tm_hour << "h_" << now->tm_min << "m_" << now->tm_sec << "s___" << rand();
1178     string RunDirectory = RunDirectoryS.str();
1179    
1180     ensure_directory_exists(RunDirectory);
1181    
1182 buchmann 1.14 TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str());
1183 buchmann 1.1 if(datafile->IsZombie()) {
1184     write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
1185     assert(!datafile->IsZombie());
1186     }
1187 buchmann 1.7 dout << "Run Directory: " << RunDirectory << endl;
1188 buchmann 1.2 TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
1189 buchmann 1.1
1190     TIter nextkey(datafile->GetListOfKeys());
1191     TKey *key;
1192     while ((key = (TKey*)nextkey()))
1193     {
1194     TH1F *obj =(TH1F*) key->ReadObj();
1195     limfile->cd();
1196     obj->Write();
1197     }
1198    
1199     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
1200    
1201     bool signalonly=true;
1202     bool dataonly=false;
1203    
1204     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
1205     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
1206     generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
1207    
1208     float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
1209     float scaledown=0,scaleup=0,scalesyst=0;
1210     float mceff=0,mcefferr=0;
1211     int flipped=0;
1212     int Neventsinfile=10000;//doesn't matter.
1213     string informalname="ShapeAnalysis";
1214    
1215     bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
1216    
1217     MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
1218 buchmann 1.6 if(mceff<0) {
1219     flipped=1;
1220     write_info(__FUNCTION__,"Doing flipping!");
1221     }
1222 buchmann 1.1 doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
1223 buchmann 1.4 float PDFuncert=0;
1224 buchmann 1.1 int NPdfs = get_npdfs(events);
1225 buchmann 1.13 if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
1226 buchmann 1.1
1227 buchmann 1.2 float JZBscale=scaledown;
1228     if(scaleup>scaledown) JZBscale=scaleup;
1229     dout << endl;
1230     dout << endl;
1231     dout << "Will use the following scalar (non-shape) systematics : " << endl;
1232     dout << " JZB Scale (JSU) : " << JZBscale << endl;
1233     dout << " PDF : " << PDFuncert << endl;
1234    
1235 buchmann 1.4 float SignalIntegral;
1236     if(flipped) {
1237     TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
1238     SignalIntegral= flisi ->Integral();
1239     delete flisi;
1240     }
1241     else {
1242     TH1F *flisi = (TH1F*)(limfile->Get("signal"));
1243     SignalIntegral= flisi ->Integral();
1244     delete flisi;
1245     }
1246    
1247 buchmann 1.2 TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
1248    
1249     TIter nnextkey(limfile->GetListOfKeys());
1250     TKey *nkey;
1251 buchmann 1.4
1252     if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
1253    
1254 buchmann 1.2 while ((nkey = (TKey*)nnextkey()))
1255     {
1256     TH1F *obj =(TH1F*) nkey->ReadObj();
1257     if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
1258     if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
1259     if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
1260     final_limfile->cd();
1261 buchmann 1.4 if(Contains(obj->GetName(),"signal")) {
1262     obj->Scale(1.0/SignalIntegral);
1263     }
1264 buchmann 1.2 obj->Write();
1265     }
1266 buchmann 1.1
1267 buchmann 1.6 prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
1268 buchmann 1.4
1269 buchmann 1.2 final_limfile->Close();
1270 buchmann 1.1 limfile->Close();
1271 buchmann 1.4 dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
1272     stringstream command;
1273 buchmann 1.12 string DropletLocation;
1274     int CreatedModelFileExitCode;
1275 buchmann 1.13 //----------------------------------------------------------------
1276     //do an asymptotic limit first
1277     command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1278     dout <<"Going to run : " << command.str() << endl;
1279     CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1280     dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1281     if(!(CreatedModelFileExitCode==0)) {
1282 buchmann 1.12 write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1283 buchmann 1.13 ShapeDroplet beta;
1284     beta.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1285     beta.SignalIntegral=1;
1286     return beta;
1287 buchmann 1.4 }
1288 buchmann 1.13 DropletLocation=RunDirectory+"/ShapeDropletResult.txt";
1289     ShapeDroplet beta;
1290     beta.readDroplet(DropletLocation);
1291     float obsLimit = beta.observed/SignalIntegral;
1292     dout << "Checking if the obtained limit (" << beta.observed << " / " << SignalIntegral << " = " << beta.observed/SignalIntegral << " is in [" << low_fullCLs << " , " << high_CLs << "]" << endl;
1293     if(obsLimit<high_CLs && obsLimit>low_fullCLs) {
1294     //if(PlottingSetup::scantype!=mSUGRA||(obsLimit<high_CLs && obsLimit>low_fullCLs)) {
1295     dout << " It is! Full CLs ahead!" << endl;
1296     DoFullCls(RunDirectory,beta.observed);
1297 buchmann 1.12 DropletLocation=RunDirectory+"/FullCLsShapeDropletResult.txt";
1298 buchmann 1.13 } else {
1299     dout << " It is not ... happy with asymptotic limits." << endl;
1300 buchmann 1.6 }
1301 buchmann 1.13 //----------------------------------------------------------------
1302 buchmann 1.2 ShapeDroplet alpha;
1303 buchmann 1.12 alpha.readDroplet(DropletLocation);
1304 buchmann 1.2 alpha.PDF=PDFuncert;
1305     alpha.JSU=JZBscale;
1306 buchmann 1.4 alpha.SignalIntegral=SignalIntegral;
1307 buchmann 1.2 dout << "Done reading limit information from droplet" << endl;
1308     dout << alpha << endl;
1309 buchmann 1.11
1310 buchmann 1.1 dout << "Everything is saved in " << RunDirectory << endl;
1311 buchmann 1.7 dout << "Cleaning up ... " << std::flush;
1312 buchmann 1.13 write_warning(__FUNCTION__,"NOT CLEANING UP");
1313 buchmann 1.14 // gSystem->Exec(("rm -rf "+RunDirectory).c_str());
1314 buchmann 1.4 dout << " ok!" << endl;
1315     delete limcan;
1316     return alpha;
1317 buchmann 1.1 }
1318    
1319 buchmann 1.11 void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1320 buchmann 1.14 dout << "Preparing data shapes ... " << endl;
1321 buchmann 1.1 map < pair<float, float>, map<string, float> > xsec;
1322     TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
1323     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
1324     TTree *faketree;
1325     bool dataonly=true;
1326     bool signalonly=false;
1327 buchmann 1.11
1328 buchmann 1.1 generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
1329 buchmann 1.11
1330     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","",xsec);
1331     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESup,"","JESUp",xsec);
1332     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESdown,"","JESDown",xsec);
1333     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffUp"),cutnJets,"","EffUp",xsec);
1334     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
1335     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
1336     generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
1337 buchmann 1.14
1338 buchmann 1.1 datafile->Close();
1339 buchmann 1.14
1340 buchmann 1.1 }
1341