29 #include "TMultiGraph.h" 37 const double pot = 8.09e20;
43 TLatex*
prelim =
new TLatex(.9, .95,
"NOvA Simulation");
44 prelim->SetTextColor(kGray+1);
46 prelim->SetTextSize(2/30.);
47 prelim->SetTextAlign(32);
48 TGaxis::SetMaxDigits(3);
53 TString
outDir2 =
"/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/NewDataSetSyst";
63 {
"/BDTG_Nbkg_Cuts_0_3_high_optimization_Myfid_cont"},
64 {
"/BDTG_NSel_Cuts_0_3_high_optimization_Myfid_cont"},
65 {
"/BDTG_NSig_Cuts_0_3_high_optimization_Myfid_cont"},
66 {
"/BDTG_NtrueSig_Cuts_0_3_high_optimization_Myfid_cont"},
67 {
"/BDTG_Cut_WithErrorBand_Nbkg_Cuts_0_3_high_optimization_Myfid_cont"},
68 {
"/BDTGcut_WithErrorBand_NSel_Cuts_0_3_high_optimization_Myfid_cont"},
69 {
"/BDTGcut_WithErrorBand_NSig_Cuts_0_3_high_optimization_Myfid_cont"},
70 {
"/BDTGcut_WithErrorBand_NtrueSig_Cuts_0_3_high_optimization_Myfid_cont"},
71 {
"/ceff03GEV_BDTG_optimization_Myfid_cont"},
72 {
"/Purity_BDTG_Myfid_cont"},
73 {
"/FOM_BDTG_Myfid_cont"},
78 {
"_07-11-18_newpngcvn_png1_BDTGoptimization_NominalwFlux"},
83 std::vector<TH1*>& ups,
84 std::vector<TH1*>& dns,
86 float headroom,
bool newaxis)
93 else if(errCol == -1) errCol = col-7;
95 nom->SetLineColor(col);
96 nom->GetXaxis()->CenterTitle();
97 nom->GetYaxis()->CenterTitle();
99 if(newaxis) nom->Draw(
"hist ][");
101 double yMax = nom->GetBinContent(nom->GetMaximumBin());
102 cout << nom->GetMaximumBin() <<
"\t" << yMax <<
endl;
103 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
105 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
106 const double y = nom->GetBinContent(binIdx);
107 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
109 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
114 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
115 double hii = (ups[systIdx]->GetBinContent(binIdx)-
y)/y;
116 double loo = (y-dns[systIdx]->GetBinContent(binIdx))/y;
117 double hi = (ups[systIdx]->GetBinContent(binIdx)-
y);
118 double lo = (y-dns[systIdx]->GetBinContent(binIdx));
120 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
126 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
129 g->SetFillColor(errCol);
131 g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
132 nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
134 nom->Draw(
"hist ][ same");
141 std::vector<TH1*>& ups,
142 std::vector<TH1*>& dns,
143 std::vector<TH1*>& upssel,
144 std::vector<TH1*>& dnssel,
146 std::vector<TH1*>& effups,
147 std::vector<TH1*>& effdns,
148 TString outname1,TString outname2,
149 TString outname3,TString outname4,
150 TString outname5, TString outname6,TString outname7,
151 TString outname8, TString outname9,
152 Bool_t printme=
false)
155 double yMax = nom->GetBinContent(nom->GetMaximumBin());
162 TH1F* hnbkg_nsel =
new TH1F(
"hnbkg_nsel",
"; BDTG Cut; #frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}",
163 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
164 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
165 TH1F* hnbkg_nsel_num =
new TH1F(
"hnbkg_nsel_num",
"; BDTG Cut value; #sqrt{N_{bkg}}",
166 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
167 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
169 TH1F* hnbkg_nsel_den =
new TH1F(
"hnbkg_nsel_den",
";BDTG Cut value; N_{sel}-N_{bkg}",
170 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
171 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
174 TH1F* hnbkg_syst =
new TH1F(
"hnbkg_syst",
";BDTG Cut; #frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}",
175 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
176 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
179 TH1F* hnsel_nsel =
new TH1F(
"hnsel_nsel",
";BDTG Cut; #frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}",
180 nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
181 nomsel->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
183 TH1F* hnsel_nsel_num =
new TH1F(
"hnsel_nsel_num",
";BDTG Cut value; #sqrt{N_{sel}}",
184 nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
185 nomsel->GetXaxis()->GetBinUpEdge(nomsel->GetNbinsX()));
187 TH1F* hnsel_syst =
new TH1F(
"hnsel_syst",
";BDTG Cut; #frac{#deltaN^{syst}_{sig}}{N_{sel}-N_{bkg}}",
188 nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
189 nomsel->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
191 TH1F* hepsilon =
new TH1F(
"hepsilon",
";BDTG Cut; #frac{#delta#varepsilon}{#varepsilon}",
192 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
193 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
195 TH1F* hsigma =
new TH1F(
"hsigma",
";BDTG Cut; #frac{#delta#sigma}{#sigma}",
196 nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
197 nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
204 for(
int binIdx = 1; binIdx < nom->GetNbinsX()+1; ++binIdx){
206 const double y = nom ->GetBinContent(binIdx);
207 const double ysel = nomsel ->GetBinContent(binIdx);
208 const double Neff = heff->GetBinContent(binIdx);
212 double deltaNbkg2 = 0;
213 double deltaNsig2 = 0;
214 double deltaeff2 = 0;
219 double delta_eff = 0;
221 double totalpresel = 0;
227 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
228 double deltaNbkg = (ups[systIdx]->GetBinContent(binIdx)-dns[systIdx]->GetBinContent(binIdx))/2;
229 double deltaNsel = (upssel[systIdx]->GetBinContent(binIdx)-dnssel[systIdx]->GetBinContent(binIdx))/2;
230 double deltaeff = (effups[systIdx]->GetBinContent(binIdx)-effdns[systIdx]->GetBinContent(binIdx))/2;
233 deltaeff2 += deltaeff*deltaeff;
235 deltaNbkg2 += deltaNbkg*deltaNbkg;
237 deltaNsig2 += deltaNsel*deltaNsel;
241 const double s = 0.23;
242 totalbkg = y*s + deltaNbkg2;
244 double sel_bkg = ysel -
y;
245 double sel_bkg2 = sel_bkg*sel_bkg;
246 double uncbkg = deltaNbkg2/sel_bkg2;
248 double uncsel = totalpresel/sel_bkg2;
249 double Neff2 = Neff*Neff;
250 double delta_eff_new = delta_eff;
251 double unceff = deltaeff2/Neff2;
252 double uncsigma2 = uncbkg + unceff + uncsel ;
255 cout <<
"#########Bin number value for 0.90: "<< nom->GetXaxis()->FindBin(0.90);
256 cout <<
"#########Bin number value for 0.91: "<< nom->GetXaxis()->FindBin(0.91);
257 cout <<
"#########Bin number value for 0.92: "<< nom->GetXaxis()->FindBin(0.92);
258 cout <<
"#########Bin number value for 0.93: "<< nom->GetXaxis()->FindBin(0.93);
259 cout <<
"#########Bin number value for 0.94: "<< nom->GetXaxis()->FindBin(0.94);
260 cout <<
"#########Bin number value for 0.95: "<< nom->GetXaxis()->FindBin(0.95);
261 cout <<
"#########Bin number value for 0.96: "<< nom->GetXaxis()->FindBin(0.96);
262 cout <<
"#########Bin number value for 0.97: "<< nom->GetXaxis()->FindBin(0.97);
263 cout <<
"#########Bin number value for 0.98: "<< nom->GetXaxis()->FindBin(0.98);
264 cout <<
"#########Bin number value for 0.99: "<< nom->GetXaxis()->FindBin(0.99);
265 cout <<
"#########Bin number value for 1: "<< nom->GetXaxis()->FindBin(1);
267 cout <<
"*********************** bin number**********************: "<< binIdx <<
endl;
268 cout <<
"Bkg number (from bin content)"<< y <<
endl;
269 cout <<
"Total bkg (sta+sys): "<< totalbkg <<
endl;
270 cout <<
"Selected Sig(denominator): "<< sel_bkg <<
endl;
271 cout <<
"Selected Sig square(denominator): "<< sel_bkg2 <<
endl;
274 cout <<
"Sys unc in bkg - the numerator" <<
sqrt(deltaNbkg2)<<
endl;
275 cout <<
"Total unc in bkg "<< uncbkg <<
endl;
277 cout <<
"selected number (from bin content): "<< ysel <<
endl;
278 cout <<
"Total selected events (sta+sys): "<< totalpresel <<
endl;
279 cout <<
"Selected Sig(denominator): "<< sel_bkg <<
endl;
280 cout <<
"Selected Sig square(denominator): "<< sel_bkg2 <<
endl;
283 cout <<
"Total unc in selected "<< uncsel <<
endl;
284 cout <<
"Efficiency bin content (from bin content): "<< Neff <<
endl;
285 cout <<
"uncertainty efficiency bin content (from bin content): "<< delta_eff_new<<
endl;
287 cout <<
"Square of Efficiency (from bin content): "<< Neff2 <<
endl;
288 cout <<
"unc Efficiency : "<< unceff <<
endl;
289 cout <<
"sigma uncertainty" << (uncsigma2)<<endl;
290 cout <<
"square root sigma uncertainty" <<
sqrt(uncsigma2)<<
endl;
292 hnbkg_syst->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(deltaNbkg2)/
sqrt(sel_bkg2));
293 hnbkg_nsel->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(y)/
sqrt(sel_bkg2));
294 hnbkg_nsel_num->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(y));
295 hnbkg_nsel_den->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sel_bkg);
297 hnsel_syst->Fill(nomsel->GetXaxis()->GetBinCenter(binIdx),
sqrt(deltaNsig2)/
sqrt(sel_bkg2));
298 hnsel_nsel_num->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(ysel));
299 hnsel_nsel->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(ysel)/
sqrt(sel_bkg2));
302 hepsilon->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(deltaeff2)/
sqrt(Neff2));
303 hsigma->Fill(nom->GetXaxis()->GetBinCenter(binIdx),
sqrt(uncsigma2));
308 hnbkg_syst->SetLineColor(
kRed);
309 hnbkg_nsel->SetLineColor(
kRed);
310 hnbkg_nsel_num->SetLineColor(
kRed);
311 hnbkg_nsel_den->SetLineColor(
kRed);
313 hnsel_syst->SetLineColor(
kRed);
314 hnsel_nsel->SetLineColor(
kRed);
315 hnsel_nsel_num->SetLineColor(
kRed);
317 hepsilon->SetLineColor(
kRed);
318 hsigma->SetLineColor(
kRed);
320 hnbkg_syst->GetXaxis()->CenterTitle();
321 hnbkg_syst->GetYaxis()->CenterTitle();
322 hnbkg_nsel->GetXaxis()->CenterTitle();
323 hnbkg_nsel->GetYaxis()->CenterTitle();
324 hnbkg_nsel_num->GetXaxis()->CenterTitle();
325 hnbkg_nsel_num->GetYaxis()->CenterTitle();
326 hnbkg_nsel_den->GetXaxis()->CenterTitle();
327 hnbkg_nsel_den->GetYaxis()->CenterTitle();
329 hnsel_syst->GetXaxis()->CenterTitle();
330 hnsel_syst->GetYaxis()->CenterTitle();
331 hnsel_nsel->GetXaxis()->CenterTitle();
332 hnsel_nsel->GetYaxis()->CenterTitle();
333 hnsel_nsel_num->GetXaxis()->CenterTitle();
334 hnsel_nsel_num->GetYaxis()->CenterTitle();
336 hepsilon->GetXaxis()->CenterTitle();
337 hepsilon->GetYaxis()->CenterTitle();
338 hsigma->GetXaxis()->CenterTitle();
339 hsigma->GetYaxis()->CenterTitle();
342 TCanvas *
c1 =
new TCanvas(
"c1",
"c1");
343 c1->SetLeftMargin(0.15);
344 hsigma->Draw(
"hist ][");
348 TCanvas *
c2 =
new TCanvas(
"c2",
"c2");
349 c2->SetLeftMargin(0.16);
350 hnbkg_syst->Draw(
"hist ][");
351 hnbkg_syst->GetYaxis()->SetTitleOffset(1.2);
355 TCanvas *c2a =
new TCanvas(
"c2a",
"c2a");
356 c2a->SetLeftMargin(0.15);
357 hnbkg_nsel->Draw(
"hist ][");
358 hnbkg_nsel->GetYaxis()->SetTitleOffset(1.1);
363 TCanvas *c2b =
new TCanvas(
"c2b",
"c2b");
364 c2b->SetLeftMargin(0.15);
365 hnbkg_nsel_num->Draw(
"hist ][");
366 hnbkg_nsel_num->GetYaxis()->SetTitleOffset(1.1);
370 TCanvas *c2c =
new TCanvas(
"c2c",
"c2c");
371 c2c->SetLeftMargin(0.15);
372 hnbkg_nsel_den->Draw(
"hist ][");
373 hnbkg_nsel_den->GetYaxis()->SetTitleOffset(1.1);
377 TCanvas *
c3 =
new TCanvas(
"c3",
"c3");
378 c3->SetLeftMargin(0.16);
379 hnsel_syst->Draw(
"hist ][");
380 hnsel_syst->GetYaxis()->SetTitleOffset(1.2);
383 TCanvas *c3a =
new TCanvas(
"c3a",
"c3a");
384 c3a->SetLeftMargin(0.15);
385 hnsel_nsel->Draw(
"hist ][");
386 hnsel_nsel->GetYaxis()->SetTitleOffset(1.1);
390 TCanvas *c3b =
new TCanvas(
"c3b",
"c3b");
391 c3b->SetLeftMargin(0.15);
392 hnsel_nsel_num->Draw(
"hist ][");
393 hnsel_nsel_num->GetYaxis()->SetTitleOffset(1.1);
398 TCanvas *
c4 =
new TCanvas(
"c4",
"c4");
399 c4->SetLeftMargin(0.15);
400 hepsilon->Draw(
"hist ][");
406 TString
outDir =
"/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/NewDataSetSyst";
407 c1->SaveAs(outDir+
"/"+outname1+CurrentDate[0].
date.c_str()+
".pdf");
408 c2->SaveAs(outDir+
"/"+outname2+CurrentDate[0].
date.c_str()+
".pdf");
409 c2a->SaveAs(outDir+
"/"+outname3+CurrentDate[0].
date.c_str()+
".pdf");
410 c2b->SaveAs(outDir+
"/"+outname4+CurrentDate[0].
date.c_str()+
".pdf");
412 c2c->SaveAs(outDir+
"/"+outname5+CurrentDate[0].
date.c_str()+
".pdf");
413 c3->SaveAs(outDir+
"/"+outname6+CurrentDate[0].
date.c_str()+
".pdf");
414 c3a->SaveAs(outDir+
"/"+outname7+CurrentDate[0].
date.c_str()+
".pdf");
415 c3b->SaveAs(outDir+
"/"+outname8+CurrentDate[0].
date.c_str()+
".pdf");
416 c4->SaveAs(outDir+
"/"+outname9+CurrentDate[0].
date.c_str()+
".pdf");
429 std::vector<ana::Spectrum*>& supdown)
432 TH1* hup = (TH1*)hnom->Clone(
"up");
433 TH1* hdown= (TH1*)hnom->Clone(
"down");
438 double mean = hnom->GetBinContent(
ibin);
443 for(
unsigned int iuniv=0; iuniv < univs.size(); iuniv++){
444 err +=
pow(univs[iuniv]->ToTH1(
pot)->GetBinContent(
ibin) - mean, 2);
446 err /= double(univs.size());
449 hup->SetBinContent(
ibin, mean+err);
450 hdown->SetBinContent(
ibin, mean-err);
462 std::string fname =
"/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/PNG1OPTISyst_newpngcvn_Nominal.root";
464 TFile
fin(fname.c_str());
468 std::vector<std::string> taxes_labels = {
"BDTG"};
469 std::vector<std::string>
cut_labels = {
"totmc",
"NSel_Cuts",
"Nsig_fid_only",
"Nbkg_Cuts"};
472 for(
int ivar = 0; ivar < nvar; ivar++){
475 std::vector<ana::Spectrum*> calibup, calibdown;
476 std::vector<ana::Spectrum*> lightup, lightdown;
477 std::vector<ana::Spectrum*> fluxup, fluxdown;
478 std::vector<ana::Spectrum*> ckv, calib_shape;
479 std::vector<ana::Spectrum*> nominal;
480 std::vector<const ana::Spectrum*> xsecup, xsecdown;
483 for(
int icut = 0; icut < 4; icut++){
485 nominal.push_back(LoadFromFile<Spectrum>(fname,
"nominal/"+cut_labels[icut]+
"_"+taxes_labels[ivar]).
release());
497 std::vector<ana::Spectrum*> flux_univs;
498 std::vector<std::unique_ptr<ana::Spectrum> > xsec_univs;
499 for(
int iuniv = 0; iuniv < 100; iuniv++){
502 "_"+taxes_labels[ivar]+
505 flux_univs.push_back(LoadFromFile<Spectrum>(fname,
"flux_"+cut_labels[icut]+
"_"+
510 std::vector<ana::Spectrum*> fupdown;
512 fluxup.push_back(fupdown[0]);
513 fluxdown.push_back(fupdown[1]);
526 TH1* Nbkg_Cuts_xsecdown = xsecdown[3]->ToTH1(
pot);
527 TH1* Nbkg_Cuts_xsecup = xsecup[3]->ToTH1(
pot);
529 TH1* Nbkg_Cuts_fluxdown = fluxdown[3]->ToTH1(
pot);
530 TH1* Nbkg_Cuts_fluxup = fluxup[3]->ToTH1(
pot);
544 std::vector<TH1*> vector_up_Nbkg_Cuts = {xsecup[3]->ToTH1(
pot),fluxup[3]->ToTH1(
pot)};
545 std::vector<TH1*> vector_dwn_Nbkg_Cuts = {xsecdown[3]->ToTH1(
pot),fluxdown[3]->ToTH1(
pot)};
548 TH1* hnom_Nbkg_Cuts = nominal[3]->ToTH1(
pot);
550 std::cout <<
"*************Bins nominal*************: " << hnom_Nbkg_Cuts->GetXaxis()->GetNbins() <<
std::endl;
552 hnom_Nbkg_Cuts->SetTitle(
";BDTG; Background Events/8.09 #times 10^{20}POT");
555 gROOT->ProcessLine(
".! mkdir -p plots");
557 TCanvas* Prong_CVN_Gamma_ID_Nbkg_Cuts =
new TCanvas(
"Prong_CVN_Gamma_ID_Nbkg_Cuts",
"Prong_CVN_Gamma_ID_Nbkg_Cuts");
558 Prong_CVN_Gamma_ID_Nbkg_Cuts->SetLeftMargin(0.15);
562 vector_dwn_Nbkg_Cuts, -1, -1, 1.3
f,
true);
566 Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs(
outDir2+CurrentFile[0].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
580 TH1* NSel_Cuts_xsecdown = xsecdown[0]->ToTH1(
pot);
581 TH1* NSel_Cuts_xsecup = xsecup[0]->ToTH1(
pot);
583 TH1* NSel_Cuts_fluxdown = fluxdown[0]->ToTH1(
pot);
584 TH1* NSel_Cuts_fluxup = fluxup[0]->ToTH1(
pot);
598 std::vector<TH1*> vector_up_NSel_Cuts = {xsecup[0]->ToTH1(
pot),fluxup[0]->ToTH1(
pot)};
599 std::vector<TH1*> vector_dwn_NSel_Cuts = {xsecdown[0]->ToTH1(
pot),fluxdown[0]->ToTH1(
pot)};
602 TH1* hnom_NSel_Cuts = nominal[0]->ToTH1(
pot);
603 hnom_NSel_Cuts->SetTitle(
";BDTG; Selected Events/8.09 #times 10^{20}POT");
605 std::cout <<
"*************Bins nominal*************: " << hnom_NSel_Cuts->GetXaxis()->GetNbins() <<
std::endl;
606 gROOT->ProcessLine(
".! mkdir -p plots");
608 TCanvas* Prong_CVN_Gamma_ID_NSel_Cuts =
new TCanvas(
"Prong_CVN_Gamma_ID_NSel_Cuts",
"Prong_CVN_Gamma_ID_NSel_Cuts");
609 Prong_CVN_Gamma_ID_NSel_Cuts ->SetLeftMargin(0.15);
613 vector_dwn_NSel_Cuts, -1, -1, 1.3
f,
true);
617 Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs(
outDir2+CurrentFile[1].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
632 TH1* signal_xsecdown = xsecdown[1]->ToTH1(
pot);
633 TH1* signal_xsecup = xsecup[1]->ToTH1(
pot);
635 TH1* signal_fluxdown = fluxdown[1]->ToTH1(
pot);
636 TH1* signal_fluxup = fluxup[1]->ToTH1(
pot);
650 std::vector<TH1*> vector_up_signal = {xsecup[1]->ToTH1(
pot),fluxup[1]->ToTH1(
pot)};
651 std::vector<TH1*> vector_dwn_signal = {xsecdown[1]->ToTH1(
pot),fluxdown[1]->ToTH1(
pot)};
653 TH1* hnom_signal = nominal[1]->ToTH1(
pot);
654 hnom_signal->SetTitle(
";BDTG; Signal events/8.09 #times 10^{20}POT");
656 std::cout <<
"*************Bins nominal*************: " << hnom_signal->GetXaxis()->GetNbins() <<
std::endl;
657 gROOT->ProcessLine(
".! mkdir -p plots");
659 TCanvas* Prong_CVN_Gamma_ID_NSig_Cuts =
new TCanvas(
"Prong_CVN_Gamma_ID_NSig_Cuts",
"Prong_CVN_Gamma_ID_NSig_Cuts");
660 Prong_CVN_Gamma_ID_NSig_Cuts->SetLeftMargin(0.15);
664 vector_dwn_signal, -1, -1, 1.3
f,
true);
668 Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs(
outDir2+CurrentFile[2].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
684 TH1* truesignal_xsecdown = xsecdown[2]->ToTH1(
pot);
685 TH1* truesignal_xsecup = xsecup[2]->ToTH1(
pot);
687 TH1* truesignal_fluxdown = fluxdown[2]->ToTH1(
pot);
688 TH1* truesignal_fluxup = fluxup[2]->ToTH1(
pot);
702 std::vector<TH1*> vector_up_truesignal = {xsecup[2]->ToTH1(
pot),fluxup[2]->ToTH1(
pot)};
703 std::vector<TH1*> vector_dwn_truesignal = {xsecdown[2]->ToTH1(
pot),fluxdown[2]->ToTH1(
pot)};
707 TH1* hnom_truesignal = nominal[2]->ToTH1(
pot);
708 hnom_truesignal->SetTitle(
";BDTG; True signal events/8.09 #times 10^{20}POT");
711 gROOT->ProcessLine(
".! mkdir -p plots");
713 TCanvas* Prong_CVN_Gamma_ID_NtrueSig_Cut=
new TCanvas(
"Prong_CVN_Gamma_ID_NtrueSig_Cut",
"Prong_CVN_Gamma_ID_NtrueSig_Cut");
714 Prong_CVN_Gamma_ID_NtrueSig_Cut->SetLeftMargin(0.15);
717 vector_up_truesignal,
718 vector_dwn_truesignal, -1, -1, 1.3
f,
true);
722 Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs(
outDir2+CurrentFile[3].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
735 TH1* hNbkg_Cuts =
new TH1F(
"Nbkg_Cuts",
"Nbkg_Cuts", 100,-1,1);
738 hNbkg_Cuts->GetXaxis()->CenterTitle();
739 hNbkg_Cuts->GetYaxis()->CenterTitle();
741 TH1* hNbkg_Cuts_xsecdown =
new TH1F(
"xsecdown",
"xsecdown",100,-1,1);
742 TH1* hNbkg_Cuts_xsecup =
new TH1F(
"xsecup",
"xsecup",100,-1,1);
743 TH1* hNbkg_Cuts_fluxdown =
new TH1F(
"fluxdown",
"fluxdown",100,-1,1);
744 TH1* hNbkg_Cuts_fluxup =
new TH1F(
"fluxup",
"fluxup",100,-1,1);
756 TH1* hPur=
new TH1F(
"hPur",
"",100,-1,1);
758 TH1* hFom=
new TH1F(
"hPur",
"",100,-1,1);
761 for(
int n = 1;
n < hnom_Nbkg_Cuts->GetNbinsX()+1; ++
n){
764 const float nNSel_Cuts = hnom_NSel_Cuts->Integral(
n, 100);
765 const float nNbkg_Cuts = hnom_Nbkg_Cuts->Integral(
n, 100);
768 float Purity_den = nNSel_Cuts;
769 float Purity_num = nNSel_Cuts-nNbkg_Cuts;
771 float Fom_den = nNSel_Cuts;
772 float Fom_num = nNSel_Cuts-nNbkg_Cuts;
774 pur=Purity_num/nNSel_Cuts;
776 cout<<
"Total Selected:"<<nNSel_Cuts<<
endl;
777 cout<<
"Background:"<<nNbkg_Cuts<<
endl;
778 cout<<
"Purity_num:"<<nNSel_Cuts-nNbkg_Cuts<<
endl;
779 cout<<
"Purity:"<<Purity_num/nNSel_Cuts<<
endl;
783 fom= Fom_num/
sqrt(nNSel_Cuts);
785 cout<<
"Total Selected:"<<nNSel_Cuts<<
endl;
786 cout<<
"Background:"<<nNbkg_Cuts<<
endl;
787 cout<<
"Fom_num:"<<nNSel_Cuts-nNbkg_Cuts<<
endl;
791 hPur->SetBinContent(
n,pur);
792 hFom->SetBinContent(
n,fom);
793 hNbkg_Cuts->SetBinContent(
n,nNbkg_Cuts);
798 hNbkg_Cuts_xsecup->SetBinContent(
n, Nbkg_Cuts_xsecup->Integral(
n,100));
799 hNbkg_Cuts_xsecdown->SetBinContent(
n, Nbkg_Cuts_xsecdown->Integral(
n,100));
804 hNbkg_Cuts_fluxup->SetBinContent(
n, Nbkg_Cuts_fluxup->Integral(
n,100));
805 hNbkg_Cuts_fluxdown->SetBinContent(
n, Nbkg_Cuts_fluxdown->Integral(
n,100));
814 std::vector<TH1*> vector_up_Nbkg_Cuts_binbybin = {hNbkg_Cuts_xsecup,hNbkg_Cuts_fluxup};
816 std::vector<TH1*> vector_dwn_Nbkg_Cuts_binbybin = {hNbkg_Cuts_xsecdown,hNbkg_Cuts_fluxdown};
822 TCanvas* cNbkg_Cuts =
new TCanvas(
"cNbkg_Cuts",
"cNbkg_Cuts");
823 cNbkg_Cuts->SetLeftMargin(0.15);
825 hNbkg_Cuts->SetTitle(
";BDTG Cut value; Background Events/8.09 #times 10^{20} POT");
830 vector_up_Nbkg_Cuts_binbybin,
831 vector_dwn_Nbkg_Cuts_binbybin, -1, -1, 1.3
f,
true);
834 cNbkg_Cuts->SaveAs(
outDir2+CurrentFile[4].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
847 TCanvas* cpur =
new TCanvas(
"cpur",
"cpur");
848 cpur->SetLeftMargin(0.15);
849 hPur->SetLineColor(
kRed);
850 hPur->GetYaxis()->SetTitle(
"Purity");
851 hPur->GetXaxis()->SetTitle(
"BDTG Cut Value");
854 cpur->SaveAs(
outDir2+CurrentFile[9].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
864 TCanvas* cfom =
new TCanvas(
"cfom",
"cfom");
865 cfom->SetLeftMargin(0.15);
866 hFom->SetLineColor(
kRed);
867 hFom->GetYaxis()->SetTitle(
"FOM");
868 hFom->GetXaxis()->SetTitle(
"BDTG Cut Value");
871 cfom->SaveAs(
outDir2+CurrentFile[10].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
881 TH1* hNSel_Cuts =
new TH1F(
"NSel_Cuts",
"NSel_Cuts", 100,-1,1 );
883 hNSel_Cuts->GetXaxis()->CenterTitle();
884 hNSel_Cuts->GetYaxis()->CenterTitle();
887 TH1* hNSel_Cuts_xsecdown =
new TH1F(
"xsecdown",
"xsecdown",100,-1,1);
888 TH1* hNSel_Cuts_xsecup =
new TH1F(
"xsecup",
"xsecup",100,-1,1);
890 TH1* hNSel_Cuts_fluxdown =
new TH1F(
"fluxdown",
"fluxdown",100,-1,1);
891 TH1* hNSel_Cuts_fluxup =
new TH1F(
"fluxup",
"fluxup",100,-1,1);
905 for(
int n = 1;
n < hnom_NSel_Cuts->GetNbinsX() + 1; ++
n){
907 const float nsel_cuts = hnom_NSel_Cuts->Integral(
n,100);
909 hNSel_Cuts->SetBinContent(
n, nsel_cuts);
914 hNSel_Cuts_xsecup->SetBinContent(
n, NSel_Cuts_xsecup->Integral(
n,100));
915 hNSel_Cuts_xsecdown->SetBinContent(
n, NSel_Cuts_xsecdown->Integral(
n,100));
920 hNSel_Cuts_fluxup->SetBinContent(
n, NSel_Cuts_fluxup->Integral(
n,100));
921 hNSel_Cuts_fluxdown->SetBinContent(
n, NSel_Cuts_fluxdown->Integral(
n,100));
930 std::vector<TH1*> vector_up_NSel_Cuts_binbybin = {hNSel_Cuts_xsecup,hNSel_Cuts_fluxup};
934 std::vector<TH1*> vector_dwn_NSel_Cuts_binbybin = {hNSel_Cuts_xsecdown,hNSel_Cuts_fluxdown};
942 TCanvas* cNSel_Cuts =
new TCanvas(
"cNSel_Cuts",
"cNSel_Cuts");
943 cNSel_Cuts->SetLeftMargin(0.15);
945 hNSel_Cuts->SetTitle(
";BDTG Cut Value;Selected Events/8.09 #times 10^{20} POT");
950 vector_up_NSel_Cuts_binbybin,
951 vector_dwn_NSel_Cuts_binbybin, -1, -1, 1.3
f,
true);
955 cNSel_Cuts->SaveAs(
outDir2+CurrentFile[5].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
970 TH1* hsignal =
new TH1F(
"signal",
"signal", 100,-1,1);
973 hsignal->GetXaxis()->CenterTitle();
974 hsignal->GetYaxis()->CenterTitle();
977 TH1* hsignal_xsecdown =
new TH1F(
"xsecdown",
"xsecdown", 100,-1,1);
978 TH1* hsignal_xsecup =
new TH1F(
"xsecup",
"xsecup", 100,-1,1);
980 TH1* hsignal_fluxdown =
new TH1F(
"fluxdown",
"fluxdown", 100,-1,1);
981 TH1* hsignal_fluxup =
new TH1F(
"fluxup",
"fluxup", 100,-1,1);
996 for(
int n = 1;
n < hnom_signal->GetNbinsX()+1; ++
n){
998 const float nsig = hnom_signal->Integral(
n,100);
1000 hsignal->SetBinContent(
n, nsig);
1005 hsignal_xsecup->SetBinContent(
n, signal_xsecup->Integral(
n,100));
1006 hsignal_xsecdown->SetBinContent(
n, signal_xsecdown->Integral(
n,100));
1011 hsignal_fluxup->SetBinContent(
n, signal_fluxup->Integral(
n,100));
1012 hsignal_fluxdown->SetBinContent(
n, signal_fluxdown->Integral(
n,100));
1020 std::vector<TH1*> vector_up_signal_binbybin = {hsignal_xsecup,hsignal_fluxup};
1022 std::vector<TH1*> vector_dwn_signal_binbybin = {hsignal_xsecdown,hsignal_fluxdown};
1029 TCanvas* csignal =
new TCanvas(
"csignal",
"csignal");
1030 csignal->SetLeftMargin(0.15);
1032 hsignal->SetTitle(
";BDTG Cut Value;Signal Events/8.09 #times 10^{20} POT");
1035 vector_up_signal_binbybin,
1036 vector_dwn_signal_binbybin, -1, -1, 1.3
f,
true);
1040 csignal->SaveAs(
outDir2+CurrentFile[6].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
1054 TH1* htruesignal =
new TH1F(
"",
"", 100,-1,1);
1057 htruesignal->GetXaxis()->CenterTitle();
1058 htruesignal->GetYaxis()->CenterTitle();
1061 TH1* htruesignal_xsecdown =
new TH1F(
"xsecdown",
"xsecdown", 100,-1,1);
1062 TH1* htruesignal_xsecup =
new TH1F(
"xsecup",
"xsecup", 100,-1,1);
1064 TH1* htruesignal_fluxdown =
new TH1F(
"fluxdown",
"fluxdown", 100,-1,1);
1065 TH1* htruesignal_fluxup =
new TH1F(
"fluxup",
"fluxup", 100,-1,1);
1080 for(
int n = 1;
n < hnom_truesignal->GetNbinsX()+1; ++
n){
1082 const float nsig = hnom_truesignal->Integral(
n,100);
1084 htruesignal->SetBinContent(
n, nsig);
1089 htruesignal_xsecup->SetBinContent(
n, truesignal_xsecup->Integral(
n,100));
1090 htruesignal_xsecdown->SetBinContent(
n, truesignal_xsecdown->Integral(
n,100));
1095 htruesignal_fluxup->SetBinContent(
n, truesignal_fluxup->Integral(
n,100));
1096 htruesignal_fluxdown->SetBinContent(
n, truesignal_fluxdown->Integral(
n,100));
1104 std::vector<TH1*> vector_up_truesignal_binbybin = {htruesignal_xsecup,htruesignal_fluxup};
1106 std::vector<TH1*> vector_dwn_truesignal_binbybin = {htruesignal_xsecdown,htruesignal_fluxdown};
1114 TCanvas* ctruesignal =
new TCanvas(
"ctruesignal",
"ctruesignal");
1118 htruesignal->GetYaxis()->SetTitle(
"True signal Events/8.09 #times 10^{20} POT");
1119 htruesignal->SetTitleOffset(0.75,
"Y");
1121 htruesignal->GetXaxis()->SetTitle(
"BDTG Cut Value");
1123 vector_up_truesignal_binbybin,
1124 vector_dwn_truesignal_binbybin, -1, -1, 1.3
f,
true);
1128 ctruesignal->SaveAs(
outDir2+CurrentFile[7].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
1141 TH1* num_xsecdown = xsecdown[1]->ToTH1(
pot);
1142 TH1* num_xsecup = xsecup[1]->ToTH1(
pot);
1144 TH1* denom_xsecdown = xsecdown[2]->ToTH1(
pot);
1145 TH1* denom_xsecup = xsecup[2]->ToTH1(
pot);
1147 TH1* num_fluxdown = fluxdown[1]->ToTH1(
pot);
1148 TH1* num_fluxup = fluxup[1]->ToTH1(
pot);
1150 TH1* denom_fluxdown = fluxdown[2]->ToTH1(
pot);
1151 TH1* denom_fluxup = fluxup[2]->ToTH1(
pot);
1173 TH1* sig = nominal[1]->ToTH1(
pot);
1174 TH1* denome = nominal[2]->ToTH1(
pot);
1176 float totsig = denome->Integral(1, 100);
1177 float vdenom_xsecdown = denom_xsecdown->Integral(1, 100);
1178 float vdenom_xsecup = denom_xsecup->Integral(1, 100);
1179 float vdenom_fluxdown = denom_fluxdown->Integral(1, 100);
1180 float vdenom_fluxup = denom_fluxup->Integral(1, 100);
1195 TH1* nomeff =
new TH1F(
"eff_",
1196 "; BDTG Cut;Selection Efficiency",
1199 TH1* eff_calibdown =
new TH1F(
"eff_calibdown",
1200 "Efficiency vs Cut Value; BDTG ;Efficiency",
1203 TH1* eff_calibup =
new TH1F(
"eff_calibup",
1204 "Efficiency vs Cut Value;BDTG ;Efficiency",
1207 TH1* eff_lightdown =
new TH1F(
"eff_lightdown",
1208 "Efficiency vs Cut Value; BDTG ;Efficiency",
1211 TH1* eff_lightup =
new TH1F(
"eff_lightup",
1212 "Efficiency vs Cut Value; BDTG ;Efficiency",
1215 TH1* eff_ckv =
new TH1F(
"eff_ckv",
1216 "Efficiency vs Cut Value;BDTG ;Efficiency",
1219 TH1* eff_calib_shape =
new TH1F(
"eff_calibshape",
1220 "Efficiency vs Cut Value; BDTG ;Efficiency",
1223 TH1* eff_fluxdown =
new TH1F(
"eff_fluxdown",
1224 "Efficiency vs Cut Value; BDTG ;Efficiency",
1227 TH1* eff_fluxup =
new TH1F(
"eff_fluxup",
1228 "Efficiency vs Cut Value; BDTG ;Efficiency",
1231 TH1* eff_xsecdown =
new TH1F(
"eff_xsecdown",
1232 "Efficiency vs Cut Value; BDTG ;Efficiency",
1235 TH1* eff_xsecup =
new TH1F(
"eff_xsecup",
1236 "Efficiency vs Cut Value; BDTG ;Efficiency",
1245 std::cout <<
"******************Total selected events (sta+sys): "<< totsig <<
std::endl;
1247 for(
int n = 1;
n < hnom_Nbkg_Cuts->GetNbinsX() +1; ++
n){
1249 const float n1sig = sig->Integral(
n, 100);
1250 const float nxsecdown = num_xsecdown->Integral(
n, 100)/vdenom_xsecdown;
1251 const float nxsecup = num_xsecup->Integral(
n, 100)/vdenom_xsecup;
1252 const float nfluxdown = num_fluxdown->Integral(
n, 100)/vdenom_fluxdown;
1253 const float nfluxup = num_fluxup->Integral(
n, 100)/vdenom_fluxup;
1262 std::cout <<
"******************Total selected events (sta+sys): "<< totsig <<
std::endl;
1265 float neff = n1sig/(totsig);
1269 nomeff->SetBinContent(
n, neff);
1272 eff_xsecdown->SetBinContent(
n, nxsecdown);
1276 eff_xsecup->SetBinContent(
n, nxsecup);
1277 eff_fluxdown->SetBinContent(
n, nfluxdown);
1278 eff_fluxup->SetBinContent(
n, nfluxup);
1288 std::vector<TH1*> rdn = {eff_xsecdown,eff_fluxdown};
1290 std::vector<TH1*> rup = {eff_xsecup,eff_fluxup};
1295 vector_up_Nbkg_Cuts_binbybin,
1296 vector_dwn_Nbkg_Cuts_binbybin,
1297 vector_up_NSel_Cuts_binbybin,
1298 vector_dwn_NSel_Cuts_binbybin,
1300 "BDTG_sigma_Systematics_optimization_Myfid_cont",
"BDTG_Nbkg_Cuts_syst_Systematics_optimization_Myfid_cont",
"BDTG_Nbkg_Cuts_stat__Systematics_optimization_Myfid_cont",
"BDTG_NSel_Cuts_syst__Systematics_optimization_Myfid_cont",
"BDTG_NSel_Cuts_stat_Systematics_optimization_Myfid_cont",
"BDTG_delta_eff__Systematics_optimization_Myfid_cont",
"c7__Systematics_BDTG_optimization_Myfid_cont",
"c8__Systematics_BDTG_optimization_Myfid_cont",
"c9__Systematics_BDTG_optimization_Myfid_cont",
true);
1303 TCanvas* ceff =
new TCanvas(
"ceff",
"ceff");
1304 ceff->SetLeftMargin(0.15);
1309 ceff->SaveAs(
outDir2+CurrentFile[8].
name.c_str()+CurrentDate[0].
date.c_str()+
".pdf");
void GetFluxError(const Spectrum *nom, std::vector< ana::Spectrum * > univs, std::vector< ana::Spectrum * > &supdown)
TSpline3 lo("lo", xlo, ylo, 12,"0")
::xsd::cxx::tree::date< char, simple_type > date
Cuts and Vars for the 2020 FD DiF Study.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
const DateUp CurrentDate[Date]
const Color_t kTotalMCErrorBandColor
void PlotWithSystErrorBandWidth(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float headroom, bool newaxis)
std::vector< std::string > cut_labels
Representation of a spectrum in any variable, with associated POT.
const FileUp CurrentFile[File]
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TSpline3 hi("hi", xhi, yhi, 18,"0")
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void PlotDeltaSigmaWithSigma(TH1 *&nom, TH1 *&nomsel, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, std::vector< TH1 * > &upssel, std::vector< TH1 * > &dnssel, TH1 *&heff, std::vector< TH1 * > &effups, std::vector< TH1 * > &effdns, TString outname1, TString outname2, TString outname3, TString outname4, TString outname5, TString outname6, TString outname7, TString outname8, TString outname9, Bool_t printme=false)
std::vector< float > Spectrum
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
const Color_t kTotalMCColor
std::string to_string(ModuleType mt)
double Livetime() const
Seconds. For informational purposes only. No calculations use this.