| 22 |
|
|
| 23 |
|
|
| 24 |
|
void rediscover_the_top(string mcjzb, string datajzb) { |
| 25 |
< |
cout << "Hi! today we are going to (try to) rediscover the top!" << endl; |
| 25 |
> |
dout << "Hi! today we are going to (try to) rediscover the top!" << endl; |
| 26 |
|
TCanvas *c3 = new TCanvas("c3","c3"); |
| 27 |
|
c3->SetLogy(1); |
| 28 |
|
vector<float> binning; |
| 62 |
|
TH1F *datapredictiont = allsamples.Draw("datapredictiont", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity); |
| 63 |
|
TH1F *datapredictiono = allsamples.Draw("datapredictiono", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity); |
| 64 |
|
datapredictiont->Add(datapredictiono,-1); |
| 65 |
< |
cout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl; |
| 65 |
> |
dout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl; |
| 66 |
|
vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10); |
| 67 |
|
datapredictiont->SetMarkerColor(kRed); |
| 68 |
|
datapredictiont->SetLineColor(kRed); |
| 141 |
|
ratio->Divide(datapredictiontb); |
| 142 |
|
|
| 143 |
|
for (int i=0;i<=ratio_binning.size();i++) { |
| 144 |
< |
cout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl; |
| 144 |
> |
dout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl; |
| 145 |
|
} |
| 146 |
|
|
| 147 |
|
// ratio->Divide(datapredictiontb); |
| 202 |
|
*/ |
| 203 |
|
} |
| 204 |
|
|
| 205 |
< |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts) { |
| 206 |
< |
cout << "Doing counting experiment ... " << endl; |
| 205 |
> |
vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, bool doobserved=false) { |
| 206 |
> |
float sigma95=0.0,sigma95A=0.0; |
| 207 |
> |
dout << "Now calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << 1<< ") " << endl; |
| 208 |
> |
sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, 1); |
| 209 |
> |
if(doobserved) { |
| 210 |
> |
dout << "Now calling : CL95A(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << 1<< ") " << endl; |
| 211 |
> |
sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], 1); |
| 212 |
> |
} |
| 213 |
> |
vector<float> sigmas; |
| 214 |
> |
sigmas.push_back(sigma95); |
| 215 |
> |
sigmas.push_back(sigma95A); |
| 216 |
> |
return sigmas; |
| 217 |
> |
} |
| 218 |
> |
|
| 219 |
> |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) { |
| 220 |
> |
dout << "Doing counting experiment ... " << endl; |
| 221 |
> |
vector<vector<string> > limits; |
| 222 |
> |
vector<vector<float> > vlimits; |
| 223 |
> |
|
| 224 |
|
|
| 225 |
|
for(int isample=0;isample<signalsamples.collection.size();isample++) { |
| 226 |
< |
cout << "Considering sample " << signalsamples.collection[isample].samplename << endl; |
| 226 |
> |
vector<string> rows; |
| 227 |
> |
vector<float> vrows; |
| 228 |
> |
dout << "Considering sample " << signalsamples.collection[isample].samplename << endl; |
| 229 |
> |
rows.push_back(signalsamples.collection[isample].samplename); |
| 230 |
|
for(int ibin=0;ibin<jzbcuts.size();ibin++) { |
| 231 |
< |
cout << " Considering bin " << jzbcuts[ibin] << endl; |
| 232 |
< |
for(int isyst=0;isyst<uncertainties[isample*jzbcuts.size()+ibin].size();isyst++) { |
| 233 |
< |
if(isyst==0) { |
| 234 |
< |
if(uncertainties[isample*jzbcuts.size()+ibin][isyst]!=jzbcuts[ibin]) cout << "WATCH OUT THERE SEEMS TO BE A PROBLEM - THE TEST BIN DOESN'T CORRESPOND TO YOUR BIN!" << endl; |
| 235 |
< |
continue; |
| 236 |
< |
} |
| 237 |
< |
cout << " Considering syst " << isyst << " : " << uncertainties[isample*jzbcuts.size()+ibin][isyst] << endl; |
| 238 |
< |
}//end of syst loop |
| 231 |
> |
dout << "_________________________________________________________________________________" << endl; |
| 232 |
> |
float JZBcutat=uncertainties[isample*jzbcuts.size()+ibin][0]; |
| 233 |
> |
float mceff=uncertainties[isample*jzbcuts.size()+ibin][1]; |
| 234 |
> |
float staterr=uncertainties[isample*jzbcuts.size()+ibin][2]; |
| 235 |
> |
float systerr=uncertainties[isample*jzbcuts.size()+ibin][3]; |
| 236 |
> |
float toterr =uncertainties[isample*jzbcuts.size()+ibin][4]; |
| 237 |
> |
float observed,null,result; |
| 238 |
> |
fill_result_histos(observed, null,null,null,null,null,null,null,mcjzb,JZBcutat,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples); |
| 239 |
> |
observed-=result;//this is the actual excess we see! |
| 240 |
> |
float expected=observed/luminosity; |
| 241 |
> |
|
| 242 |
> |
dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl; |
| 243 |
> |
vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,doobserved); |
| 244 |
> |
|
| 245 |
> |
if(doobserved) { |
| 246 |
> |
rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")"); |
| 247 |
> |
vrows.push_back(sigmas[0]); |
| 248 |
> |
vrows.push_back(sigmas[1]); |
| 249 |
> |
vrows.push_back(expected); |
| 250 |
> |
} |
| 251 |
> |
else { |
| 252 |
> |
rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")"); |
| 253 |
> |
vrows.push_back(sigmas[0]); |
| 254 |
> |
vrows.push_back(expected); |
| 255 |
> |
} |
| 256 |
|
}//end of bin loop |
| 257 |
+ |
limits.push_back(rows); |
| 258 |
+ |
vlimits.push_back(vrows); |
| 259 |
|
}//end of sample loop |
| 260 |
+ |
dout << endl << endl << "PAS table 3: " << endl << endl; |
| 261 |
+ |
dout << "\t"; |
| 262 |
+ |
for (int irow=0;irow<jzbcuts.size();irow++) { |
| 263 |
+ |
dout << jzbcuts[irow] << "\t"; |
| 264 |
+ |
} |
| 265 |
+ |
dout << endl; |
| 266 |
+ |
for(int irow=0;irow<limits.size();irow++) { |
| 267 |
+ |
for(int ientry=0;ientry<limits[irow].size();ientry++) { |
| 268 |
+ |
dout << limits[irow][ientry] << "\t"; |
| 269 |
+ |
} |
| 270 |
+ |
dout << endl; |
| 271 |
+ |
} |
| 272 |
+ |
|
| 273 |
+ |
if(!doobserved) { |
| 274 |
+ |
dout << endl << endl << "LIMITS: " << endl; |
| 275 |
+ |
dout << "\t"; |
| 276 |
+ |
for (int irow=0;irow<jzbcuts.size();irow++) { |
| 277 |
+ |
dout << jzbcuts[irow] << "\t"; |
| 278 |
+ |
} |
| 279 |
+ |
dout << endl; |
| 280 |
+ |
for(int irow=0;irow<limits.size();irow++) { |
| 281 |
+ |
dout << limits[irow][0] << "\t"; |
| 282 |
+ |
for(int ientry=0;ientry<jzbcuts.size();ientry++) { |
| 283 |
+ |
dout << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "\t"; |
| 284 |
+ |
} |
| 285 |
+ |
dout << endl; |
| 286 |
+ |
} |
| 287 |
+ |
}//do observed |
| 288 |
+ |
|
| 289 |
+ |
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; |
| 290 |
+ |
dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t Prediction [pb]" << endl; |
| 291 |
+ |
for(int icut=0;icut<jzbcuts.size();icut++) { |
| 292 |
+ |
dout << "Region with JZB>" << jzbcuts[icut] << endl; |
| 293 |
+ |
for(int isample=0;isample<signalsamples.collection.size();isample++) { |
| 294 |
+ |
dout << limits[icut][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],1) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],1) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],1) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl; |
| 295 |
+ |
} |
| 296 |
+ |
dout << endl; |
| 297 |
+ |
} |
| 298 |
|
} |
| 299 |
|
|
| 300 |
|
void susy_scan_axis_labeling(TH2F *histo) { |
| 308 |
|
TCanvas *c3 = new TCanvas("c3","c3"); |
| 309 |
|
vector<float> binning; |
| 310 |
|
binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800); |
| 234 |
– |
/* |
| 235 |
– |
binning.push_back(50); |
| 236 |
– |
binning.push_back(75); |
| 237 |
– |
binning.push_back(100); |
| 238 |
– |
binning.push_back(150); |
| 239 |
– |
binning.push_back(200); |
| 240 |
– |
binning.push_back(500); |
| 241 |
– |
*/ |
| 311 |
|
float arrbinning[binning.size()]; |
| 312 |
|
for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i]; |
| 313 |
|
TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity); |
| 363 |
|
puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data |
| 364 |
|
stringstream saveas; |
| 365 |
|
saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23; |
| 366 |
< |
cout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl; |
| 366 |
> |
dout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl; |
| 367 |
|
// TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4")); |
| 368 |
|
// TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4")); |
| 369 |
|
// TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity); |
| 371 |
|
// signalpred->Add(signalpredro); |
| 372 |
|
// puresignal->Add(signalpred,-1);//subtracting signal contamination |
| 373 |
|
//--------------------- |
| 374 |
< |
// cout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl; |
| 374 |
> |
// dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl; |
| 375 |
|
// TH1F *puresignal = allsamples.Draw("puresignal",mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4")); |
| 376 |
|
vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile); |
| 377 |
|
if(results.size()==0) { |
| 383 |
|
exclusionmap2s->Fill(m23,m0,results[2]); |
| 384 |
|
exclusionmap3s->Fill(m23,m0,results[3]); |
| 385 |
|
delete puresignal; |
| 386 |
< |
cout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl; |
| 386 |
> |
dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl; |
| 387 |
|
} |
| 388 |
|
}//end of model scan for loop |
| 389 |
|
|
| 390 |
< |
cout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl; |
| 390 |
> |
dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl; |
| 391 |
|
c3->cd(); |
| 392 |
|
exclusionmap->Draw("CONTZ"); |
| 393 |
|
CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values"); |
