1 |
#include <iostream>
|
2 |
#include <vector>
|
3 |
#include <sys/stat.h>
|
4 |
|
5 |
#include <TCut.h>
|
6 |
#include <TROOT.h>
|
7 |
#include <TCanvas.h>
|
8 |
#include <TMath.h>
|
9 |
#include <TColor.h>
|
10 |
#include <TPaveText.h>
|
11 |
#include <TRandom.h>
|
12 |
#include <TH1.h>
|
13 |
#include <TH2.h>
|
14 |
#include <TF1.h>
|
15 |
#include <TSQLResult.h>
|
16 |
#include <TProfile.h>
|
17 |
#include <TPaveStats.h>
|
18 |
|
19 |
#include "GeneratorLevelStudyModule.C"
|
20 |
|
21 |
//#include "TTbar_stuff.C"
|
22 |
using namespace std;
|
23 |
|
24 |
using namespace PlottingSetup;
|
25 |
|
26 |
|
27 |
void do_experimental_pred_obs_calculation(float cut ,string mcjzb,string datajzb, int mcordata) {
|
28 |
dout << "Crunching the numbers for JZB>" << cut << endl;
|
29 |
string xlabel="JZB [GeV] -- for algoritm internal use only!";
|
30 |
TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
|
31 |
TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
|
32 |
TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
|
33 |
TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
|
34 |
|
35 |
TH1F *SBOSSFP;
|
36 |
TH1F *SBOSOFP;
|
37 |
TH1F *SBOSSFN;
|
38 |
TH1F *SBOSOFN;
|
39 |
|
40 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
41 |
if(PlottingSetup::RestrictToMassPeak) {
|
42 |
SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
|
43 |
SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
|
44 |
SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
|
45 |
SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
|
46 |
}
|
47 |
|
48 |
|
49 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
50 |
if(PlottingSetup::RestrictToMassPeak) {
|
51 |
dout << " Observed : " << ZOSSFP->Integral() << endl;
|
52 |
dout << " Predicted: " << ZOSSFN->Integral() << " + (1/3)*(" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<") + (1/3)*(" << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<") + (1/3)*(" << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<")" << endl;
|
53 |
dout << " P(ZJets ) \t " << ZOSSFN->Integral() << endl;
|
54 |
dout << " P(eμ]) \t " << ZOSOFP->Integral() << "-" << ZOSOFN->Integral() << " = " << ZOSOFP->Integral()-ZOSOFN->Integral()<<endl;
|
55 |
dout << " P(ossf,sb]) \t " << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<" = "<<SBOSSFP->Integral()-SBOSSFN->Integral()<<endl;
|
56 |
dout << " P(osof,sb]) \t " << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<" = "<<SBOSOFP->Integral()-SBOSOFN->Integral()<<endl;
|
57 |
} else {
|
58 |
dout << " Observed : " << ZOSSFP->Integral() << endl;
|
59 |
dout << " Predicted: " << ZOSSFN->Integral() << " + (" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<")" << endl;
|
60 |
dout << " P(ZJets ) \t " << ZOSSFN->Integral() << endl;
|
61 |
dout << " P(eμ]) \t " << ZOSOFP->Integral() << "-" << ZOSOFN->Integral() << " = " << ZOSOFP->Integral()-ZOSOFN->Integral()<<endl;
|
62 |
}
|
63 |
|
64 |
|
65 |
delete ZOSSFP;
|
66 |
delete ZOSOFP;
|
67 |
delete ZOSSFN;
|
68 |
delete ZOSOFN;
|
69 |
|
70 |
delete SBOSSFP;
|
71 |
delete SBOSOFP;
|
72 |
delete SBOSSFN;
|
73 |
delete SBOSOFN;
|
74 |
}
|
75 |
|
76 |
void look_at_sidebands(string mcjzb, string datajzb, bool includejetcut, float cutat=0) {
|
77 |
|
78 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak -- this funciton is meaningless for the offpeak case
|
79 |
if(!PlottingSetup::RestrictToMassPeak) return;
|
80 |
dout << "Looking at sidebands ... " << endl;
|
81 |
int mcordata=data;//data // you can perform this study for mc or data ...
|
82 |
|
83 |
TCut specialjetcut;
|
84 |
if(includejetcut) specialjetcut=cutnJets&&basiccut;
|
85 |
else specialjetcut="mll>0";
|
86 |
dout << "The datajzb variable is defined as " << datajzb << endl;
|
87 |
stringstream addcut;
|
88 |
addcut<<"(pt>"<<cutat<<")";
|
89 |
TCut additionalcut=addcut.str().c_str();
|
90 |
|
91 |
int nbins=75; float min=51;float max=201;string xlabel="mll [GeV]";
|
92 |
TCanvas *c1 = new TCanvas("c1","c1");
|
93 |
c1->SetLogy(1);
|
94 |
TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&specialjetcut&&additionalcut,data,luminosity);
|
95 |
THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&specialjetcut&&additionalcut,mc,luminosity);
|
96 |
|
97 |
datahistoOSSF->SetMinimum(1);
|
98 |
datahistoOSSF->Draw();
|
99 |
mcstackOSSF.Draw("same");
|
100 |
datahistoOSSF->Draw("same");
|
101 |
TLegend *kinleg = allsamples.allbglegend();
|
102 |
kinleg->AddEntry(datahistoOSSF,"OSSF (data)","p");
|
103 |
kinleg->Draw();
|
104 |
if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSSF");
|
105 |
else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSSF_nojetcut");
|
106 |
|
107 |
TH1F *datahistoOSOF = allsamples.Draw("datahistoOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&specialjetcut&&additionalcut,data,luminosity);
|
108 |
THStack mcstackOSOF = allsamples.DrawStack("mcstackOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&specialjetcut&&additionalcut,mc,luminosity);
|
109 |
// datahistoOSOF->SetMinimum(0.4);
|
110 |
datahistoOSOF->Draw();
|
111 |
mcstackOSOF.Draw("same");
|
112 |
datahistoOSOF->Draw("same");
|
113 |
TLegend *kinleg2 = allsamples.allbglegend();
|
114 |
kinleg2->AddEntry(datahistoOSOF,"OSOF (data)","p");
|
115 |
kinleg2->Draw();
|
116 |
if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSOF");
|
117 |
else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSOF_nojetcut");
|
118 |
|
119 |
TH1F *rawmlleemmData = allsamples.Draw("rawmlleemmData","mll",200,0,200, "mll [GeV]", "events", cutOSSF&&specialjetcut&&sidebandcut&&additionalcut,data, luminosity);
|
120 |
TH1F *rawmllemData = allsamples.Draw("rawmllemData" ,"mll",200,0,200, "mll [GeV]", "events", cutmass&&cutOSOF&&specialjetcut&&additionalcut,data, luminosity);
|
121 |
dout << "Number of events in peak for OSOF: " << rawmllemData->GetEntries() << endl;
|
122 |
dout << "Number of events in SB for OSSF: " << rawmlleemmData->GetEntries() << endl;
|
123 |
|
124 |
TH1F *SFttbarZpeak = allsamples.Draw("SFttbarZpeak",mcjzb,100,-200,400, "JZB [GeV]", "events",cutmass&&cutOSSF&&specialjetcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
|
125 |
TH1F *OFttbarZpeak = allsamples.Draw("OFttbarZpeak",mcjzb,100,-200,400, "JZB [GeV]", "events",cutmass&&cutOSOF&&specialjetcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
|
126 |
TH1F *SFttbarsideb = allsamples.Draw("SFttbarsideb",mcjzb,100,-200,400, "JZB [GeV]", "events",cutOSSF&&specialjetcut&&sidebandcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
|
127 |
|
128 |
SFttbarZpeak->SetLineColor(kBlack);
|
129 |
OFttbarZpeak->SetLineColor(kBlue);
|
130 |
OFttbarZpeak->SetMarkerColor(kBlue);
|
131 |
SFttbarsideb->SetLineColor(kPink);
|
132 |
SFttbarsideb->SetMarkerColor(kPink);
|
133 |
|
134 |
SFttbarZpeak->Draw("histo");
|
135 |
OFttbarZpeak->Draw("histo,same");
|
136 |
SFttbarsideb->Draw("histo,same");
|
137 |
|
138 |
TLegend *leg3 = new TLegend(0.6,0.8,0.89,0.89);
|
139 |
leg3->AddEntry(SFttbarZpeak,"SF ttbar Z peak","l");
|
140 |
leg3->AddEntry(OFttbarZpeak,"OF ttbar Z peak","l");
|
141 |
leg3->AddEntry(SFttbarsideb,"SF ttbar SB","l");
|
142 |
leg3->SetFillColor(kWhite);
|
143 |
leg3->SetLineColor(kWhite);
|
144 |
leg3->SetBorderSize(0);
|
145 |
|
146 |
leg3->Draw();
|
147 |
if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison");
|
148 |
else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_nojetcut");
|
149 |
|
150 |
|
151 |
c1->SetLogy(0);
|
152 |
|
153 |
SFttbarsideb->SetFillColor(TColor::GetColor("#F5A9A9"));
|
154 |
SFttbarsideb->SetLineColor(TColor::GetColor("#F5A9A9"));
|
155 |
SFttbarsideb->SetFillStyle(3004);
|
156 |
OFttbarZpeak->SetFillColor(TColor::GetColor("#819FF7"));
|
157 |
OFttbarZpeak->SetLineColor(TColor::GetColor("#819FF7"));
|
158 |
OFttbarZpeak->SetFillColor(kBlue);
|
159 |
OFttbarZpeak->SetFillStyle(3005);
|
160 |
|
161 |
OFttbarZpeak->Rebin(2);
|
162 |
SFttbarZpeak->Rebin(2);
|
163 |
SFttbarsideb->Rebin(2);
|
164 |
OFttbarZpeak->Divide(SFttbarZpeak);
|
165 |
SFttbarsideb->Divide(SFttbarZpeak);
|
166 |
OFttbarZpeak->GetYaxis()->SetRangeUser(0,5);
|
167 |
OFttbarZpeak->GetYaxis()->SetTitle("ratio");
|
168 |
TF1 *centralfitO = new TF1("centralfitO","pol1",-40,120);
|
169 |
TF1 *centralfit1 = new TF1("centralfit1","pol1",-200,400);
|
170 |
// TF1 *centralfitS = new TF1("centralfitS","pol1",-40,120);
|
171 |
SFttbarsideb->Fit(centralfitO,"R");
|
172 |
//OFttbarZpeak->Fit(centralfitO,"R");
|
173 |
centralfit1->SetParameters(centralfitO->GetParameters());
|
174 |
// SFttbarZpeak->Fit(centralfitS,"R");
|
175 |
OFttbarZpeak->Draw("e5");
|
176 |
SFttbarsideb->Draw("e5,same");
|
177 |
OFttbarZpeak->Draw("same");
|
178 |
SFttbarsideb->Draw("same");
|
179 |
centralfit1->SetLineColor(kOrange);
|
180 |
// centralfitS->SetLineColor(kOrange);
|
181 |
centralfit1->SetLineWidth(2);
|
182 |
// centralfitS->SetLineWidth(2);
|
183 |
centralfit1->Draw("same");
|
184 |
// centralfitS->Draw("same");
|
185 |
TLine *oneline = new TLine(-200,1,400,1);
|
186 |
oneline->SetLineColor(kBlue);
|
187 |
TLine *point5 = new TLine(-200,0.5,400,0.5);
|
188 |
point5->SetLineStyle(2);
|
189 |
point5->SetLineColor(kGreen);
|
190 |
TLine *op5 = new TLine(-200,1.5,400,1.5);
|
191 |
op5->SetLineStyle(2);
|
192 |
op5->SetLineColor(kGreen);
|
193 |
TLine *point7 = new TLine(-200,0.7,400,0.7);
|
194 |
point7->SetLineStyle(2);
|
195 |
point7->SetLineColor(kBlack);
|
196 |
TLine *op7 = new TLine(-200,1.3,400,1.3);
|
197 |
op7->SetLineStyle(2);
|
198 |
op7->SetLineColor(kBlack);
|
199 |
oneline->Draw("same");
|
200 |
point5->Draw("same");
|
201 |
point7->Draw("same");
|
202 |
op5->Draw("same");
|
203 |
op7->Draw("same");
|
204 |
TLegend *leg4 = new TLegend(0.6,0.65,0.89,0.89);
|
205 |
leg4->AddEntry(OFttbarZpeak,"OF ttbar Z peak / truth","l");
|
206 |
leg4->AddEntry(SFttbarsideb,"SF ttbar SB / truth","l");
|
207 |
leg4->AddEntry(centralfit1,"Fit to [-40,120] GeV region (OF)","l");
|
208 |
leg4->AddEntry(point5,"50% systematic envelope","l");
|
209 |
leg4->AddEntry(point7,"30% systematic envelope","l");
|
210 |
// leg4->AddEntry(centralfitS,"Fit to [-40,120] GeV region (SF)","l");
|
211 |
leg4->SetFillColor(kWhite);
|
212 |
leg4->SetLineColor(kWhite);
|
213 |
leg4->SetBorderSize(0);
|
214 |
leg4->Draw("same");
|
215 |
if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_ratio");
|
216 |
else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_ratio_nojetcut");
|
217 |
dout << "Moving on to predicted / observed yields! " << endl;
|
218 |
dout << "Sideband definition: " << (const char*) sidebandcut << endl;
|
219 |
/*
|
220 |
do_experimental_pred_obs_calculation(50,mcjzb,datajzb,mcordata);
|
221 |
do_experimental_pred_obs_calculation(75,mcjzb,datajzb,mcordata);
|
222 |
do_experimental_pred_obs_calculation(100,mcjzb,datajzb,mcordata);
|
223 |
do_experimental_pred_obs_calculation(125,mcjzb,datajzb,mcordata);
|
224 |
do_experimental_pred_obs_calculation(150,mcjzb,datajzb,mcordata);
|
225 |
*/
|
226 |
|
227 |
delete rawmlleemmData;
|
228 |
delete rawmllemData;
|
229 |
delete SFttbarZpeak;
|
230 |
delete OFttbarZpeak;
|
231 |
delete SFttbarsideb;
|
232 |
delete datahistoOSOF;
|
233 |
delete datahistoOSSF;
|
234 |
|
235 |
|
236 |
}
|
237 |
|
238 |
void look_at_sidebands(string mcjzb, string datajzb) {
|
239 |
// for (int i=0;i<100;i+=10) {
|
240 |
int i=0;
|
241 |
{
|
242 |
look_at_sidebands(mcjzb,datajzb, true,i);
|
243 |
look_at_sidebands(mcjzb,datajzb, false,i);
|
244 |
}
|
245 |
}
|
246 |
|
247 |
void find_sideband_definition() {
|
248 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
249 |
if(!PlottingSetup::RestrictToMassPeak) return; // this function is meaningless for the offpeak analysis
|
250 |
|
251 |
TH1F *mllttbar = allsamples.Draw("mllttbar","mll",145,55,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&!cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
|
252 |
TH1F *mllttbarz = allsamples.Draw("mllttbarz","mll",1,50,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
|
253 |
float leftstop=0;
|
254 |
float rightstop=0;
|
255 |
int pos90=mllttbar->FindBin(90);
|
256 |
float leftsum=0; float rightsum=0;
|
257 |
float peaksum=mllttbarz->Integral();
|
258 |
for(int i=0;i<mllttbar->GetNbinsX()&&!(leftstop&&rightstop);i++) {
|
259 |
if(pos90-i<1) leftstop=mllttbar->GetBinLowEdge(1);
|
260 |
if(mllttbar->GetBinLowEdge(pos90+i)+mllttbar->GetBinWidth(pos90+i)>190) rightstop=190;
|
261 |
if(!leftstop) leftsum+=mllttbar->GetBinContent(pos90-i);
|
262 |
if(!rightstop) rightsum+=mllttbar->GetBinContent(pos90+i);
|
263 |
if(leftsum+rightsum>peaksum) {
|
264 |
if(!leftstop) leftstop=mllttbar->GetBinLowEdge(pos90-i);
|
265 |
if(!rightstop) rightstop=mllttbar->GetBinLowEdge(pos90+i)+mllttbar->GetBinWidth(pos90+i);
|
266 |
dout << "Found the boundaries! on the left: " << leftstop << " and on the right " << rightstop << endl;
|
267 |
dout << "Total sum : " << leftsum+rightsum << " which supposedly corresponds ~ to " << peaksum << endl;
|
268 |
}
|
269 |
}
|
270 |
TH1F *mllttbart = allsamples.Draw("mllttbart","mll",1,55,155, "mll [GeV]", "events",cutOSSF&&cutnJets&&!cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
|
271 |
dout << mllttbart->Integral() << endl;
|
272 |
|
273 |
|
274 |
}
|
275 |
|
276 |
void calculate_upper_limits(string mcjzb, string datajzb) {
|
277 |
write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
|
278 |
// write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
|
279 |
// rediscover_the_top(mcjzb,datajzb);
|
280 |
/*
|
281 |
TCanvas *c3 = new TCanvas("c3","c3");
|
282 |
c3->SetLogy(1);
|
283 |
vector<float> binning;
|
284 |
//binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
|
285 |
binning.push_back(50);
|
286 |
binning.push_back(100);
|
287 |
binning.push_back(150);
|
288 |
binning.push_back(200);
|
289 |
binning.push_back(500);
|
290 |
TH1F *datapredictiona = allsamples.Draw("datapredictiona", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
|
291 |
TH1F *datapredictionb = allsamples.Draw("datapredictionb", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
|
292 |
TH1F *datapredictionc = allsamples.Draw("datapredictionc", datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
|
293 |
TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
|
294 |
dataprediction->Add(datapredictionb,-1);
|
295 |
dataprediction->Add(datapredictionc);
|
296 |
TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
|
297 |
TH1F *signalpred = allsamples.Draw("signalpred", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
|
298 |
TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
|
299 |
TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
|
300 |
TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
|
301 |
signalpred->Add(signalpredlo,-1);
|
302 |
signalpred->Add(signalpredro);
|
303 |
puresignal->Add(signalpred,-1);//subtracting signal contamination
|
304 |
ofstream myfile;
|
305 |
myfile.open ("ShapeFit_log.txt");
|
306 |
establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
|
307 |
myfile.close();
|
308 |
*/
|
309 |
}
|
310 |
|
311 |
TH1F *runcheckhisto(string cut) {
|
312 |
string histoname=GetNumericHistoName();
|
313 |
TH1F *histo = new TH1F(histoname.c_str(),histoname.c_str(),100,163000,168000);
|
314 |
(allsamples.collection)[0].events->Draw(("runNum>>"+histoname).c_str(),cut.c_str());
|
315 |
return histo;
|
316 |
}
|
317 |
|
318 |
void run_check() {
|
319 |
gROOT->SetStyle("Plain");
|
320 |
TCanvas *c1 = new TCanvas("runnum","runnum",800,1000);
|
321 |
c1->Divide(2,4);
|
322 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
323 |
if(PlottingSetup::RestrictToMassPeak) c1->Divide(2,2); // there are only four regions ...
|
324 |
|
325 |
c1->cd(1);
|
326 |
TH1F *ossfp = runcheckhisto((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
|
327 |
ossfp->Draw();
|
328 |
TText *t1 = write_title("OSSF,P");t1->Draw();
|
329 |
|
330 |
c1->cd(2);
|
331 |
TH1F *ossfn = runcheckhisto((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
|
332 |
ossfn->Draw();
|
333 |
TText *t2 = write_title("OSSF,N");t2->Draw();
|
334 |
|
335 |
c1->cd(3);
|
336 |
TH1F *osofp = runcheckhisto((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
|
337 |
osofp->Draw();
|
338 |
TText *t3 = write_title("OSOF,P");t3->Draw();
|
339 |
c1->cd(4);
|
340 |
TH1F *osofn = runcheckhisto((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
|
341 |
osofn->Draw();
|
342 |
TText *t4 = write_title("OSOF,N");t4->Draw();
|
343 |
|
344 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
345 |
if(PlottingSetup::RestrictToMassPeak) {
|
346 |
c1->cd(5);
|
347 |
TH1F *sbofp = runcheckhisto((const char*)(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
|
348 |
sbofp->Draw();
|
349 |
TText *t5 = write_title("SB,OSOF,P");t5->Draw();
|
350 |
c1->cd(6);
|
351 |
TH1F *sbofn = runcheckhisto((const char*)(cutOSOF&&cutnJets&&basiccut&&sidebandcut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
|
352 |
sbofn->Draw();
|
353 |
TText *t6 = write_title("SB,OSOF,N");t6->Draw();
|
354 |
|
355 |
c1->cd(7);
|
356 |
TH1F *sbsfp = runcheckhisto((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
|
357 |
sbsfp->Draw();
|
358 |
TText *t7 = write_title("SB,OSSF,P");t7->Draw();
|
359 |
c1->cd(8);
|
360 |
TH1F *sbsfn = runcheckhisto((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
|
361 |
sbsfn->Draw();
|
362 |
TText *t8 = write_title("SB,OSSF,N");t8->Draw();
|
363 |
}
|
364 |
|
365 |
c1->SaveAs("runNumber.png");
|
366 |
}
|
367 |
|
368 |
TH1F *give_boson_pred(TCut bcut,string mcjzb) {
|
369 |
int nbins=50;
|
370 |
TH1F *jzbn = allsamples.Draw("jzbn","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
371 |
TH1F *jzbno = allsamples.Draw("jzbno","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSOF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
372 |
TH1F *jzbpo = allsamples.Draw("jzbp",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSOF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
373 |
|
374 |
//Sidebands
|
375 |
TH1F *jzbnos;
|
376 |
TH1F *jzbpos;
|
377 |
TH1F *jzbnss;
|
378 |
TH1F *jzbpss;
|
379 |
|
380 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
381 |
if(PlottingSetup::RestrictToMassPeak) {
|
382 |
jzbnos = allsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
383 |
jzbpos = allsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
384 |
jzbnss = allsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
385 |
jzbpss = allsamples.Draw("jzbpss",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
386 |
}
|
387 |
|
388 |
|
389 |
TH1F *pred = (TH1F*)jzbn->Clone("pred");
|
390 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
391 |
if(PlottingSetup::RestrictToMassPeak) {
|
392 |
pred->Add(jzbno,-1.0/3);
|
393 |
pred->Add(jzbpo,1.0/3);
|
394 |
pred->Add(jzbnos,-1.0/3);
|
395 |
pred->Add(jzbpos,1.0/3);
|
396 |
pred->Add(jzbnss,-1.0/3);
|
397 |
pred->Add(jzbpss,1.0/3);
|
398 |
} else {
|
399 |
pred->Add(jzbno,-1.0);
|
400 |
pred->Add(jzbpo,1.0);
|
401 |
}
|
402 |
pred->SetLineColor(kRed);
|
403 |
pred->SetMinimum(0);
|
404 |
delete jzbn;
|
405 |
delete jzbpo;
|
406 |
delete jzbno;
|
407 |
delete jzbpos;
|
408 |
delete jzbnos;
|
409 |
delete jzbpss;
|
410 |
delete jzbnss;
|
411 |
|
412 |
return pred;
|
413 |
}
|
414 |
|
415 |
|
416 |
void show_dibosons(string datajzb, string mcjzb) {
|
417 |
TCut WW("(abs(genMID1)==24&&abs(genMID2)==24)||(abs(genGMID1)==24&&abs(genGMID2)==24)");
|
418 |
TCut ZZ("(abs(genMID1)==23&&abs(genMID2)==23)||(abs(genGMID1)==23&&abs(genGMID2)==23)");
|
419 |
TCut WZ("((abs(genMID1)==23&&abs(genMID2)==24)||(abs(genGMID1)==23&&abs(genGMID2)==24))||((abs(genMID1)==24&&abs(genMID2)==23)||(abs(genGMID1)==24&&abs(genGMID2)==23))");
|
420 |
|
421 |
TCanvas *dibs = new TCanvas("dibs","dibs",900,900);
|
422 |
dibs->Divide(2,2);
|
423 |
|
424 |
dibs->cd(1);
|
425 |
TH1F *wwjzbp = allsamples.Draw("wwjzbp",mcjzb,70,0,350, "JZB [GeV]", "events", WW&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
426 |
TH1F *wwpred = (TH1F*) give_boson_pred(WW,mcjzb);
|
427 |
wwpred->Draw("histo");
|
428 |
wwjzbp->Draw("histo,same");
|
429 |
TLegend *leg = make_legend("WW");
|
430 |
leg->SetFillColor(kWhite);
|
431 |
leg->SetLineColor(kWhite);
|
432 |
leg->SetHeader("WW (MC)");
|
433 |
leg->AddEntry(wwjzbp,"Observed","l");
|
434 |
leg->AddEntry(wwpred,"Predicted","l");
|
435 |
leg->Draw("same");
|
436 |
|
437 |
dibs->cd(2);
|
438 |
TH1F *wzjzbp = allsamples.Draw("wzjzbp",mcjzb,70,0,350, "JZB [GeV]", "events",WZ&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
439 |
TH1F *wzpred = (TH1F*) give_boson_pred(WZ,mcjzb);
|
440 |
wzpred->Draw("histo");
|
441 |
wzjzbp->Draw("same,histo");
|
442 |
TLegend *leg2 = (TLegend*)leg->Clone("leg2");
|
443 |
leg2->SetHeader("WZ (MC)");
|
444 |
leg2->Draw("same");
|
445 |
DrawPrelim();
|
446 |
|
447 |
dibs->cd(3);
|
448 |
TH1F *zzjzbp = allsamples.Draw("zzjzbp",mcjzb,70,0,350, "JZB [GeV]", "events",ZZ&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
449 |
TH1F *zzpred = (TH1F*) give_boson_pred(ZZ,mcjzb);
|
450 |
zzpred->Draw("histo");
|
451 |
zzjzbp->Draw("same,histo");
|
452 |
TLegend *leg3 = (TLegend*)leg->Clone("leg2");
|
453 |
leg3->SetHeader("ZZ (MC)");
|
454 |
leg3->Draw("same");
|
455 |
leg3->Draw("same");
|
456 |
DrawPrelim();
|
457 |
|
458 |
dibs->cd(4);
|
459 |
TH1F *alljzbp = allsamples.Draw("alljzbp",mcjzb,70,0,350, "JZB [GeV]", "events",(WW||WZ||ZZ)&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
|
460 |
TH1F *allpred = (TH1F*) give_boson_pred((WW||WZ||ZZ),mcjzb);
|
461 |
allpred->Draw("histo");
|
462 |
alljzbp->Draw("same,histo");
|
463 |
TLegend *leg4 = (TLegend*)leg->Clone("leg2");
|
464 |
leg4->SetHeader("All dibosons (MC)");
|
465 |
leg4->Draw("same");
|
466 |
DrawPrelim();
|
467 |
|
468 |
CompleteSave(dibs,"Studies/Dibosons");
|
469 |
|
470 |
}
|
471 |
|
472 |
void show_rare_samples(string datajzb, string mcjzb) {
|
473 |
TCanvas *rares = new TCanvas("rares","Rare Samples",900,900);
|
474 |
|
475 |
int nbins=50;
|
476 |
TH1F *jzbp = raresample.Draw("jzbp", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity);
|
477 |
|
478 |
TH1F *jzbn = raresample.Draw("jzbn", "-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity);
|
479 |
TH1F *jzbno = raresample.Draw("jzbno", "-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity);
|
480 |
TH1F *jzbpo = raresample.Draw("jzbp", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity);
|
481 |
|
482 |
//Sidebands
|
483 |
TH1F *jzbnos = raresample.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity);
|
484 |
TH1F *jzbpos = raresample.Draw("jzbpos", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity);
|
485 |
TH1F *jzbnss = raresample.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity);
|
486 |
TH1F *jzbpss = raresample.Draw("jzbpss", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity);
|
487 |
|
488 |
TH1F *pred = (TH1F*)jzbn->Clone("pred");
|
489 |
pred->Add(jzbno,-1.0/3);
|
490 |
pred->Add(jzbpo,1.0/3);
|
491 |
pred->Add(jzbnos,-1.0/3);
|
492 |
pred->Add(jzbpos,1.0/3);
|
493 |
pred->Add(jzbnss,-1.0/3);
|
494 |
pred->Add(jzbpss,1.0/3);
|
495 |
pred->SetLineColor(kRed);
|
496 |
pred->SetMinimum(0);
|
497 |
delete jzbn;
|
498 |
delete jzbpo;
|
499 |
delete jzbno;
|
500 |
delete jzbpos;
|
501 |
delete jzbnos;
|
502 |
delete jzbpss;
|
503 |
delete jzbnss;
|
504 |
|
505 |
pred->Draw("histo");
|
506 |
pred->GetYaxis()->SetTitleOffset(1.3);
|
507 |
jzbp->Draw("histo,same");
|
508 |
TLegend *leg = make_legend("WW");
|
509 |
leg->SetFillColor(kWhite);
|
510 |
leg->SetLineColor(kWhite);
|
511 |
leg->SetHeader("Rare Samples (MC)");
|
512 |
leg->AddEntry(jzbp,"Observed","l");
|
513 |
leg->AddEntry(pred,"Predicted","l");
|
514 |
leg->Draw("same");
|
515 |
|
516 |
|
517 |
CompleteSave(rares,"Studies/Rare_Samples");
|
518 |
delete jzbp;
|
519 |
delete pred;
|
520 |
delete leg;
|
521 |
delete rares;
|
522 |
}
|
523 |
|
524 |
|
525 |
class signature {
|
526 |
public:
|
527 |
int runNum;
|
528 |
int eventNum;
|
529 |
int lumi;
|
530 |
};
|
531 |
|
532 |
vector<signature> get_list_of_events(string cut) {
|
533 |
float jzb;
|
534 |
int runNum,lumi,eventNum;
|
535 |
(allsamples.collection)[0].events->SetBranchAddress("jzb",&jzb);
|
536 |
(allsamples.collection)[0].events->SetBranchAddress("eventNum",&eventNum);
|
537 |
(allsamples.collection)[0].events->SetBranchAddress("lumi",&lumi);
|
538 |
(allsamples.collection)[0].events->SetBranchAddress("runNum",&runNum);
|
539 |
|
540 |
TTreeFormula *select = new TTreeFormula("select", cut.c_str()&&essentialcut, (allsamples.collection)[0].events);
|
541 |
vector<signature> allevents;
|
542 |
for (Int_t entry = 0 ; entry < (allsamples.collection)[0].events->GetEntries() ; entry++) {
|
543 |
(allsamples.collection)[0].events->LoadTree(entry);
|
544 |
if (select->EvalInstance()) {
|
545 |
(allsamples.collection)[0].events->GetEntry(entry);
|
546 |
signature newevent;
|
547 |
newevent.runNum=runNum;
|
548 |
newevent.eventNum=eventNum;
|
549 |
newevent.lumi=lumi;
|
550 |
allevents.push_back(newevent);
|
551 |
}
|
552 |
}
|
553 |
cout << "Done looping!" << endl;
|
554 |
return allevents;
|
555 |
}
|
556 |
|
557 |
void make_double_plot(string variable, int nbins, float min, float max, float ymax, bool logscale,string xlabel, string filename,TCut observed, TCut predicted, bool is_data, bool noscale = false) {
|
558 |
TCut ibasiccut=basiccut;
|
559 |
|
560 |
//Step 2: Refine the cut
|
561 |
TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
|
562 |
ckin->SetLogy(logscale);
|
563 |
|
564 |
TH1F *datahistoa = allsamples.Draw("datahistoa",variable,nbins,min,max, xlabel, "events",observed,is_data,luminosity);
|
565 |
TH1F *datahistob = allsamples.Draw("datahistob",variable,nbins,min,max, xlabel, "events",predicted,is_data,luminosity);
|
566 |
|
567 |
datahistoa->SetLineColor(kBlue);
|
568 |
datahistob->SetLineColor(kRed);
|
569 |
|
570 |
TLegend *kinleg = new TLegend(0.6,0.7,0.8,0.89);
|
571 |
kinleg->SetFillColor(kWhite);
|
572 |
kinleg->SetLineColor(kWhite);
|
573 |
kinleg->SetBorderSize(0);
|
574 |
kinleg->AddEntry(datahistoa,"Observed "+TString(is_data?"Data":"MC"),"l");
|
575 |
kinleg->AddEntry(datahistob,"Predicted "+TString(is_data?"Data":"MC"),"l");
|
576 |
|
577 |
datahistoa->SetMaximum(ymax);
|
578 |
datahistoa->Draw("histo");
|
579 |
datahistob->SetLineStyle(2);
|
580 |
if ( !noscale ) datahistob->Scale(0.3);
|
581 |
datahistob->Draw("histo,sames");
|
582 |
TVirtualPad::Pad()->Update();
|
583 |
|
584 |
TPaveStats *sa = (TPaveStats*)datahistoa->GetListOfFunctions()->FindObject("stats");
|
585 |
TPaveStats *sb = (TPaveStats*)datahistob->GetListOfFunctions()->FindObject("stats");
|
586 |
if ( sa && sb ) {
|
587 |
sa->SetTextColor(datahistoa->GetLineColor());
|
588 |
sb->SetTextColor(datahistob->GetLineColor());
|
589 |
sb->SetY1NDC(sb->GetY1NDC()-0.25);
|
590 |
sb->SetY2NDC(sb->GetY2NDC()-0.25);
|
591 |
TVirtualPad::Pad()->Update();
|
592 |
}
|
593 |
kinleg->Draw();
|
594 |
TText* write_cut = write_title(variable);
|
595 |
write_cut->Draw();
|
596 |
CompleteSave(ckin,"special_kin/"+filename);
|
597 |
datahistoa->Delete();
|
598 |
datahistob->Delete();
|
599 |
delete ckin;
|
600 |
}
|
601 |
|
602 |
TH1F *give_lm0_pred(TCut bcut,string mcjzb) {
|
603 |
int nbins=50;
|
604 |
TH1F *jzbn = signalsamples.Draw("jzbn","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
|
605 |
TH1F *jzbno = signalsamples.Draw("jzbno","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
|
606 |
TH1F *jzbpo = signalsamples.Draw("jzbp",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSOF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
|
607 |
|
608 |
//Sidebands
|
609 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
610 |
TH1F *jzbnos;
|
611 |
TH1F *jzbpos;
|
612 |
TH1F *jzbnss;
|
613 |
TH1F *jzbpss;
|
614 |
|
615 |
if(PlottingSetup::RestrictToMassPeak) {
|
616 |
jzbnos = signalsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
|
617 |
jzbpos = signalsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
|
618 |
jzbnss = signalsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
|
619 |
jzbpss = signalsamples.Draw("jzbpss",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
|
620 |
}
|
621 |
|
622 |
TH1F *pred = (TH1F*)jzbn->Clone("pred");
|
623 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
624 |
if(PlottingSetup::RestrictToMassPeak) {
|
625 |
pred->Add(jzbno,-1.0/3);
|
626 |
pred->Add(jzbpo,1.0/3);
|
627 |
pred->Add(jzbnos,-1.0/3);
|
628 |
pred->Add(jzbpos,1.0/3);
|
629 |
pred->Add(jzbnss,-1.0/3);
|
630 |
pred->Add(jzbpss,1.0/3);
|
631 |
} else {
|
632 |
pred->Add(jzbno,-1.0);
|
633 |
pred->Add(jzbpo,1.0);
|
634 |
}
|
635 |
|
636 |
pred->SetLineColor(kRed);
|
637 |
pred->SetMinimum(0);
|
638 |
|
639 |
delete jzbn;
|
640 |
delete jzbpo;
|
641 |
delete jzbno;
|
642 |
delete jzbpos;
|
643 |
delete jzbnos;
|
644 |
delete jzbpss;
|
645 |
delete jzbnss;
|
646 |
cout << "pred contains " << pred->Integral() << "entries" << endl;
|
647 |
return pred;
|
648 |
}
|
649 |
|
650 |
|
651 |
void lm0_illustration() {
|
652 |
TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
|
653 |
can->SetLogy(1);
|
654 |
|
655 |
int sbg_nbins=130;
|
656 |
float sbg_min=-500; //-110;
|
657 |
float sbg_max=800; //jzbHigh;
|
658 |
|
659 |
float simulatedlumi=1000;//in pb please
|
660 |
|
661 |
TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("DYJetsToLL"));
|
662 |
TH1F *JZBplotLM0 = signalsamples.Draw("JZBplotLM0",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,signalsamples.FindSample("LM0"));
|
663 |
TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("TTJets"));
|
664 |
|
665 |
JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
|
666 |
JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
|
667 |
JZBplotZJETs->SetLineColor(kBlack);
|
668 |
JZBplotLM0->SetLineStyle(2);
|
669 |
JZBplotZJETs->SetMaximum(JZBplotZJETs->GetMaximum()*5);
|
670 |
JZBplotZJETs->SetMinimum(1);
|
671 |
|
672 |
JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
|
673 |
JZBplotTtbar->SetMinimum(0.01);
|
674 |
JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
|
675 |
JZBplotTtbar->DrawClone("histo");
|
676 |
JZBplotZJETs->Draw("histo,same");
|
677 |
JZBplotTtbar->SetFillColor(0);
|
678 |
JZBplotTtbar->DrawClone("histo,same");
|
679 |
JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
|
680 |
JZBplotLM0->Draw("histo,same");
|
681 |
|
682 |
|
683 |
TLegend *signal_bg_comparison_leg2 = make_legend("",0.55,0.75,false);
|
684 |
signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
|
685 |
signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
|
686 |
signal_bg_comparison_leg2->AddEntry(JZBplotLM0,"LM0","f");
|
687 |
signal_bg_comparison_leg2->Draw();
|
688 |
TText* title = write_title("CMS MC simulation, #sqrt{s}= 7 TeV @ L=1 fb^{-1}");
|
689 |
title->Draw();
|
690 |
CompleteSave(can,"LM0/jzb_bg_vs_signal_distribution__LM0");
|
691 |
can->SetLogy(0);
|
692 |
TH1F *lm0jzbp = signalsamples.Draw("lm0jzb",jzbvariablemc,70,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
|
693 |
TH1F *lm0pred = (TH1F*) give_lm0_pred("mll>0",jzbvariablemc);
|
694 |
lm0pred->Draw("histo");
|
695 |
lm0jzbp->Draw("histo,same");
|
696 |
TLegend *leg = make_legend("LM0");
|
697 |
leg->SetFillColor(kWhite);
|
698 |
leg->SetLineColor(kWhite);
|
699 |
leg->SetHeader("LM0 (MC)");
|
700 |
leg->AddEntry(lm0jzbp,"Observed","l");
|
701 |
leg->AddEntry(lm0pred,"Predicted","l");
|
702 |
leg->Draw("same");
|
703 |
CompleteSave(can,"LM0/LeftRight");
|
704 |
|
705 |
delete lm0jzbp;
|
706 |
delete lm0pred;
|
707 |
}
|
708 |
|
709 |
void kinematic_dist_of_pred_and_obs() {//former plot_list
|
710 |
gStyle->SetOptStat("oueMri");
|
711 |
TCut observed=(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)");
|
712 |
|
713 |
TCut predicted = (cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
|
714 |
|| (sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
|
715 |
|| (sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)");
|
716 |
TCut predictedMC = (cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)")
|
717 |
|| (sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)")
|
718 |
|| (sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)");
|
719 |
|
720 |
flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
|
721 |
if(!PlottingSetup::RestrictToMassPeak) {
|
722 |
predicted = TCut((cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
|
723 |
predictedMC = TCut((cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)"));
|
724 |
}
|
725 |
// TCut predicted=((cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
|
726 |
// ||(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
|
727 |
// ||(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
|
728 |
// ||(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
|
729 |
// ||(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
|
730 |
// ||(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
|
731 |
// ||(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
|
732 |
|
733 |
bool doPF=false;
|
734 |
bool dolog=true;
|
735 |
bool nolog=false;
|
736 |
|
737 |
// Mll: do not scale
|
738 |
make_double_plot("mll",20,50,150,11,nolog,"m_{ll} [GeV]","mll",observed,predicted,data,true);
|
739 |
make_double_plot("mll",20,50,150,11,nolog,"m_{ll} [GeV]","mll_MC",observed,predictedMC,mc,true);
|
740 |
make_double_plot("met[4]",20,0,400,11,nolog,"pfMET [GeV]","pfMET",observed,predicted,data);
|
741 |
make_double_plot("met[4]",20,0,400,11,nolog,"pfMET [GeV]","pfMET_MC",observed,predictedMC,mc);
|
742 |
make_double_plot("jetpt[0]",10,0,400,11,nolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0",observed,predicted,data);
|
743 |
make_double_plot("jetpt[0]",10,0,400,11,nolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0_MC",observed,predictedMC,mc);
|
744 |
make_double_plot("jeteta[0]",10,-5,5,11,nolog,"leading jet #eta","pfJetGoodEta_0",observed,predicted,data);
|
745 |
make_double_plot("jeteta[0]",10,-5,5,11,nolog,"leading jet #eta","pfJetGoodEta_0_MC",observed,predictedMC,mc);
|
746 |
make_double_plot("pt",10,0,300,11,nolog,"Z p_{T} [GeV]","Zpt",observed,predicted,data);
|
747 |
make_double_plot("pt",10,0,300,11,nolog,"Z p_{T} [GeV]","Zpt_MC",observed,predictedMC,mc);
|
748 |
make_double_plot("pt1",10,0,200,15,nolog,"p_{T} [GeV]","pt1",observed,predicted,data);
|
749 |
make_double_plot("pt1",10,0,200,15,nolog,"p_{T} [GeV]","pt1_MC",observed,predictedMC,mc);
|
750 |
make_double_plot("pt2",10,0,200,25,nolog,"p_{T} [GeV]","pt2",observed,predicted,data);
|
751 |
make_double_plot("pt2",10,0,200,25,nolog,"p_{T} [GeV]","pt2_MC",observed,predictedMC,mc);
|
752 |
make_double_plot("eta1",10,-5,5,11,nolog,"#eta_{1,l}","eta_1",observed,predicted,data);
|
753 |
make_double_plot("eta1",10,-5,5,11,nolog,"#eta_{1,l}","eta_1_MC",observed,predictedMC,mc);
|
754 |
make_double_plot("eta2",10,-5,5,11,nolog,"#eta_{2,l}","eta_2",observed,predicted,data);
|
755 |
make_double_plot("eta2",10,-5,5,11,nolog,"#eta_{2,l}","eta_2_MC",observed,predictedMC,mc);
|
756 |
make_double_plot("phi1-phi2",10,-6.0,6.0,11,nolog,"#phi_{1}-#phi_{2}","dphi",observed,predicted,data);
|
757 |
make_double_plot("phi1-phi2",10,-6.0,6.0,11,nolog,"#phi_{1}-#phi_{2}","dphi_MC",observed,predictedMC,mc);
|
758 |
make_double_plot("pfJetGoodNum",8,0.5,8.5,20,nolog,"nJets","nJets",observed,predicted,data);
|
759 |
make_double_plot("pfJetGoodNum",8,0.5,8.5,20,nolog,"nJets","nJets_MC",observed,predictedMC,mc);
|
760 |
|
761 |
}
|
762 |
|
763 |
|
764 |
void generator_study_plots() {
|
765 |
|
766 |
//uncomment whichever one you want to see :-)
|
767 |
|
768 |
/*
|
769 |
GenLevelStudy::X_vs_generation_lm4();
|
770 |
GenLevelStudy::compare_sms_lm4();
|
771 |
*/
|
772 |
|
773 |
///GenLevelStudy::MomentumFraction();
|
774 |
///GenLevelStudy::AngleLSPLSP();
|
775 |
///GenLevelStudy::DeltaLSPmomentum();
|
776 |
///GenLevelStudy::ZDecayIllustration();
|
777 |
/*
|
778 |
GenLevelStudy::AngleMETsumLSP();
|
779 |
GenLevelStudy::RatioMETsumLSP();
|
780 |
GenLevelStudy::AngleLSPZ();
|
781 |
GenLevelStudy::WidthIllustration();
|
782 |
|
783 |
|
784 |
GenLevelStudy::pStarIllustration(0.5);
|
785 |
//GenLevelStudy::pStarIllustration(0.25);
|
786 |
//GenLevelStudy::pStarIllustration(0.75);
|
787 |
GenLevelStudy::DrawJetBand();
|
788 |
GenLevelStudy::DrawOnly100to150inJetBand();*/
|
789 |
GenLevelStudy::ImpactOfGluinoChi2MassDifference();
|
790 |
|
791 |
}
|
792 |
|
793 |
|
794 |
void jzb_negative_generator_study() {
|
795 |
write_warning(__FUNCTION__,"We are going to use a t5zz scan file, and \033[1;31m WON'T \033[1;35m cut on MassGlu/MassLSP in order to improve statistics. This is ok for small studies, for a real study you'll need to look at points individually ...");
|
796 |
|
797 |
|
798 |
scansample.AddSample("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v7/SMS-T5zz_x-05_Mgluino-150to1200_mLSP-50to1150_7TeV-Pythia6Z__Summer11-PU_START42_V11_FastSim-v2__AODSIM___inclindex_v2.root","SMST5zz",1,1,false,false,0,kRed);
|
799 |
TCanvas *jcan = new TCanvas("jcan","jcan");
|
800 |
scansample.collection[scansample.collection.size()-1].events->Draw("(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
|
801 |
TH1F *h = new TH1F("h","h",100,-500,500);
|
802 |
h->SetLineColor(kBlack);
|
803 |
TProfile *p = (TProfile*)jcan->GetPrimitive("htemp");
|
804 |
p->GetXaxis()->SetTitle("Generator JZB");
|
805 |
p->GetXaxis()->CenterTitle();
|
806 |
p->GetYaxis()->SetTitle("( LSP p_{T} ) / ( LSP mother p_{T} )");
|
807 |
p->GetYaxis()->CenterTitle();
|
808 |
p->SetLineColor(kBlue);
|
809 |
TLegend* leg = make_legend("", 0.40, 0.75, false);
|
810 |
leg->AddEntry(p,"LSP1 pt / LSP1 Mother pt","l");
|
811 |
leg->AddEntry(h,"Z pt / Z Mother pt","l");
|
812 |
leg->Draw();
|
813 |
TText *title = write_title("JZB as a function of the first LSP's momentum transfer");
|
814 |
title->Draw();
|
815 |
scansample.collection[scansample.collection.size()-1].events->Draw("(genZPt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF,same");
|
816 |
|
817 |
CompleteSave(jcan,"NegativeJZBStudy/LSPpt_LSPMopt");
|
818 |
|
819 |
scansample.collection[scansample.collection.size()-1].events->Draw("angleLSPLSP:pureGeneratorJZB","genNjets>2","PROF");
|
820 |
TProfile *p1 = (TProfile*)jcan->GetPrimitive("htemp");
|
821 |
p1->GetXaxis()->SetTitle("Generator JZB");
|
822 |
p1->GetXaxis()->CenterTitle();
|
823 |
p1->GetYaxis()->SetTitle("#angle(LSP1,LSP2)");
|
824 |
p1->GetYaxis()->CenterTitle();
|
825 |
TText *title1 = write_title("JZB as a function of the angle between the two LSPs");
|
826 |
title1->Draw();
|
827 |
CompleteSave(jcan,"NegativeJZBStudy/AngleLSPLSP");
|
828 |
|
829 |
TH1F *jzbdistributionvsz[5];
|
830 |
THStack zstack;
|
831 |
jcan->SetLogy(1);
|
832 |
TLegend* leg2 = make_legend("", 0.15, 0.75, false);
|
833 |
leg2->SetX2(0.4);
|
834 |
TLegend* leg3 = make_legend("", 0.15, 0.75, false);
|
835 |
leg3->SetX2(0.4);
|
836 |
for(int z=0;z<5;z++) {
|
837 |
stringstream specialcut;
|
838 |
if(z<4) specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2 << "&&(LSP1pt/LSP1Mopt)<" << (z+1)*0.2;
|
839 |
else specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2;
|
840 |
stringstream histtitle;
|
841 |
histtitle << "splitup_" << z;
|
842 |
stringstream ntitle;
|
843 |
if(z<4) ntitle << z*0.2 << " < z < " << (z+1)*0.2;
|
844 |
else ntitle << z*0.2 << " < z";
|
845 |
jzbdistributionvsz[z] = scansample.Draw(histtitle.str().c_str(),"pureGeneratorJZB",100,-500,500, "generator JZB (GeV)", "events",specialcut.str().c_str(),mc,1.0);
|
846 |
jzbdistributionvsz[z]->SetLineColor(z+1);
|
847 |
jzbdistributionvsz[z]->SetMarkerSize(0);
|
848 |
zstack.Add(jzbdistributionvsz[z]);
|
849 |
leg2->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"f");
|
850 |
leg3->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"l");
|
851 |
}
|
852 |
|
853 |
jzbdistributionvsz[0]->GetYaxis()->SetRangeUser(2,800);
|
854 |
jzbdistributionvsz[0]->DrawNormalized("histo");
|
855 |
for(int z=0;z<5;z++) {
|
856 |
jzbdistributionvsz[z]->DrawNormalized("histo,same");
|
857 |
}
|
858 |
|
859 |
// zstack.Draw("nostack,histo");
|
860 |
leg3->Draw("same");
|
861 |
CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionIndividual");
|
862 |
|
863 |
for(int z=0;z<5;z++) {
|
864 |
jzbdistributionvsz[z]->SetFillColor(z+1);
|
865 |
}
|
866 |
|
867 |
zstack.Draw("histo");
|
868 |
leg2->Draw("same");
|
869 |
CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionStacked");
|
870 |
|
871 |
//varangle vasysyn
|
872 |
|
873 |
// scansample.collection[scansample.collection.size()-1].events->Draw("(LSPPromptnessLevel[0]/4.0)*angleLSPLSP/(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
|
874 |
|
875 |
|
876 |
|
877 |
}
|
878 |
|
879 |
void compare_lm4_sms_variable(TTree *eventsLM4, TTree *eventsSMS, string variable, int nbins, float xmin, float xmax, TCut cut, string saveas, bool dology=false) {
|
880 |
TCanvas *can = new TCanvas("can","can");
|
881 |
can->SetLogy(dology);
|
882 |
TH1F *hlm4 = new TH1F("hlm4","hlm4",nbins,xmin,xmax);
|
883 |
TH1F *hsms = new TH1F("hsms","hsms",nbins,xmin,xmax);
|
884 |
eventsLM4->Draw((variable+">>hlm4").c_str(),cut,"goff");
|
885 |
eventsSMS->Draw((variable+">>hsms").c_str(),cut,"goff");
|
886 |
hlm4->SetLineColor(kBlue);
|
887 |
hsms->SetLineColor(kRed);
|
888 |
|
889 |
if(hlm4->Integral()>0) hlm4->Scale(1.0/hlm4->Integral());
|
890 |
else write_warning(__FUNCTION__,"Watch out, LM4 histo is empty!");
|
891 |
if(hsms->Integral()>0) hsms->Scale(1.0/hsms->Integral());
|
892 |
else write_warning(__FUNCTION__,"Watch out, SMS histo is empty!");
|
893 |
|
894 |
float min=get_nonzero_minimum(hlm4);
|
895 |
float max=hlm4->GetMaximum();
|
896 |
if(get_nonzero_minimum(hsms)<min) min=get_nonzero_minimum(hsms);
|
897 |
if(hsms->GetMaximum()>max) max=hsms->GetMaximum();
|
898 |
if(dology) max*=5;
|
899 |
else max*=2;
|
900 |
|
901 |
hlm4->GetYaxis()->SetRangeUser(min,max);
|
902 |
hlm4->GetXaxis()->SetTitle(variable.c_str());
|
903 |
hlm4->GetXaxis()->CenterTitle();
|
904 |
hlm4->Draw("histo");
|
905 |
hsms->Draw("histo,same");
|
906 |
TLegend *leg = make_legend("",0.2,0.98,false);
|
907 |
leg->SetY2(1.0);
|
908 |
leg->SetNColumns(2);
|
909 |
leg->AddEntry(hlm4,"LM4","l");
|
910 |
leg->AddEntry(hsms,"\"LM4\" SMS","l");
|
911 |
leg->Draw();
|
912 |
stringstream saveas2;
|
913 |
saveas2 << "ComparingLM4_SMS/" << removefunnystring(saveas);
|
914 |
CompleteSave(can,saveas2.str());
|
915 |
delete can;
|
916 |
delete hlm4;
|
917 |
delete hsms;
|
918 |
}
|
919 |
|
920 |
|
921 |
void compare_LM4_and_SMS() {
|
922 |
TFile *f1 = new TFile("/shome/lbaeni/jzb/LM4_SMS/SMS_LM4_JZB.root");
|
923 |
TTree *LM4events = (TTree*)f1->Get("events");
|
924 |
TFile *f2 = new TFile("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v5/LM4_SUSY_sftsht_7TeV-pythia6__Summer11-PU_S4_START42_V11-v2__withIndex.root");
|
925 |
TTree *SMSevents = (TTree*)f2->Get("events");
|
926 |
|
927 |
compare_lm4_sms_variable(LM4events, SMSevents, "mll",100,50,150,cutOSSF&&cutnJets&&basiccut,"mll",true);
|
928 |
compare_lm4_sms_variable(LM4events, SMSevents, "jzb[1]+0.04*pt-2.32212",100,-300,700,cutmass&&cutOSSF&&cutnJets&&basiccut,"jzb",true);
|
929 |
compare_lm4_sms_variable(LM4events, SMSevents, "pureGeneratorJZB",100,-300.0,700.0,cutOSSF&&basiccut,"pureGeneratorJZB",true);
|
930 |
compare_lm4_sms_variable(LM4events, SMSevents, "pfJetGoodNum",10,-0.5,9.5,cutmass&&cutOSSF&&basiccut,"pfJetGoodNum",true);
|
931 |
compare_lm4_sms_variable(LM4events, SMSevents, "pt",100,15.0,200.0,cutOSSF&&basiccut,"pt",false);
|
932 |
compare_lm4_sms_variable(LM4events, SMSevents, "pt1",100,15.0,100.0,cutOSSF&&basiccut,"pt1",false);
|
933 |
compare_lm4_sms_variable(LM4events, SMSevents, "pt2",100,15.0,100.0,cutOSSF&&basiccut,"pt2",false);
|
934 |
compare_lm4_sms_variable(LM4events, SMSevents, "met[4]",100,0.0,600.0,cutOSSF&&basiccut,"met",false);
|
935 |
compare_lm4_sms_variable(LM4events, SMSevents, "genMET",100,0.0,600.0,cutOSSF&&basiccut,"genMET",false);
|
936 |
compare_lm4_sms_variable(LM4events, SMSevents, "genNjets",10,-0.5,9.5,cutOSSF&&basiccut,"genNjets",true);
|
937 |
}
|
938 |
|
939 |
|
940 |
void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb, int sampleindex) {
|
941 |
cout << (signalsamples.collection)[sampleindex].filename << endl;
|
942 |
}
|
943 |
|
944 |
void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb) {
|
945 |
for(int i=0;i<signalsamples.collection.size();i++) compare_onpeak_offpeak_signal_distributions(mcjzb,datajzb,i);
|
946 |
}
|
947 |
|
948 |
bool load_adequate_npred_nobs(string code){
|
949 |
|
950 |
if(code=="201540") {
|
951 |
Nobs[0]=1108.67;
|
952 |
Npred[0]=1106.69;
|
953 |
Nprederr[0]=275.424852964089;
|
954 |
Nobs[1]=291.679;
|
955 |
Npred[1]=288.979;
|
956 |
Nprederr[1]=74.0248355630055;
|
957 |
Nobs[2]=76.0776;
|
958 |
Npred[2]=76.8332;
|
959 |
Nprederr[2]=20.9768252604726;
|
960 |
Nobs[3]=22.4581;
|
961 |
Npred[3]=22.5385;
|
962 |
Nprederr[3]=7.46758781734772;
|
963 |
Nobs[4]=5.41088;
|
964 |
Npred[4]=9.51903;
|
965 |
Nprederr[4]=3.86595196996807;
|
966 |
return true;
|
967 |
}
|
968 |
|
969 |
if(code=="201040") {
|
970 |
Nobs[0]=1210.02;
|
971 |
Npred[0]=1228.62;
|
972 |
Nprederr[0]=305.642034977193;
|
973 |
Nobs[1]=324.008;
|
974 |
Npred[1]=323.042;
|
975 |
Nprederr[1]=82.6094009556418;
|
976 |
Nobs[2]=84.8866;
|
977 |
Npred[2]=86.4243;
|
978 |
Nprederr[2]=23.4841694480367;
|
979 |
Nobs[3]=25.4339;
|
980 |
Npred[3]=25.6437;
|
981 |
Nprederr[3]=8.32541119033168;
|
982 |
Nobs[4]=6.24946;
|
983 |
Npred[4]=9.97731;
|
984 |
Nprederr[4]=4.07510841593202;
|
985 |
return true;
|
986 |
}
|
987 |
|
988 |
if(code=="151040") {
|
989 |
Nobs[0]=1218.67;
|
990 |
Npred[0]=1239.47;
|
991 |
Nprederr[0]=308.268321394268;
|
992 |
Nobs[1]=326.245;
|
993 |
Npred[1]=325.731;
|
994 |
Nprederr[1]=83.2795749634327;
|
995 |
Nobs[2]=85.3292;
|
996 |
Npred[2]=87.239;
|
997 |
Nprederr[2]=23.6875859996856;
|
998 |
Nobs[3]=25.4339;
|
999 |
Npred[3]=25.9456;
|
1000 |
Nprederr[3]=8.40219371514963;
|
1001 |
Nobs[4]=6.24946;
|
1002 |
Npred[4]=9.97731;
|
1003 |
Nprederr[4]=4.07510841593202;
|
1004 |
return true;
|
1005 |
}
|
1006 |
|
1007 |
if(code=="151550") {
|
1008 |
Nobs[0]=1009.57;
|
1009 |
Npred[0]=1015.24;
|
1010 |
Nprederr[0]=252.027368634063;
|
1011 |
Nobs[1]=266.851;
|
1012 |
Npred[1]=264.99;
|
1013 |
Nprederr[1]=67.9766170833765;
|
1014 |
Nobs[2]=68.2042;
|
1015 |
Npred[2]=69.3185;
|
1016 |
Nprederr[2]=19.0998669414868;
|
1017 |
Nobs[3]=20.873;
|
1018 |
Npred[3]=20.1324;
|
1019 |
Nprederr[3]=6.85046567121535;
|
1020 |
Nobs[4]=4.67275;
|
1021 |
Npred[4]=8.42002;
|
1022 |
Nprederr[4]=3.55560434037591;
|
1023 |
return true;
|
1024 |
}
|
1025 |
|
1026 |
if(code=="151540") {
|
1027 |
Nobs[0]=1112.51;
|
1028 |
Npred[0]=1111.01;
|
1029 |
Nprederr[0]=276.447968514945;
|
1030 |
Nobs[1]=292.622;
|
1031 |
Npred[1]=290.731;
|
1032 |
Nprederr[1]=74.461488188526;
|
1033 |
Nobs[2]=76.1949;
|
1034 |
Npred[2]=77.1915;
|
1035 |
Nprederr[2]=21.0662269219811;
|
1036 |
Nobs[3]=22.4581;
|
1037 |
Npred[3]=22.8404;
|
1038 |
Nprederr[3]=7.54484166514447;
|
1039 |
Nobs[4]=5.41088;
|
1040 |
Npred[4]=9.51903;
|
1041 |
Nprederr[4]=3.86595196996807;
|
1042 |
return true;
|
1043 |
}
|
1044 |
|
1045 |
if(code=="151530") {
|
1046 |
Nobs[0]=1182.29;
|
1047 |
Npred[0]=1177.48;
|
1048 |
Nprederr[0]=293.397396547567;
|
1049 |
Nobs[1]=310.125;
|
1050 |
Npred[1]=308.279;
|
1051 |
Nprederr[1]=78.8363941324691;
|
1052 |
Nobs[2]=83.376;
|
1053 |
Npred[2]=80.804;
|
1054 |
Nprederr[2]=21.9686562734479;
|
1055 |
Nobs[3]=24.9171;
|
1056 |
Npred[3]=23.4465;
|
1057 |
Nprederr[3]=7.69983840500565;
|
1058 |
Nobs[4]=6.151;
|
1059 |
Npred[4]=9.7932;
|
1060 |
Nprederr[4]=3.94255080734542;
|
1061 |
return true;
|
1062 |
}
|
1063 |
|
1064 |
if(code=="101040") {
|
1065 |
Nobs[0]=1219.78;
|
1066 |
Npred[0]=1241.65;
|
1067 |
Nprederr[0]=308.795385452892;
|
1068 |
Nobs[1]=326.577;
|
1069 |
Npred[1]=327.172;
|
1070 |
Nprederr[1]=83.6388357552878;
|
1071 |
Nobs[2]=85.3292;
|
1072 |
Npred[2]=87.5409;
|
1073 |
Nprederr[2]=23.7629564995099;
|
1074 |
Nobs[3]=25.4339;
|
1075 |
Npred[3]=26.2475;
|
1076 |
Nprederr[3]=8.47896339395919;
|
1077 |
Nobs[4]=6.24946;
|
1078 |
Npred[4]=9.97731;
|
1079 |
Nprederr[4]=4.07510841593202;
|
1080 |
return true;
|
1081 |
}
|
1082 |
|
1083 |
if(code=="101030") {
|
1084 |
Nobs[0]=1335.36;
|
1085 |
Npred[0]=1342.37;
|
1086 |
Nprederr[0]=334.276222135228;
|
1087 |
Nobs[1]=355.554;
|
1088 |
Npred[1]=355.065;
|
1089 |
Nprederr[1]=90.5352675557984;
|
1090 |
Nobs[2]=94.9863;
|
1091 |
Npred[2]=93.7287;
|
1092 |
Nprederr[2]=25.3082516747483;
|
1093 |
Nobs[3]=28.3785;
|
1094 |
Npred[3]=28.3783;
|
1095 |
Nprederr[3]=9.02022017545026;
|
1096 |
Nobs[4]=6.98957;
|
1097 |
Npred[4]=10.8503;
|
1098 |
Nprederr[4]=4.31460183234792;
|
1099 |
return true;
|
1100 |
}
|
1101 |
|
1102 |
if(code=="101020") {
|
1103 |
Nobs[0]=1418.83;
|
1104 |
Npred[0]=1421.23;
|
1105 |
Nprederr[0]=353.214726699284;
|
1106 |
Nobs[1]=387.594;
|
1107 |
Npred[1]=380.483;
|
1108 |
Nprederr[1]=96.9195932538411;
|
1109 |
Nobs[2]=106.565;
|
1110 |
Npred[2]=99.4694;
|
1111 |
Nprederr[2]=26.7421279631969;
|
1112 |
Nobs[3]=32.9309;
|
1113 |
Npred[3]=31.6286;
|
1114 |
Nprederr[3]=9.84399940108186;
|
1115 |
Nobs[4]=9.25395;
|
1116 |
Npred[4]=12.5158;
|
1117 |
Nprederr[4]=4.76585512033255;
|
1118 |
return true;
|
1119 |
}
|
1120 |
|
1121 |
if(code=="202040") {
|
1122 |
Nobs[0]=960.506;
|
1123 |
Npred[0]=954.77;
|
1124 |
Nprederr[0]=238.5811858172;
|
1125 |
Nobs[1]=251.453;
|
1126 |
Npred[1]=250.838;
|
1127 |
Nprederr[1]=64.2967104354;
|
1128 |
Nobs[2]=64.4903;
|
1129 |
Npred[2]=66.6518;
|
1130 |
Nprederr[2]=18.4320150537;
|
1131 |
Nobs[3]=19.1953;
|
1132 |
Npred[3]=20.1075;
|
1133 |
Nprederr[3]=6.844063617;
|
1134 |
Nobs[4]=5.09139;
|
1135 |
Npred[4]=8.67175;
|
1136 |
Nprederr[4]=3.6271903178;
|
1137 |
return true;
|
1138 |
}
|
1139 |
|
1140 |
if(code=="onpeak") {
|
1141 |
Nobs[0]=387.268;
|
1142 |
Npred[0]=377.146;
|
1143 |
Nprederr[0]=55.4000603791;
|
1144 |
Nobs[1]=93.6071;
|
1145 |
Npred[1]=89.2911;
|
1146 |
Nprederr[1]=13.7946793213;
|
1147 |
Nobs[2]=22.4222;
|
1148 |
Npred[2]=22.2462;
|
1149 |
Nprederr[2]=4.252924039;
|
1150 |
Nobs[3]=7.49126;
|
1151 |
Npred[3]=6.98267;
|
1152 |
Nprederr[3]=1.913797227;
|
1153 |
Nobs[4]=1.80426;
|
1154 |
Npred[4]=3.00437;
|
1155 |
Nprederr[4]=1.2381643019;
|
1156 |
return true;
|
1157 |
}
|
1158 |
|
1159 |
if(code=="151550") {
|
1160 |
Nobs[0]=1009.57;
|
1161 |
Npred[0]=1015.24;
|
1162 |
Nprederr[0]=252.0273686341;
|
1163 |
Nobs[1]=266.851;
|
1164 |
Npred[1]=264.99;
|
1165 |
Nprederr[1]=67.9766170834;
|
1166 |
Nobs[2]=68.2042;
|
1167 |
Npred[2]=69.3185;
|
1168 |
Nprederr[2]=19.0998669415;
|
1169 |
Nobs[3]=20.873;
|
1170 |
Npred[3]=20.1324;
|
1171 |
Nprederr[3]=6.8504656712;
|
1172 |
Nobs[4]=4.67275;
|
1173 |
Npred[4]=8.42002;
|
1174 |
Nprederr[4]=3.5556043404;
|
1175 |
return true;
|
1176 |
}
|
1177 |
|
1178 |
if(code=="101040") {
|
1179 |
Nobs[0]=1108.67;
|
1180 |
Npred[0]=1106.69;
|
1181 |
Nprederr[0]=275.4248529641;
|
1182 |
Nobs[1]=291.679;
|
1183 |
Npred[1]=288.979;
|
1184 |
Nprederr[1]=74.024835563;
|
1185 |
Nobs[2]=76.0776;
|
1186 |
Npred[2]=76.8332;
|
1187 |
Nprederr[2]=20.9768252605;
|
1188 |
Nobs[3]=22.4581;
|
1189 |
Npred[3]=22.5385;
|
1190 |
Nprederr[3]=7.4675878173;
|
1191 |
Nobs[4]=5.41088;
|
1192 |
Npred[4]=9.51903;
|
1193 |
Nprederr[4]=3.86595197;
|
1194 |
return true;
|
1195 |
}
|
1196 |
|
1197 |
if(code=="151040") {
|
1198 |
Nobs[0]=1218.67;
|
1199 |
Npred[0]=1239.47;
|
1200 |
Nprederr[0]=308.2683213943;
|
1201 |
Nobs[1]=326.245;
|
1202 |
Npred[1]=325.731;
|
1203 |
Nprederr[1]=83.2795749634;
|
1204 |
Nobs[2]=85.3292;
|
1205 |
Npred[2]=87.239;
|
1206 |
Nprederr[2]=23.6875859997;
|
1207 |
Nobs[3]=25.4339;
|
1208 |
Npred[3]=25.9456;
|
1209 |
Nprederr[3]=8.4021937151;
|
1210 |
Nobs[4]=6.24946;
|
1211 |
Npred[4]=9.97731;
|
1212 |
Nprederr[4]=4.0751084159;
|
1213 |
return true;
|
1214 |
}
|
1215 |
|
1216 |
if(code=="201040") {
|
1217 |
Nobs[0]=1210.02;
|
1218 |
Npred[0]=1228.62;
|
1219 |
Nprederr[0]=305.6420349772;
|
1220 |
Nobs[1]=324.008;
|
1221 |
Npred[1]=323.042;
|
1222 |
Nprederr[1]=82.6094009556;
|
1223 |
Nobs[2]=84.8866;
|
1224 |
Npred[2]=86.4243;
|
1225 |
Nprederr[2]=23.484169448;
|
1226 |
Nobs[3]=25.4339;
|
1227 |
Npred[3]=25.6437;
|
1228 |
Nprederr[3]=8.3254111903;
|
1229 |
Nobs[4]=6.24946;
|
1230 |
Npred[4]=9.97731;
|
1231 |
Nprederr[4]=4.0751084159;
|
1232 |
return true;
|
1233 |
}
|
1234 |
|
1235 |
if(code=="151540") {
|
1236 |
Nobs[0]=1112.51;
|
1237 |
Npred[0]=1111.01;
|
1238 |
Nprederr[0]=276.4479685149;
|
1239 |
Nobs[1]=292.622;
|
1240 |
Npred[1]=290.731;
|
1241 |
Nprederr[1]=74.4614881885;
|
1242 |
Nobs[2]=76.1949;
|
1243 |
Npred[2]=77.1915;
|
1244 |
Nprederr[2]=21.066226922;
|
1245 |
Nobs[3]=22.4581;
|
1246 |
Npred[3]=22.8404;
|
1247 |
Nprederr[3]=7.5448416651;
|
1248 |
Nobs[4]=5.41088;
|
1249 |
Npred[4]=9.51903;
|
1250 |
Nprederr[4]=3.86595197;
|
1251 |
return true;
|
1252 |
}
|
1253 |
|
1254 |
|
1255 |
return false;
|
1256 |
}
|
1257 |
|
1258 |
vector<float> do_simulate_upper_limits(string mcjzb, string datajzb, float MCPeakError, vector<float> jzb_cut, string code) {
|
1259 |
vector<vector<float> > all_systematics;
|
1260 |
vector<float> bestUL;
|
1261 |
//we need to set the correct npred, nprederr, and nobs.
|
1262 |
bool loaded_info = load_adequate_npred_nobs(code);
|
1263 |
if(!loaded_info) {
|
1264 |
write_warning(__FUNCTION__,"obs/pred/prederr could not be loaded for this configuration.");
|
1265 |
cout << "Configuration " << code << " (pt1/pt2/mll) caused problems. will be skipped." << endl;
|
1266 |
} else {
|
1267 |
cout << "Loaded configuration for code " << code << " successfully, e.g. Npred[0] is " << Npred[0] << endl;
|
1268 |
bool doobserved=false;
|
1269 |
all_systematics=compute_systematics(mcjzb,MCPeakError,alwaysflip,datajzb,signalsamples,jzb_cut,requireZ);
|
1270 |
bestUL = compute_upper_limits_from_counting_experiment(all_systematics,jzb_cut,mcjzb,doobserved,alwaysflip);
|
1271 |
|
1272 |
}
|
1273 |
return bestUL;
|
1274 |
}
|
1275 |
|
1276 |
void decipherM0M12sample(string samplename, float &M0, float &M12) {
|
1277 |
int position = samplename.find("M0_");
|
1278 |
string interestingpart = samplename.substr(position+3,4);
|
1279 |
position = interestingpart.find("_");
|
1280 |
if(position>0&&position<5) interestingpart=interestingpart.substr(0,position);
|
1281 |
stringstream M0a;
|
1282 |
M0a << interestingpart;
|
1283 |
M0a>>M0;
|
1284 |
position = samplename.find("M12_");
|
1285 |
interestingpart = samplename.substr(position+4,4);
|
1286 |
position = interestingpart.find("_");
|
1287 |
if(position>0&&position<5) interestingpart=interestingpart.substr(0,position);
|
1288 |
stringstream M12a;
|
1289 |
M12a << interestingpart;
|
1290 |
M12a>>M12;
|
1291 |
}
|
1292 |
|
1293 |
vector<TGraph*> create_limit_graphs(vector<float> limits) {
|
1294 |
float x[limits.size()];
|
1295 |
float y[limits.size()];
|
1296 |
int nx=0;
|
1297 |
float x2[limits.size()];
|
1298 |
float y2[limits.size()];
|
1299 |
int nx2=0;
|
1300 |
for(int i=0;i<limits.size();i++) {
|
1301 |
string samplename=signalsamples.collection[i].samplename;
|
1302 |
float m0,m12;
|
1303 |
decipherM0M12sample(samplename,m0,m12);
|
1304 |
if(m12==300) {
|
1305 |
x[nx]=m0;
|
1306 |
y[nx]=limits[i];
|
1307 |
nx++;
|
1308 |
}
|
1309 |
if(m12==400) {
|
1310 |
x2[nx2]=m0;
|
1311 |
y2[nx2]=limits[i];
|
1312 |
nx2++;
|
1313 |
}
|
1314 |
}
|
1315 |
|
1316 |
TGraph *gra = new TGraph(nx,x,y);
|
1317 |
TGraph *grb = new TGraph(nx2,x2,y2);
|
1318 |
vector<TGraph*> graphs;
|
1319 |
graphs.push_back(gra);
|
1320 |
graphs.push_back(grb);
|
1321 |
return graphs;
|
1322 |
}
|
1323 |
|
1324 |
|
1325 |
|
1326 |
void interpret_bestULs_for_cuts(vector<pair<vector<float>,string> > bestUL) {
|
1327 |
TCanvas *can = new TCanvas("can","can",1800,900);
|
1328 |
can->Divide(2,1);
|
1329 |
can->cd(1);
|
1330 |
TGraph *graphs[2][bestUL.size()];
|
1331 |
TLegend *leg = make_legend();
|
1332 |
for(int i=0;i<bestUL.size();i++) {
|
1333 |
vector<TGraph*> grs=create_limit_graphs(bestUL[i].first);
|
1334 |
graphs[0][i]=grs[0];
|
1335 |
graphs[1][i]=grs[1];
|
1336 |
graphs[0][i]->SetLineColor(GenLevelStudy::diversehistocolor(i));
|
1337 |
graphs[1][i]->SetLineColor(GenLevelStudy::diversehistocolor(i));
|
1338 |
graphs[0][i]->SetMarkerColor(GenLevelStudy::diversehistocolor(i));
|
1339 |
graphs[1][i]->SetMarkerColor(GenLevelStudy::diversehistocolor(i));
|
1340 |
leg->AddEntry(graphs[0][i],bestUL[i].second.c_str(),"lp");
|
1341 |
}
|
1342 |
|
1343 |
|
1344 |
for(int i=0;i<bestUL.size();i++) {
|
1345 |
graphs[0][i]->GetXaxis()->SetTitle("M_{0}");
|
1346 |
graphs[0][i]->GetYaxis()->SetTitle("Upper Limit (M_{1/2}=300)");
|
1347 |
graphs[0][i]->GetXaxis()->CenterTitle();
|
1348 |
graphs[0][i]->GetYaxis()->CenterTitle();
|
1349 |
|
1350 |
graphs[1][i]->GetXaxis()->SetTitle("M_{0}");
|
1351 |
graphs[1][i]->GetYaxis()->SetTitle("Upper Limit (M_{1/2}=400)");
|
1352 |
graphs[1][i]->GetXaxis()->CenterTitle();
|
1353 |
graphs[1][i]->GetYaxis()->CenterTitle();
|
1354 |
|
1355 |
can->cd(1);
|
1356 |
if(i==0) graphs[0][i]->Draw("ALP");
|
1357 |
else graphs[0][i]->Draw("LP");
|
1358 |
can->cd(2);
|
1359 |
if(i==0) graphs[1][i]->Draw("ALP");
|
1360 |
else graphs[1][i]->Draw("LP");
|
1361 |
}
|
1362 |
can->cd(1);
|
1363 |
leg->Draw();
|
1364 |
can->cd(2);
|
1365 |
leg->Draw();
|
1366 |
CompleteSave(can,"mSUGRAlimits");
|
1367 |
delete can;
|
1368 |
}
|
1369 |
|
1370 |
void do_simulate_upper_limits(string mcjzb, string datajzb, float MCPeakError, vector<float> jzb_cut) {
|
1371 |
|
1372 |
cout << "Going to simulate upper limits (i.e. expected limits) USING MONTE CARLO!!!!" << endl;
|
1373 |
while(jzb_cut.size()>Nobs.size()||jzb_cut.size()>Npred.size()||jzb_cut.size()>Nprederr.size()) {
|
1374 |
Nobs.push_back(0);
|
1375 |
Npred.push_back(0);
|
1376 |
Nprederr.push_back(0);
|
1377 |
}
|
1378 |
|
1379 |
vector<pair<vector<float>,string> > bestUL;
|
1380 |
|
1381 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201540"),"201540"));
|
1382 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201040"),"201040"));/*
|
1383 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151040"),"151040"));
|
1384 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151550"),"151550"));
|
1385 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151540"),"151540"));
|
1386 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151530"),"151530"));
|
1387 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101040"),"101040"));
|
1388 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101030"),"101030"));
|
1389 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101020"),"101020"));
|
1390 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151540"),"151540"));
|
1391 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151550"),"151550"));
|
1392 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"201040"),"201040"));
|
1393 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"151040"),"151040"));
|
1394 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"101040"),"101040"));
|
1395 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"onpeak"),"onpeak"));
|
1396 |
bestUL.push_back(make_pair(do_simulate_upper_limits(mcjzb,datajzb,MCPeakError,jzb_cut,"202040"),"202040"));
|
1397 |
*/
|
1398 |
/*
|
1399 |
//-------------------
|
1400 |
write_warning(__FUNCTION__,"Just testing the components now ...");
|
1401 |
vector<float> limits;
|
1402 |
for(int i=0;i<signalsamples.collection.size();i++) limits.push_back(1/(0.01*(i+1)));
|
1403 |
vector<float> limits2;
|
1404 |
for(int i=0;i<signalsamples.collection.size();i++) limits2.push_back(1/(0.01*(2*i+2)));
|
1405 |
bestUL.push_back(make_pair(limits,"abcd"));
|
1406 |
bestUL.push_back(make_pair(limits2,"efgh"));
|
1407 |
|
1408 |
//-------------------
|
1409 |
*/
|
1410 |
|
1411 |
|
1412 |
interpret_bestULs_for_cuts(bestUL);
|
1413 |
}
|