ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/LimitCalculation.C
Revision: 1.1
Committed: Mon Jan 30 14:46:25 2012 UTC (13 years, 3 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Log Message: 
Initial commit of Ice Cream versions

File Contents

# User Rev Content
1 buchmann 1.1 /****
2    
3     Off peak status (RestrictToMassPeak) :
4    
5     x Necessary adaptations identified
6     x Started working on necessary adaptations
7     x Necessary adaptations implemented
8     x Necessary adaptations tested
9    
10     DONE!
11    
12    
13     ****/
14     #include <iostream>
15     #include <vector>
16     #include <sys/stat.h>
17     #include <fstream>
18    
19     #include <TCut.h>
20     #include <TROOT.h>
21     #include <TCanvas.h>
22     #include <TMath.h>
23     #include <TColor.h>
24     #include <TPaveText.h>
25     #include <TRandom.h>
26     #include <TH1.h>
27     #include <TH2.h>
28     #include <TF1.h>
29     #include <TSQLResult.h>
30     #include <TProfile.h>
31     #include <TSystem.h>
32     #include "LimitDroplet.C"
33    
34     //#include "TTbar_stuff.C"
35     using namespace std;
36    
37     using namespace PlottingSetup;
38    
39    
40     void rediscover_the_top(string mcjzb, string datajzb) {
41     dout << "Hi! today we are going to (try to) rediscover the top!" << endl;
42     TCanvas *c3 = new TCanvas("c3","c3");
43     c3->SetLogy(1);
44     vector<float> binning;
45     //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
46     /*
47     binning.push_back(50);
48     binning.push_back(100);
49     binning.push_back(150);
50     binning.push_back(200);
51     binning.push_back(500);
52    
53    
54     TH1F *dataprediction = allsamples.Draw("dataprediction", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
55     TH1F *puresignal = allsamples.Draw("puresignal", datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
56     // TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("TTJets"));
57     TH1F *observed = allsamples.Draw("observed", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
58     /*
59     ofstream myfile;
60     TH1F *ratio = (TH1F*)observed->Clone();
61     ratio->Divide(dataprediction);
62     ratio->GetYaxis()->SetTitle("Ratio obs/pred");
63     ratio->Draw();
64     c3->SaveAs("testratio.png");
65     myfile.open ("ShapeFit_log.txt");
66     establish_upper_limits(observed,dataprediction,puresignal,"LM4",myfile);
67     myfile.close();
68     */
69    
70    
71     int nbins=100;
72     float low=0;
73     float hi=500;
74     TCanvas *c4 = new TCanvas("c4","c4",900,900);
75     c4->Divide(2,2);
76     c4->cd(1);
77     c4->cd(1)->SetLogy(1);
78     TH1F *datapredictiont = allsamples.Draw("datapredictiont", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
79     TH1F *datapredictiono = allsamples.Draw("datapredictiono", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
80     datapredictiont->Add(datapredictiono,-1);
81     dout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
82     vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10);
83     datapredictiont->SetMarkerColor(kRed);
84     datapredictiont->SetLineColor(kRed);
85     datapredictiont->Draw();
86     functions[1]->Draw("same");
87     TText *title1 = write_title("Top Background Prediction (JZB<0, with osof subtr)");
88     title1->Draw();
89    
90     c4->cd(2);
91     c4->cd(2)->SetLogy(1);
92     TH1F *observedt = allsamples.Draw("observedt", datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
93     observedt->Draw();
94     datapredictiont->Draw("histo,same");
95     functions[1]->Draw("same");
96     TText *title2 = write_title("Observed and predicted background");
97     title2->Draw();
98    
99     c4->cd(3);
100     c4->cd(3)->SetLogy(1);
101     // TH1F *ratio = (TH1F*)observedt->Clone();
102    
103     TH1F *analytical_background_prediction= new TH1F("analytical_background_prediction","",nbins,low,hi);
104     for(int i=0;i<=nbins;i++) {
105     analytical_background_prediction->SetBinContent(i+1,functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5)));
106     analytical_background_prediction->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5))));
107     }
108     analytical_background_prediction->GetYaxis()->SetTitle("JZB [GeV]");
109     analytical_background_prediction->GetYaxis()->CenterTitle();
110     TH1F *analyticaldrawonly = (TH1F*)analytical_background_prediction->Clone();
111     analytical_background_prediction->SetFillColor(TColor::GetColor("#3399FF"));
112     analytical_background_prediction->SetMarkerSize(0);
113     analytical_background_prediction->Draw("e5");
114     analyticaldrawonly->Draw("histo,same");
115     functions[1]->Draw("same");
116     TText *title3 = write_title("Analytical bg pred histo");
117     title3->Draw();
118    
119     c4->cd(4);
120     // c4->cd(4)->SetLogy(1);
121     vector<float> ratio_binning;
122     ratio_binning.push_back(0);
123     ratio_binning.push_back(5);
124     ratio_binning.push_back(10);
125     ratio_binning.push_back(20);
126     ratio_binning.push_back(50);
127     // ratio_binning.push_back(60);
128     /*
129     ratio_binning.push_back(51);
130     ratio_binning.push_back(52);
131     ratio_binning.push_back(53);
132     ratio_binning.push_back(54);
133     ratio_binning.push_back(55);
134     ratio_binning.push_back(56);
135     ratio_binning.push_back(57);
136     ratio_binning.push_back(58);
137     ratio_binning.push_back(59);
138     ratio_binning.push_back(60);
139     // ratio_binning.push_back(70);*/
140     // ratio_binning.push_back(80);
141     // ratio_binning.push_back(90);
142     ratio_binning.push_back(80);
143     // ratio_binning.push_back(110);
144     ratio_binning.push_back(500);
145    
146     TH1F *observedtb = allsamples.Draw("observedtb", datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
147     TH1F *datapredictiontb = allsamples.Draw("datapredictiontb", "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
148     TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo", "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
149     datapredictiontb->Add(datapredictiontbo,-1);
150     TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data, luminosity);
151     for(int i=0;i<=ratio_binning.size();i++) {
152     analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
153     analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
154     }
155    
156     TH1F *ratio = (TH1F*) observedtb->Clone();
157     ratio->Divide(datapredictiontb);
158    
159     for (int i=0;i<=ratio_binning.size();i++) {
160     dout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
161     }
162    
163     // ratio->Divide(datapredictiontb);
164     // ratio->Divide(analytical_background_predictionb);
165     // TGraphAsymmErrors *JZBratio= histRatio(observedtb,analytical_background_predictionb,data,ratio_binning);
166     // ratio->Divide(analytical_background_prediction);
167     // ratio->Divide(datapredictiont);
168     // ratio->GetYaxis()->SetTitle("obs/pred");
169     // JZBratio->Draw("AP");
170     ratio->GetYaxis()->SetRangeUser(0,10);
171     ratio->Draw();
172     //analytical_background_predictionb->Draw();
173     // JZBratio->SetTitle("");
174     TText *title4 = write_title("Ratio of observed to predicted");
175     title4->Draw();
176    
177     // CompleteSave(c4,"test/ttbar_discovery_dataprediction___analytical_function");
178     CompleteSave(c4,"test/ttbar_discovery_dataprediction__analytical__new_binning_one_huge_bin_from_80");
179    
180    
181    
182    
183    
184     }
185    
186     vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected, int flipped) {
187     float sigma95=-9.9,sigma95A=-9.9;
188     /*
189     USAGE OF ROOSTATS_CL95
190     " Double_t limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
191     " LimitResult expected_limit = roostats_clm(ilum, slum, eff, seff, bck, sbck, ntoys, nuisanceModel, method, seed); \n"
192     " Double_t average_limit = roostats_cla(ilum, slum, eff, seff, bck, sbck, nuisanceModel, method, seed); \n"
193     " \n"
194     "
195     " Double_t obs_limit = limit.GetObservedLimit(); \n"
196     " Double_t exp_limit = limit.GetExpectedLimit(); \n"
197     " Double_t exp_up = limit.GetOneSigmaHighRange(); \n"
198     " Double_t exp_down = limit.GetOneSigmaLowRange(); \n"
199     " Double_t exp_2up = limit.GetTwoSigmaHighRange(); \n"
200     " Double_t exp_2down = limit.GetTwoSigmaLowRange(); \n"
201     */
202     if(mceff<=0) {
203     write_warning(__FUNCTION__,"Cannot compute upper limit in this configuration as the efficiency is negative:");
204     dout << "mc efficiency=" << mceff << " +/- " << mcefferr;
205     vector<float> sigmas;
206     sigmas.push_back(-1);
207     sigmas.push_back(-1);
208     return sigmas;
209     } else {
210     int nlimittoysused=1;
211    
212     ///------------------------------------------ < NEW > ----------------------------------------------------------
213    
214     int secondssince1970=time(NULL);
215     stringstream repname;
216     repname << PlottingSetup::cbafbasedir << "/exchange/report_" << secondssince1970 << "_"<<plotfilename<< "__"<< ".txt";
217    
218     /* - report filename [1]
219     - luminosity [2]
220     - lumi uncert [3]
221     - MC efficiency [4]
222     - MC efficiency error [5]
223     - Npred [6]
224     - Nprederr [7]
225     - Nobs [8]
226     - JZB cut [9]
227     - plot name [10]*/
228    
229     if(flipped==0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
230     if(flipped>0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << flippedNpred[ibin] << "," << flippedNprederr[ibin] << "," << flippedNobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
231    
232     stringstream command;
233     if(flipped==0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected;
234     if(flipped>0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << flippedNpred[ibin] << " " << flippedNprederr[ibin] << " " << flippedNobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected;
235     dout << command.str() << endl;
236    
237     int retval = 256;
238     int attempts=0;
239     while(!(retval==0||attempts>=3)) {//try up to 3 times
240     attempts++;
241     dout << "Starting limit calculation (TimedLimitCapsule) now : Attempt " << attempts << endl;
242     retval=gSystem->Exec(command.str().c_str());
243     }
244     char hostname[1023];
245     gethostname(hostname,1023);
246     if((!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))&&retval==256) {
247     //running via CRAB and encountered the same problem too often: place a problem file to mark this problem!
248     stringstream markproblem;
249     markproblem << "touch " << PlottingSetup::basedirectory << "/exchange/problemswhilesettinglimits.txt";
250     gSystem->Exec(markproblem.str().c_str());
251     }
252     LimitDroplet limres;
253     limres.readDroplet(repname.str());
254     dout << limres << endl;
255     remove(repname.str().c_str());
256     sigma95=limres.observed;
257    
258    
259     ///------------------------------------------ < /NEW > ----------------------------------------------------------
260     vector<float> sigmas;
261     sigmas.push_back(sigma95);
262     if(doexpected) {
263     sigmas.push_back(limres.expected);
264     sigmas.push_back(limres.upper68);
265     sigmas.push_back(limres.lower68);
266     sigmas.push_back(limres.upper95);
267     sigmas.push_back(limres.lower95);
268     }
269    
270     return sigmas;
271    
272    
273     }//end of mc efficiency is ok
274     }
275    
276     void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected, int flipped) {
277     dout << "Doing counting experiment ... " << endl;
278     vector<vector<string> > limits;
279     vector<vector<float> > vlimits;
280    
281    
282     for(int isample=0;isample<signalsamples.collection.size();isample++) {
283     vector<string> rows;
284     vector<float> vrows;
285     dout << "Considering sample " << signalsamples.collection[isample].samplename << endl;
286     rows.push_back(signalsamples.collection[isample].samplename);
287     for(int ibin=0;ibin<jzbcuts.size();ibin++) {
288     dout << "_________________________________________________________________________________" << endl;
289     float JZBcutat=uncertainties[isample*jzbcuts.size()+ibin][0];
290     float mceff=uncertainties[isample*jzbcuts.size()+ibin][1];
291     float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
292     float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
293     float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
294     float observed,observederr,null,result;
295    
296     string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
297     dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
298     vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped);
299    
300     tripple LibraryUpperLimits;
301     LibraryUpperLimits.name=signalsamples.collection[isample].samplename;
302     LibraryUpperLimits.first=mceff*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
303     LibraryUpperLimits.second=staterr*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
304     int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
305     (allresults.predictions[resultindex]).SignalYield.push_back(LibraryUpperLimits);
306    
307     if(doexpected) {
308     rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
309     vrows.push_back(sigmas[0]);
310     vrows.push_back(sigmas[1]);
311     vrows.push_back(signalsamples.collection[isample].xs);
312     }
313     else {
314     rows.push_back(any2string(sigmas[0]));
315     vrows.push_back(sigmas[0]);
316     vrows.push_back(signalsamples.collection[isample].xs);
317     }
318     }//end of bin loop
319     limits.push_back(rows);
320     vlimits.push_back(vrows);
321     }//end of sample loop
322    
323     dout << endl << endl << "_______________________________________________________________________________________" << endl;
324     dout << "Going to store upper limit on event yield in result library: " << endl;
325     for(int ibin=0;ibin<jzbcuts.size();ibin++) {
326     int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
327     vector<float> Normsigmas = compute_one_upper_limit(1.0,0.0, resultindex, mcjzb, "UPPERLIMIT", false, 0);
328     (allresults.predictions[resultindex]).UpperLimit=Normsigmas[0]*PlottingSetup::luminosity;
329     }
330     dout << endl << "_______________________________________________________________________________________" << endl;
331    
332     dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
333     dout << endl << endl << "PAS table 3: (notation: limit [95%CL])" << endl << endl;
334     dout << "\t";
335     for (int irow=0;irow<jzbcuts.size();irow++) {
336     dout << jzbcuts[irow] << "\t";
337     }
338     dout << endl;
339     for(int irow=0;irow<limits.size();irow++) {
340     for(int ientry=0;ientry<limits[irow].size();ientry++) {
341     if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
342     else dout << " (N/A) \t";
343     }
344     dout << endl;
345     }
346    
347     if(!doexpected) {
348     dout << endl << endl << "LIMITS: (Tex)" << endl;
349     tout << "\\begin{table}[hbtp]" << endl;
350     tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
351     tout << "\\begin{center}" << endl;
352     tout << "\\caption{Observed upper limits on the cross section of different LM benchmark points " << (ConsiderSignalContaminationForLimits?" (accounting for signal contamination)":" (not accounting for signal contamination)") << "}\\label{tab:lmresults}" << endl;
353     tout << "" << endl;
354     tout << "\\begin{tabular}{ | l | ";
355     for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
356     tout << "} " << endl << " \\hline " << endl << "& \t ";
357     for (int irow=0;irow<jzbcuts.size();irow++) {
358     tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
359     }
360     tout << " \\\\ \\hline " << endl;
361     for(int irow=0;irow<limits.size();irow++) {
362     tout << limits[irow][0] << " \t";
363     for(int ientry=0;ientry<jzbcuts.size();ientry++) {
364     if(vlimits[irow][2*ientry]>0) tout << " & " << Round(vlimits[irow][2*ientry],2) << " \t (" << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "x \\sigma ) \t";
365     else tout << " & ( N / A ) \t";
366     // dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
367     }
368     tout << " \\\\ \\hline " << endl;
369     }
370     tout << "\\end{tabular}" << endl;
371     tout << " \\end{tabular}"<< endl;
372     tout << "\\end{center}"<< endl;
373     tout << "\\end{table} "<< endl;
374    
375     }//do observed
376    
377     dout << endl << endl << "Final selection efficiencies with total statistical and systematic errors, and corresponding observed and expected upper limits (UL) on ($\\sigma\\times$ BR $\\times$ acceptance) for the LM4 and LM8 scenarios, in the different regions. The last column contains the predicted ($\\sigma \\times $BR$\\times$ acceptance) at NLO obtained from Monte Carlo simulation." << endl;
378     dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
379     for(int icut=0;icut<jzbcuts.size();icut++) {
380     dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?" (accounting for signal contamination)":" (not accounting for signal contamination)") << endl;
381     for(int isample=0;isample<signalsamples.collection.size();isample++) {
382     dout << limits[isample][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],3) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],3) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],3) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl;
383     }
384     dout << endl;
385     }
386     allresults.Print();
387    
388    
389     }
390    
391    
392    
393     /********************************************************************** new : Limits using SHAPES ***********************************
394    
395    
396     SSSSSSSSSSSSSSS hhhhhhh
397     SS:::::::::::::::Sh:::::h
398     S:::::SSSSSS::::::Sh:::::h
399     S:::::S SSSSSSSh:::::h
400     S:::::S h::::h hhhhh aaaaaaaaaaaaa ppppp ppppppppp eeeeeeeeeeee ssssssssss
401     S:::::S h::::hh:::::hhh a::::::::::::a p::::ppp:::::::::p ee::::::::::::ee ss::::::::::s
402     S::::SSSS h::::::::::::::hh aaaaaaaaa:::::ap:::::::::::::::::p e::::::eeeee:::::eess:::::::::::::s
403     SS::::::SSSSS h:::::::hhh::::::h a::::app::::::ppppp::::::pe::::::e e:::::es::::::ssss:::::s
404     SSS::::::::SS h::::::h h::::::h aaaaaaa:::::a p:::::p p:::::pe:::::::eeeee::::::e s:::::s ssssss
405     SSSSSS::::S h:::::h h:::::h aa::::::::::::a p:::::p p:::::pe:::::::::::::::::e s::::::s
406     S:::::S h:::::h h:::::h a::::aaaa::::::a p:::::p p:::::pe::::::eeeeeeeeeee s::::::s
407     S:::::S h:::::h h:::::ha::::a a:::::a p:::::p p::::::pe:::::::e ssssss s:::::s
408     SSSSSSS S:::::S h:::::h h:::::ha::::a a:::::a p:::::ppppp:::::::pe::::::::e s:::::ssss::::::s
409     S::::::SSSSSS:::::S h:::::h h:::::ha:::::aaaa::::::a p::::::::::::::::p e::::::::eeeeeeee s::::::::::::::s
410     S:::::::::::::::SS h:::::h h:::::h a::::::::::aa:::ap::::::::::::::pp ee:::::::::::::e s:::::::::::ss
411     SSSSSSSSSSSSSSS hhhhhhh hhhhhhh aaaaaaaaaa aaaap::::::pppppppp eeeeeeeeeeeeee sssssssssss
412     p:::::p
413     p:::::p
414     p:::::::p
415     p:::::::p
416     p:::::::p
417     ppppppppp
418    
419    
420     *********************************************************************** new : Limits using SHAPES ***********************************/
421    
422    
423     void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
424     dout << "Creatig shape templates ... ";
425     if(identifier!="") dout << "for systematic called "<<identifier;
426     dout << endl;
427     int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
428     if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
429    
430     TCut limitnJetcut;
431     if(JES==noJES) limitnJetcut=cutnJets;
432     else {
433     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
434     if(JES==JESup) limitnJetcut=cutnJetsJESup;
435     }
436     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
437     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
438     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
439     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
440    
441     TH1F *SBOSSFP;
442     TH1F *SBOSOFP;
443     TH1F *SBOSSFN;
444     TH1F *SBOSOFN;
445    
446     TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
447     TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
448     TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
449     TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
450    
451     TH1F *LSBOSSFP;
452     TH1F *LSBOSOFP;
453     TH1F *LSBOSSFN;
454     TH1F *LSBOSOFN;
455    
456     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
457     if(PlottingSetup::RestrictToMassPeak) {
458     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
459     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
460     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
461     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
462    
463     LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
464     LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
465     LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
466     LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
467     }
468    
469     string obsname="data_obs";
470     string predname="background";
471     string signalname="signal";
472     if(identifier!="") {
473     obsname=("data_"+identifier);
474     predname=("background_"+identifier);
475     signalname="signal_"+identifier;
476     }
477    
478     TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
479     obs->SetName(obsname.c_str());
480     obs->Write();
481     TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
482     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
483     if(PlottingSetup::RestrictToMassPeak) {
484     pred->Add(ZOSOFP,1.0/3);
485     pred->Add(ZOSOFN,-1.0/3);
486     pred->Add(SBOSSFP,1.0/3);
487     pred->Add(SBOSSFN,-1.0/3);
488     pred->Add(SBOSOFP,1.0/3);
489     pred->Add(SBOSOFN,-1.0/3);
490     } else {
491     pred->Add(ZOSOFP,1.0);
492     pred->Add(ZOSOFN,-1.0);
493     }
494    
495     pred->SetName(predname.c_str());
496     pred->Write();
497    
498     // TH1F *Lobs = (TH1F*)LZOSSFP->Clone();
499     // TH1F *Lpred = (TH1F*)LZOSSFN->Clone();
500    
501     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
502     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
503     Lobs->Add(LZOSSFP);
504     Lpred->Add(LZOSSFN);
505     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
506     if(PlottingSetup::RestrictToMassPeak) {
507     Lpred->Add(LZOSOFP,1.0/3);
508     Lpred->Add(LZOSOFN,-1.0/3);
509     Lpred->Add(LSBOSSFP,1.0/3);
510     Lpred->Add(LSBOSSFN,-1.0/3);
511     Lpred->Add(LSBOSOFP,1.0/3);
512     Lpred->Add(LSBOSOFN,-1.0/3);
513     } else {
514     Lpred->Add(LZOSOFP,1.0);
515     Lpred->Add(LZOSOFN,-1.0);
516     }
517    
518     TH1F *signal = (TH1F*)Lobs->Clone();
519     signal->Add(Lpred,-1);
520     signal->SetName(signalname.c_str());
521     signal->Write();
522    
523     delete Lobs;
524     delete Lpred;
525    
526     delete ZOSSFP;
527     delete ZOSOFP;
528     delete ZOSSFN;
529     delete ZOSOFN;
530    
531     if(PlottingSetup::RestrictToMassPeak) {
532     delete SBOSSFP;
533     delete SBOSOFP;
534     delete SBOSSFN;
535     delete SBOSOFN;
536     }
537    
538     delete LZOSSFP;
539     delete LZOSOFP;
540     delete LZOSSFN;
541     delete LZOSOFN;
542    
543     if(PlottingSetup::RestrictToMassPeak) {
544     delete LSBOSSFP;
545     delete LSBOSOFP;
546     delete LSBOSSFN;
547     delete LSBOSOFN;
548     }
549    
550     }
551    
552     void prepare_datacard(TFile *f) {
553     TH1F *dataob = (TH1F*)f->Get("data_obs");
554     TH1F *signal = (TH1F*)f->Get("signal");
555     TH1F *background = (TH1F*)f->Get("background");
556    
557     ofstream datacard;
558     ensure_directory_exists(get_directory()+"/limits");
559     datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
560     datacard << "Writing this to a file.\n";
561     datacard << "imax 1\n";
562     datacard << "jmax 1\n";
563     datacard << "kmax *\n";
564     datacard << "---------------\n";
565     datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
566     datacard << "---------------\n";
567     datacard << "bin 1\n";
568     datacard << "observation "<<dataob->Integral()<<"\n";
569     datacard << "------------------------------\n";
570     datacard << "bin 1 1\n";
571     datacard << "process signal background\n";
572     datacard << "process 0 1\n";
573     datacard << "rate "<<signal->Integral()<<" "<<background->Integral()<<"\n";
574     datacard << "--------------------------------\n";
575     datacard << "lumi lnN 1.10 1.0\n";
576     datacard << "bgnorm lnN 1.00 1.4 uncertainty on our prediction (40%)\n";
577     datacard << "JES shape 1 1 uncertainty on background shape and normalization\n";
578     datacard << "peak shape 1 1 uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
579     datacard << "# so divide the unit gaussian by 2 before doing the interpolation\n";
580     datacard.close();
581     }
582    
583    
584     void prepare_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
585     ensure_directory_exists(get_directory()+"/limits");
586     TFile *limfile = new TFile((get_directory()+"/limits/limitfile.root").c_str(),"RECREATE");
587     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
588     limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
589     limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
590     limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
591     limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
592     limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
593    
594     prepare_datacard(limfile);
595     limfile->Close();
596     write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
597    
598     }