11 #include "3FlavorAna/Cuts/NueCutsSecondAna.h" 75 #include "TGraphErrors.h" 78 #include "TDecompSVD.h" 83 const double pot = 8.09e20;
105 void fcn(Int_t &npar, Double_t *gin, Double_t &
f, Double_t *
par, Int_t iflag){
107 double chisquare = 0;
155 FitVector2[
i][0] = 0;
160 FitVector[0][
i] = ( data - (par[0]*signal + par[1]*ccbkg + par[2]*ncbkg));
161 FitVector2[
i][0] = ( data -(par[0]*signal+ par[1]*ccbkg + par[2]*ncbkg));
166 float secondterm = Chi2(0,0);
168 chisquare = secondterm;
172 void fcn2Var(Int_t &npar, Double_t *gin, Double_t &
f, Double_t *
par, Int_t iflag){
174 double chisquare = 0;
224 FitVector2[
i][0] = 0;
229 FitVector[0][
i] = ( data - (par[0]*signal + par[1]*(ccbkg+ncbkg)));
230 FitVector2[
i][0] = ( data -(par[0]*signal+ par[1]*(ccbkg+ncbkg)));
236 float secondterm = Chi2(0,0);
239 chisquare = secondterm;
249 TLatex*
prelim =
new TLatex(.55, .80, str);
250 prelim->SetTextColor(kBlack);
252 prelim->SetTextSize(0.05);
259 TLatex*
prelim =
new TLatex(.20, .75, str);
260 prelim->SetTextColor(kBlack);
262 prelim->SetTextSize(0.05);
268 std::vector<TH1F*>& ups,
269 std::vector<TH1F*>& dns,
271 float headroom,
bool newaxis)
277 else if(errCol == -1) errCol = col-9;
279 nom->SetLineColor(col);
280 nom->GetXaxis()->CenterTitle();
281 nom->GetYaxis()->CenterTitle();
284 double yMax = nom->GetBinContent(nom->GetMaximumBin());
286 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
288 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
289 const double y = nom->GetBinContent(binIdx);
290 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
292 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
297 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
298 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
299 double lo = y-dns[systIdx]->GetBinContent(binIdx);
302 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
310 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
313 g->SetFillColor(errCol);
315 for(TH1* up: ups)
delete up;
316 for(TH1* dn: dns)
delete dn;
321 std::vector<TH2F*>
mc,
322 TH1F* sig_weights, TH1F* cc_weights,
324 std::vector<TH2F*>shifts_up,
325 std::vector<TH2F*>shifts_down,
330 for(
int x = 1;
x <= data->GetYaxis()->GetNbins();
x++){
332 std::stringstream low;
333 low<<std::fixed<<std::setprecision(2)<<
334 data->GetYaxis()->GetBinLowEdge(
x);
335 std::stringstream high;
336 high<<std::fixed<<std::setprecision(2)<<
337 data->GetYaxis()->GetBinUpEdge(
x);
340 TH1F* dataproj = (TH1F*)data->ProjectionX(
"dataproj",
x,
x);
342 std::vector<TH1F*> shift_up_proj;
343 std::vector<TH1F*> shift_down_proj;
344 std::vector<TH1F*> shift_up_proj_ratio;
345 std::vector<TH1F*> shift_down_proj_ratio;
349 for(
uint j = 0;
j < shifts_up.size();
j++){
352 sprintf(name_up,
"%s_%s_%i",
"proj",
"up",
j);
353 sprintf(name_down,
"%s_%s_%i",
"proj",
"down",
j);
354 shift_up_proj.push_back((TH1F*)shifts_up[
j]->
ProjectionX(name_up,
x,
x));
355 shift_down_proj.push_back((TH1F*)shifts_down[
j]->
ProjectionX(name_down,
x,
x));
356 sprintf(name_up,
"%s_%s_%i",
"proj",
"up_ratio",
j);
357 sprintf(name_down,
"%s_%s_%i",
"proj",
"down_ratio",
j);
358 shift_up_proj_ratio.push_back((TH1F*)shifts_up[
j]->
ProjectionX(name_up,
x,
x));
359 shift_down_proj_ratio.push_back((TH1F*)shifts_down[
j]->
ProjectionX(name_down,
x,
x));
361 shift_up_proj[
j]->Rebin(2);
362 shift_down_proj[
j]->Rebin(2);
363 shift_up_proj_ratio[
j]->Rebin(2);
364 shift_down_proj_ratio[
j]->Rebin(2);
367 std::vector<TH1F*> hMC;
368 std::vector<TH1F*> hMC_wgt;
369 for(
uint i = 0;
i < mc.size();
i++){
373 sprintf(name,
"%s_%i",
"nominal",
i);
374 sprintf(name_wgt,
"%s_%i",
"weight",
i);
379 hMC_wgt[
i]->Rebin(2);
382 hMC_wgt[1]->Scale(sig_weights->GetBinContent(
x));
383 hMC_wgt[6]->Scale(cc_weights->GetBinContent(
x));
384 hMC_wgt[7]->Scale(nc_weights->GetBinContent(
x));
386 TH1F *hTotalMC_wgt = (TH1F*)hMC_wgt[1]->
Clone(
"hTotalMC_wgt");
387 hTotalMC_wgt->Add(hMC_wgt[6]);
388 hTotalMC_wgt->Add(hMC_wgt[7]);
389 hTotalMC_wgt->SetLineColor(
kBlue);
392 TH1F* hBkgNom = (TH1F*)hMC[6]->
Clone(
"hBkgNom");
393 hBkgNom->Add(hMC[7]);
394 TH1F* hBkgScl = (TH1F*)hMC_wgt[6]->
Clone(
"hBkgScl");
395 hBkgScl->Add(hMC_wgt[7]);
397 TH1F* hSigNom = (TH1F*)hMC[1]->
Clone(
"hSigNom");
398 TH1F* hSigScl = (TH1F*)hMC_wgt[1]->
Clone(
"hSigScl");
400 hSigNom->SetLineColor(30);
401 hBkgNom->SetLineColor(
kOrange-3);
403 hSigScl->SetLineColor(3);
404 hBkgScl->SetLineColor(
kOrange+8);
407 TH1F* hDataRatio = (TH1F*) hMC[0]->
Clone(
"data_ratio");
408 TH1F* hDataRatio_Wgt = (TH1F*)hTotalMC_wgt->Clone(
"data_ratio_wgt");
410 hDataRatio_Wgt->SetTitle(
"");
411 hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
413 hDataRatio->Divide(dataproj);
414 hDataRatio_Wgt->Divide(dataproj);
416 TH1F* hSignalRatio = (TH1F*)hSigScl->Clone(
"hSignalRatio");
417 hSignalRatio->Divide(hSigNom);
419 TH1F* hBkgRatio = (TH1F*)hBkgScl->Clone(
"hBkgRatio");
420 hBkgRatio->Divide(hBkgNom);
423 shift_up_proj_ratio[
ibin]->Divide(dataproj);
424 shift_down_proj_ratio[
ibin]->Divide(dataproj);
430 for(
int i = 1;
i <= dataproj->GetXaxis()->GetNbins();
i++){
431 float d_i = dataproj->GetBinContent(
i);
433 float m_i = hTotalMC_wgt->GetBinContent(
i);
434 float error = dataproj->GetBinError(
i)*dataproj->GetBinError(
i) +
435 hTotalMC_wgt->GetBinError(
i)*hTotalMC_wgt->GetBinError(
i);
437 if(error == 0)
continue;
439 chisquared += (d_i - m_i) * (d_i - m_i)/
error;
443 std::cout <<
"******ChiSq in bin :*****" <<
x <<
"\t" << chisquared <<
"\t" <<
"******And ChiSq/DoF is******* : " << chisquared/(dataproj->GetXaxis()->GetNbins() -3) <<
std::endl;
446 std::stringstream low_a;
447 low_a<<std::fixed<<std::setprecision(2)<<
448 chisquared/(dataproj->GetXaxis()->GetNbins() -3) ;
449 std::stringstream high_a;
450 high_a<<std::fixed<<std::setprecision(2)<<
451 dataproj->GetXaxis()->GetNbins() - 3;
458 hDataRatio->SetLineColor(
kRed);
459 hDataRatio_Wgt->SetLineColor(
kBlue);
461 dataproj->GetYaxis()->SetTitle(
"Events/8.09#times10^{20}POT");
462 dataproj->SetMarkerStyle(20);
465 TLegend *
leg =
new TLegend(0.4,0.65,0.89,0.89);
469 leg->AddEntry(hMC[0],
"Nominal MC",
"PL");
470 leg->AddEntry(hTotalMC_wgt,
"Scaled MC",
"PL");
471 leg->AddEntry(hSigNom,
"Nominal Sig",
"PL");
472 leg->AddEntry(hBkgNom,
"Nominal Bkg",
"PL");
473 leg->AddEntry(hSigScl,
"Scaled Sig",
"PL");
474 leg->AddEntry(hBkgScl,
"Scaled Bkg",
"PL");
475 leg->AddEntry(dataproj,
"Fake Data",
"PL");
478 TCanvas *cPID1 =
new TCanvas(
"cPID1",
"cPID1",700,700);
479 gStyle->SetPadBorderMode(0);
480 gStyle->SetFrameBorderMode(0);
481 Float_t small = 1
e-5;
482 cPID1->Divide(1,2,small,small);
484 gPad->SetBottomMargin(small);
490 dataproj->GetXaxis()->SetRangeUser(-1,0.92);
492 dataproj->SetTitle(caption.c_str());
493 dataproj->SetTitleOffset(0.08);
494 dataproj->GetXaxis()->SetTitle(
"");
496 hMC[0]->Draw(
"hist same");
497 hTotalMC_wgt->Draw(
"hist same");
498 hSigNom->Draw(
"hist same");
499 hBkgNom->Draw(
"hist same");
500 hSigScl->Draw(
"hist same");
501 hBkgScl->Draw(
"hist same");
502 TGraphAsymmErrors* asymm =
507 asymm->Draw(
"e2 same");
508 dataproj->Draw(
"same");
509 hMC[0]->Draw(
"hist same");
510 hTotalMC_wgt->Draw(
"hist same");
511 hSigNom->Draw(
"hist same");
512 hBkgNom->Draw(
"hist same");
513 hSigScl->Draw(
"hist same");
514 hBkgScl->Draw(
"hist same");
520 gPad->SetTopMargin(small);
521 hDataRatio->GetYaxis()->SetTitle(
"MC/Fake Data");
522 hDataRatio->GetXaxis()->SetRangeUser(-1,0.92);
523 hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
524 hDataRatio->GetYaxis()->CenterTitle();
525 hDataRatio->GetXaxis()->CenterTitle();
526 hDataRatio->Draw(
"hist");
527 TGraphAsymmErrors* asymm2 =
529 shift_down_proj_ratio,
531 asymm2->Draw(
"e2 same");
532 hDataRatio->Draw(
"hist same");
533 hDataRatio_Wgt->Draw(
"hist same");
536 sprintf(outname,
"%s%s_bin_%i.pdf",
"PlotsFinal/",
"bdt",
x);
537 cPID1->SaveAs(outname);
540 TLegend *legA =
new TLegend(0.4,0.6,0.89,0.89);
541 legA->AddEntry(hSignalRatio,
"Sclaed Sig/Nom. Sig",
"PL");
542 legA->AddEntry(hBkgRatio,
"Sclaed Bkg/Nom. Bkg",
"PL");
548 TH1F* data_y = (TH1F*)data->ProjectionY(
"data_y");
550 data_y->SetMarkerStyle(20);
551 data_y->GetXaxis()->SetRangeUser(0.2,1.5);
553 std::vector<TH1F*> mc_y;
554 std::vector<TH1F*> mc_y_wgt;
555 for(
uint i = 0;
i < mc.size();
i++){
558 sprintf(name,
"%s_%i",
"yproj",
i);
559 sprintf(name_wgt,
"%s_%i",
"yproj_wgt",
i);
561 mc_y_wgt.push_back((TH1F*)mc[
i]->
ProjectionY(name_wgt));
564 std::vector<TH1F*> shift_up_y;
565 std::vector<TH1F*> shift_down_y;
566 std::vector<TH1F*> shift_up_y_ratio;
567 std::vector<TH1F*> shift_down_y_ratio;
568 for(
uint i = 0;
i < shifts_up.size();
i++){
571 sprintf(name,
"%s_%i",
"syst",
i);
572 sprintf(name_wgt,
"%s_%i",
"syst_down",
i);
574 shift_up_y.push_back((TH1F*)shifts_up[
i]->
ProjectionY(name));
575 shift_down_y.push_back((TH1F*)shifts_up[
i]->
ProjectionY(name_wgt));
577 sprintf(name,
"%s_%i_ratio",
"syst",
i);
578 sprintf(name_wgt,
"%s_%i_ratio",
"syst_down",
i);
580 shift_up_y_ratio.push_back((TH1F*)shifts_up[
i]->
ProjectionY(name));
581 shift_down_y_ratio.push_back((TH1F*)shifts_up[
i]->
ProjectionY(name_wgt));
585 mc_y_wgt[1]->Multiply(sig_weights);
586 mc_y_wgt[6]->Multiply(cc_weights);
587 mc_y_wgt[7]->Multiply(nc_weights);
589 TH1F* hMCWeightedTot = (TH1F*)mc_y_wgt[1]->
Clone(
"hMCWeightedTotal");
590 hMCWeightedTot->Add(mc_y_wgt[6]);
591 hMCWeightedTot->Add(mc_y_wgt[7]);
592 hMCWeightedTot->SetLineColor(
kBlue);
594 TH1F* hDataRatio = (TH1F*)mc_y[0]->
Clone(
"hDataRatio");
595 TH1F* hDataRatio_Wgt = (TH1F*)hMCWeightedTot->Clone(
"hDataRatio_Wgt");
597 hDataRatio->SetLineColor(
kRed);
598 hDataRatio_Wgt->SetLineColor(
kBlue);
599 hDataRatio->GetYaxis()->SetRangeUser(0.5,1.5);
601 hDataRatio->Divide(data_y);
602 hDataRatio_Wgt->Divide(data_y);
604 for(
uint i = 0;
i < shifts_up.size();
i++){
605 shift_up_y_ratio[
i]->Divide(data_y);
606 shift_down_y_ratio[
i]->Divide(data_y);
609 TLegend *
leg =
new TLegend(0.45,0.5,0.80,0.80);
610 leg->AddEntry(data_y,
"Fake Data",
"PL");
611 leg->AddEntry(mc_y[0],
"Nominal MC",
"PL");
612 leg->AddEntry(hMCWeightedTot,
"Scaled MC",
"PL");
616 for(
int i = data_y->GetXaxis()->FindBin(0.21);
617 i <= data_y->GetXaxis()->FindBin(1.45);
i++){
618 float d_i = data_y->GetBinContent(
i);
619 float m_i = hMCWeightedTot->GetBinContent(
i);
620 float error = data_y->GetBinError(
i)*data_y->GetBinError(
i) +
621 hMCWeightedTot->GetBinError(
i)*hMCWeightedTot->GetBinError(
i);
623 if(error == 0)
continue;
625 chisquared += (d_i - m_i) * (d_i - m_i)/d_i;
628 std::stringstream low_a;
629 low_a<<std::fixed<<std::setprecision(2)<<
631 std::stringstream high_a;
632 high_a<<std::fixed<<std::setprecision(2)<<
642 TCanvas *cPID2 =
new TCanvas(
"cPID2",
"cPID2");
643 cPID2->SetBottomMargin(0.15);
644 gStyle->SetPadBorderMode(0);
645 gStyle->SetFrameBorderMode(0);
646 Float_t small = 1
e-5;
647 cPID2->Divide(1,2,small,small);
650 data_y->GetYaxis()->SetTitle(
"Events/8.09#times10^{20}POT");
651 data_y->GetYaxis()->SetRangeUser(0,50000);
652 data_y->GetXaxis()->SetRangeUser(0.2,1.5);
653 data_y->GetYaxis()->CenterTitle();
654 data_y->GetXaxis()->SetTitle(
"");
655 data_y->GetXaxis()->CenterTitle();
657 TGraphAsymmErrors* asymm3 =
660 asymm3->Draw(
"e2 same");
661 data_y->Draw(
"same");
662 mc_y[0]->Draw(
"hist same");
663 hMCWeightedTot->Draw(
"hist same");
669 gPad->SetTopMargin(small);
670 hDataRatio->GetYaxis()->SetTitle(
"Fake Data/MC");
671 hDataRatio->GetYaxis()->SetRangeUser(0.4,1.6);
672 hDataRatio->GetYaxis()->CenterTitle();
673 hDataRatio->GetXaxis()->CenterTitle();
674 hDataRatio->Draw(
"hist");
675 TGraphAsymmErrors* asymm4 =
679 asymm4->Draw(
"e2 same");
680 hDataRatio->Draw(
"hist same");
681 hDataRatio_Wgt->Draw(
"hist same");
683 cPID2->SaveAs(
"PlotsFinal/RecoKE.pdf");
688 TH1F* sig_weights, TH1F* cc_weights,
692 TH1F* data_y_signal = (TH1F*)data_sig->ProjectionY(
"data_y_signal");
694 std::vector<TH1F*> hMC_y;
695 std::vector<TH1F*> hMC_y_wgt;
696 for(
uint i = 0;
i < mc.size();
i++){
698 sprintf(name,
"mc_%i",
i);
701 sprintf(name_wgt,
"mc_wgt_%i",
i);
702 hMC_y_wgt.push_back((TH1F*)mc[
i]->
ProjectionY(name_wgt));
706 hMC_y_wgt[1]->Multiply(sig_weights);
707 hMC_y_wgt[6]->Multiply(cc_weights);
708 hMC_y_wgt[7]->Multiply(nc_weights);
710 for(
int i = 1;
i < hMC_y_wgt[1]->GetXaxis()->GetNbins();
i++){
712 float fSignal = hMC_y[1]->GetBinContent(
i);
713 float fSignalError = hMC_y[1]->GetBinError(
i);
714 float fTotal = hMC_y[0]->GetBinContent(
i) - fSignal;
715 float fWeight = sig_weights->GetBinContent(
i);
716 float fError = sig_weights->GetBinError(
i);
718 float total_error = (fWeight * fSignal)*
pow(
pow((fSignalError/fSignal),2) +
pow(fError/fWeight,2),0.5);
720 hMC_y_wgt[1]->SetBinError(
i,total_error);
725 data_y_signal->SetLineColor(
kOrange+7);
726 data_y_signal->SetMarkerColor(
kOrange+7);
729 hMC_y[1]->SetLineColor(
kRed);
730 hMC_y[1]->SetMarkerColor(
kRed);
733 hMC_y_wgt[1]->SetLineColor(
kBlue);
734 hMC_y_wgt[1]->SetMarkerColor(
kBlue);
737 TLegend *
leg =
new TLegend(0.40,0.50,0.80,0.90);
738 leg->AddEntry(data_y_signal,
"Fake Data True Signal",
"PL");
739 leg->AddEntry(hMC_y[1],
"Nominal Sig estimate",
"PL");
740 leg->AddEntry(hMC_y_wgt[1],
"Scaled Sig estimate",
"PL");
745 for(
int i = data_y_signal->GetXaxis()->FindBin(0.21);
746 i <= data_y_signal->GetXaxis()->FindBin(1.45);
i++){
747 float d_i = data_y_signal->GetBinContent(
i);
748 float m_i = hMC_y_wgt[1]->GetBinContent(
i);
749 float e_d_i = data_y_signal->GetBinError(
i);
750 float e_m_i = hMC_y_wgt[1]->GetBinError(
i);
751 if(d_i > max_value) max_value = d_i;
752 if(d_i == 0)
continue;
754 chisquared += (d_i - m_i) * (d_i - m_i)/d_i;
757 std::stringstream low_a;
758 low_a<<std::fixed<<std::setprecision(2)<<
760 std::stringstream high_a;
761 high_a<<std::fixed<<std::setprecision(2)<<
770 TCanvas *
c3 =
new TCanvas(
"c3",
"c3");
771 c3->SetBottomMargin(0.15);
772 gStyle->SetPadBorderMode(0);
773 gStyle->SetFrameBorderMode(0);
774 Float_t small = 1
e-5;
775 c3->Divide(1,2,small,small);
777 gPad->SetBottomMargin(small);
778 data_y_signal->GetXaxis()->SetRangeUser(0.2,1.5);
781 data_y_signal->GetYaxis()->SetTitle(
"Signal events/8.09#times10^{20}POT");
782 data_y_signal->GetYaxis()->CenterTitle();
783 data_y_signal->SetMarkerStyle(20);
784 data_y_signal->SetMarkerColor(kBlack);
785 data_y_signal->SetLineColor(kBlack);
786 data_y_signal->Draw();
787 hMC_y_wgt[1]->Draw(
"hist same");
788 hMC_y[1]->Draw(
"hist same");
789 data_y_signal->Draw(
"same");
795 gPad->SetTopMargin(small);
796 TH1F* hRatio_mc = (TH1F*)data_y_signal->Clone(
"hRatio_mc");
797 TH1F* hRatio_wgt = (TH1F*)data_y_signal->Clone(
"hRatio_wgt");
798 hRatio_mc->Divide(hMC_y[1]);
799 hRatio_wgt->Divide(hMC_y_wgt[1]);
800 hRatio_mc->SetLineColor(
kRed);
801 hRatio_wgt->SetLineColor(
kBlue);
802 hRatio_mc->GetYaxis()->SetRangeUser(0.8,1.2);
803 hRatio_mc->GetYaxis()->SetTitle(
"Data/MC");
804 hRatio_mc->GetYaxis()->CenterTitle();
805 hRatio_mc->GetXaxis()->CenterTitle();
806 hRatio_mc->Draw(
"hist");
807 hRatio_wgt->Draw(
"hist same");
809 sprintf(outname,
"%s%s_%s_%s.pdf",
"PlotsFinal/",
"signal_prediction",
810 dataname.c_str(), weightname.c_str());
811 c3->SaveAs(
"PlotsFinal/SignalPrediction.pdf");
821 "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/PPFX_universes_bdte_energy_new_bdt_bins.root";
822 TFile *flux_file =
new TFile(outname_flux.c_str(),
"read");
825 "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/GENIE_universes_bdte_energy_new_bdt_bins.root";
826 TFile *genie_file =
new TFile(outname_genie.c_str(),
"read");
830 "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Systs_samples_bdte_energy_new_bdt_bins.root";
831 TFile *systs_file =
new TFile(outname_systs.c_str(),
"read");
834 "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Nominal_MC_bdte_energy_new_bdt_bins.root";
835 TFile *nominal_file =
new TFile(outname_nominal.c_str(),
"read");
838 "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Data_bdte_energy_new_bdt_bins.root";
839 TFile *
data_file =
new TFile(outname_data.c_str(),
"read");
888 std::vector<std::unique_ptr<ana::Spectrum>>genie_spects;
889 for(
int i = 0;
i < 100;
i++){
891 sprintf(name,
"%s_%il",
"genie_universe",
i);
902 TH2F* hGenie_Nominal = (TH2F*)sGenie_Nominal->
ToTH2(
pot);
903 TH2F* hGenie_Upper = (TH2F*)sGenie_Nominal->
ToTH2(
pot);
904 TH2F* hGenie_Lower = (TH2F*)sGenie_Nominal->
ToTH2(
pot);
906 TH1F* hGenie_Upper_1d = (TH1F*)sGenie_Upper->
ToTH1(
pot);
907 TH1F* hGenie_Lower_1d = (TH1F*)sGenie_Lower->
ToTH1(
pot);
909 std::vector<float> histvalue_up_genie;
910 std::vector<float> histvalue_low_genie;
911 for(
int i = 1;
i < hGenie_Upper_1d->GetXaxis()->GetNbins();
i++){
912 histvalue_up_genie.push_back(hGenie_Upper_1d->GetBinContent(
i));
913 histvalue_low_genie.push_back(hGenie_Lower_1d->GetBinContent(
i));
917 for(
int i = 1;
i <= hGenie_Nominal->GetXaxis()->GetNbins();
i++){
918 for(
int j = 1;
j <= hGenie_Nominal->GetYaxis()->GetNbins();
j++){
919 hGenie_Upper->SetBinContent(
i,
j,histvalue_up_genie[counter]);
920 hGenie_Lower->SetBinContent(
i,
j,histvalue_low_genie[counter]);
927 std::vector<ana::Spectrum> systs_spects;
928 std::vector<ana::Spectrum> systs_spects_up;
929 std::vector<ana::Spectrum> systs_spects_down;
931 systs_spects_up.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"upcal_systs")));
932 systs_spects_down.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"dwncal_systs")));
933 systs_spects_up.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"lightup_systs")));
934 systs_spects_down.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"lightdwn_systs")));
936 systs_spects.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"calibshape_systs")));
937 systs_spects.push_back(*
Spectrum::LoadFrom(systs_file->GetDirectory(
"cherenkov_systs")));
939 std::vector<TH2F*> systs_hists;
940 std::vector<TH2F*> systs_hists_up;
941 std::vector<TH2F*> systs_hists_down;
943 for(
uint i = 0;
i < systs_spects.size();
i++){
944 systs_hists.push_back((TH2F*)systs_spects[
i].
ToTH2(
pot));
946 for(
uint i = 0;
i < systs_spects_up.size();
i++){
947 systs_hists_up.push_back((TH2F*)systs_spects_up[
i].
ToTH2(
pot));
948 systs_hists_down.push_back((TH2F*)systs_spects_down[
i].
ToTH2(
pot));
951 systs_hists_up.push_back(hGenie_Upper);
952 systs_hists_down.push_back(hGenie_Lower);
956 std::vector<ana::Spectrum> nominal_spects;
959 sprintf(name,
"%s_%s",
"template",
chns[
i].name.c_str());
963 std::vector<TH2F*> nominal_hists;
964 for(
uint i = 0;
i < nominal_spects.size();
i++){
965 nominal_hists.push_back((TH2F*)nominal_spects[
i].
ToTH2(
pot));
976 sprintf(grabData,
"%s_TotMC", dataname.c_str());
980 TH2F* data_hist = (TH2F*)data_spect.ToTH2(
pot);
982 char grabData_signal[50];
983 sprintf(grabData_signal,
"%s_Sig", dataname.c_str() );
987 TH2F* data_hist_signal = (TH2F*)data_signal_spect.ToTH2(
pot);
992 TH1F* sig_weights = (TH1F*)data_hist->ProjectionY(
"sig_weights");
993 TH1F* cc_weights = (TH1F*)data_hist->ProjectionY(
"cc_weights");
994 TH1F* nc_weights = (TH1F*)data_hist->ProjectionY(
"nc_weights");
1002 for(
int i = 1;
i <= data_hist->GetYaxis()->GetNbins();
i++){
1005 std::stringstream low;
1006 low<<std::fixed<<std::setprecision(2)<<
1007 data_hist->GetYaxis()->GetBinLowEdge(
i);
1008 std::stringstream high;
1009 high<<std::fixed<<std::setprecision(2)<<
1010 data_hist->GetYaxis()->GetBinUpEdge(
i);
1016 TH1F* hData1D = (TH1F*)data_hist->ProjectionX(
"hData1D",
i,
i);
1019 std::vector<TH1D*>hMC;
1020 for(
int j = 0;
j < (
int)nominal_hists.size();
j++){
1022 sprintf(name,
"mc_bin_%i_sample_%i",
i,
j);
1031 std::vector<TH1D*>hSysts;
1032 for(
int j = 0;
j < (
int)systs_hists.size();
j++){
1034 sprintf(name,
"syst_bin_%i_sample_%i",
i,
j);
1035 TH1D* Systs_holder = systs_hists[
j]->ProjectionX(name,
i,
i);
1036 for(
int bin = 1;
bin <= Systs_holder->GetXaxis()->GetNbins();
bin++){
1037 double err = Systs_holder->GetBinContent(
bin);
1038 double cv = hMC[0]->GetBinContent(
bin);
1039 Systs_holder->SetBinContent(
bin,
fabs(err-cv));
1041 hSysts.push_back(Systs_holder);
1042 hSysts.back()->Rebin(2);
1046 for(
int j = 0;
j < (
int) systs_hists_up.size();
j++){
1048 sprintf(name,
"syst_up_bin_%i,smaple_%i",
i,
j);
1050 sprintf(name2,
"syst_down_bin_%i,smaple_%i",
i,
j);
1051 TH1D* Systs_holder_up = systs_hists_up[
j]->ProjectionX(name,
i,
i);
1052 TH1D* Systs_holder_down = systs_hists_down[
j]->ProjectionX(name,
i,
i);
1054 for(
int bin =1;
bin <=Systs_holder_up->GetXaxis()->GetNbins();
bin++){
1055 double hi = Systs_holder_up->GetBinContent(
bin);
1056 double lo = Systs_holder_down->GetBinContent(
bin);
1057 double cv = hMC[0]->GetBinContent(
bin);
1059 Systs_holder_up->SetBinContent(
bin, value);
1062 hSysts.push_back(Systs_holder_up);
1063 hSysts.back()->Rebin(2);
1068 for(
int j = 1;
j <= hSysts[0]->GetXaxis()->GetNbins();
j++){
1069 for(
int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
1071 double multiple = 0;
1075 for(
int iSyst = 0; iSyst < (
int)hSysts.size(); iSyst++){
1076 double val_j = hMC[0]->GetBinContent(
j) +
1077 hSysts[iSyst]->GetBinContent(
j);
1078 double val_k = hMC[0]->GetBinContent(k) +
1079 hSysts[iSyst]->GetBinContent(k);
1080 multiple += val_j * val_k;
1096 for(
int j = 1;
j <= hSysts[0]->GetXaxis()->GetNbins();
j++){
1097 for(
int k = 1; k <=hSysts[0]->GetXaxis()->GetNbins(); k++){
1100 hData1D->GetBinError(
j)*hData1D->GetBinError(
j) +
1101 hMC[0]->GetBinError(k)*hMC[0]->GetBinError(k);
1114 for(
int x = 1;
x <= hData1D->GetXaxis()->GetNbins();
x++)
1116 float data = hData1D->GetBinContent(
x);
1117 float total = hMC[0]->GetBinContent(
x);
1118 float signal = hMC[1]->GetBinContent(
x);
1119 float ccbkg = hMC[6]->GetBinContent(
x);
1120 float ncbkg = hMC[7]->GetBinContent(
x);
1129 TH2F* covariance_hist =
new TH2F(
"covariance_hist",
";NC#pi^{0}ID;NC#pi^{0}ID",
1131 TH2F* stat_covariance_hist =
new TH2F(
"stat_covariance_hist",
1132 ";NC#pi^{0}ID;NC#pi^{0}ID",
1135 TH2F* correlation_hist =
new TH2F(
"correlation_hist",
1136 ";NC#pi^{0}ID;NC#pi^{0}ID",
1139 covariance_hist->GetXaxis()->CenterTitle();
1140 covariance_hist->GetYaxis()->CenterTitle();
1141 covariance_hist->GetZaxis()->CenterTitle();
1143 stat_covariance_hist->GetXaxis()->CenterTitle();
1144 stat_covariance_hist->GetYaxis()->CenterTitle();
1145 stat_covariance_hist->GetZaxis()->CenterTitle();
1147 correlation_hist->GetXaxis()->CenterTitle();
1148 correlation_hist->GetYaxis()->CenterTitle();
1149 correlation_hist->GetZaxis()->CenterTitle();
1160 correlation_hist->SetBinContent(
i+1,j+1,
1166 stat_covariance_hist->Add(covariance_hist);
1168 TCanvas *c9 =
new TCanvas(
"c9",
"c9");
1169 gStyle->SetTitleAlign(32);
1170 c9->SetRightMargin(0.15);
1171 covariance_hist->SetTitle(caption.c_str());
1172 covariance_hist->Draw(
"COLZ");
1175 sprintf(cov_name,
"%s%s_bin_%i.pdf",
"PlotsFinal/",
"error_matrix",
i);
1176 c9->SaveAs(cov_name);
1179 TCanvas *cCor =
new TCanvas(
"cCor",
"cCor");
1180 gStyle->SetTitleAlign(32);
1181 cCor->SetRightMargin(0.15);
1182 correlation_hist->GetZaxis()->SetTitle(
"Correlation");
1183 correlation_hist->SetTitle(caption.c_str());
1184 correlation_hist->Draw(
"COLZ");
1187 sprintf(cor_name,
"%s%s_bin_%i.pdf",
"PlotsFinal/",
"CorMatrix",
i);
1188 cCor->SaveAs(cor_name);
1191 TMinuit *gMinuit =
new TMinuit(3);
1192 gMinuit->SetPrintLevel(-1);
1193 gMinuit->SetFCN(
fcn);
1194 Double_t arglist[4];
1197 gMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
1200 double vstart[3] = {1.0,1.0,1.0};
1201 double verror[3] = {0,0,0};
1202 double vstep[3] = {0.1,0.1,0.1};
1204 gMinuit->mnparm(0,
"Signal scaling", vstart[0], vstep[0], 0, 0, ierflg);
1205 gMinuit->mnparm(1,
"CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1206 gMinuit->mnparm(2,
"NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1211 gMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
1213 double fParNewStart;
1214 double fParamNewErr;
1215 for(
int ll = 0; ll < 3; ll++){
1216 gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1217 vstart[ll] = fParNewStart;
1221 gMinuit->mnexcm(
"MIGRAD", arglist, 3, ierflg);
1222 gMinuit->mnexcm(
"SEEk", arglist, 3, ierflg);
1223 for(
int ll = 0; ll < 3; ll++){
1224 gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1225 vstart[ll] = fParNewStart;
1230 TMinuit *hMinuit =
new TMinuit(3);
1232 hMinuit->SetFCN(
fcn);
1234 hMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
1236 hMinuit->mnparm(0,
"Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
1237 hMinuit->mnparm(1,
"CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1238 hMinuit->mnparm(2,
"NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1245 hMinuit->mnexcm(
"SET STR",arglist,1,ierflg);
1247 hMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
1248 hMinuit->mnexcm(
"MIGRAD", arglist, 5, ierflg);
1249 hMinuit->mnexcm(
"SEEk", arglist, 4, ierflg);
1250 hMinuit->mnexcm(
"MINOS", arglist, 4, ierflg);
1252 for(
int ll = 0; ll < 3; ll++){
1253 hMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1254 vstart[ll] = fParNewStart;
1255 verror[ll] = fParamNewErr;
1258 bool troubleFitting =
false;
1262 Double_t errdef = 0;
1266 hMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat);
1268 troubleFitting =
true;
1270 if(!troubleFitting){
1271 hMinuit->Command(
"SHOw FCNvalue");
1276 TMinuit *jMinuit =
new TMinuit(2);
1277 jMinuit->SetPrintLevel(-1);
1282 jMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
1284 jMinuit->mnparm(0,
"Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
1285 jMinuit->mnparm(1,
"TotBkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1290 jMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
1291 jMinuit->mnexcm(
"MIGRAD", arglist, 3, ierflg);
1292 jMinuit->mnexcm(
"SEEk", arglist, 3, ierflg);
1294 jMinuit->GetParameter(0,fParNewStart,fParamNewErr);
1295 vstart[0] = fParNewStart;
1296 verror[0] = fParamNewErr;
1299 jMinuit->GetParameter(1,fParNewStart,fParamNewErr);
1300 vstart[1] = fParNewStart;
1301 verror[1] = fParamNewErr;
1302 vstart[2] = fParNewStart;
1303 verror[2] = fParamNewErr;
1305 jMinuit->Command(
"SHOw FCNvalue");
1314 hMinuit->GetParameter(0,fParamVal,fParamErr);
1317 hMinuit->GetParameter(1,fParamVal_2,fParamErr_2);
1320 hMinuit->GetParameter(2,fParamVal_3,fParamErr_3);
1323 std::cout <<
"Signal Scaling Factor= " << vstart[0] <<
" +- " <<
1325 sig_weights->SetBinContent(
i, vstart[0]);
1326 sig_weights->SetBinError(
i, verror[0]);
1328 std::cout <<
"CC-Bkg Scaling Factor= " << vstart[1] <<
" +- " <<
1330 cc_weights->SetBinContent(
i, vstart[1]);
1331 cc_weights->SetBinError(
i, verror[1]);
1332 std::cout <<
"NC-Bkg Factor= " << vstart[2] <<
" +- " <<
1334 nc_weights->SetBinContent(
i, vstart[2]);
1335 nc_weights->SetBinError(
i, verror[2]);
1345 MakePlots(data_hist,nominal_hists,sig_weights,cc_weights,
1346 nc_weights, systs_hists_up,systs_hists_down,dataname,weightname);
1348 sig_weights, cc_weights,
1349 nc_weights,dataname,weightname);
1352 TCanvas *cWeights =
new TCanvas(
"cWeights",
"cWeights");
1353 sig_weights->GetYaxis()->SetTitle(
"Signal Weight");
1354 sig_weights->GetXaxis()->SetTitleOffset(0.75);
1355 sig_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1356 sig_weights->GetYaxis()->SetRangeUser(0,1.5);
1357 sig_weights->GetXaxis()->CenterTitle();
1358 sig_weights->GetYaxis()->CenterTitle();
1359 sig_weights->Draw(
"e1");
1363 cWeights->SaveAs(
"PlotsFinal/SigWts.pdf");
1368 TCanvas *cWeights3 =
new TCanvas(
"cWeights3",
"cWeights3");
1369 cc_weights->GetYaxis()->SetTitle(
"CCBkg Weight");
1370 cc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1371 cc_weights->GetYaxis()->SetRangeUser(0,1.5);
1372 cc_weights->GetXaxis()->SetTitleOffset(0.75);
1373 cc_weights->GetXaxis()->CenterTitle();
1374 cc_weights->GetYaxis()->CenterTitle();
1375 cc_weights->Draw(
"e1");
1378 cWeights3->SaveAs(
"PlotsFinal/CCBkgWts.pdf");
1383 TCanvas *cWeights5 =
new TCanvas(
"cWeights5",
"cWeights5");
1384 nc_weights->GetYaxis()->SetTitle(
"NC Weight");
1385 nc_weights->GetXaxis()->SetTitleOffset(0.75);
1386 nc_weights->GetXaxis()->CenterTitle();
1387 nc_weights->GetYaxis()->CenterTitle();
1388 nc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1389 nc_weights->GetYaxis()->SetRangeUser(0,1.5);
1390 nc_weights->Draw(
"e1");
1392 cWeights5->SaveAs(
"PlotsFinal/NCWtd.pdf");
const std::vector< int > kNueCCColorDef
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
TSpline3 lo("lo", xlo, ylo, 12,"0")
float covMatrix[xbins][xbins]
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Cuts and Vars for the 2020 FD DiF Study.
fvar< T > fabs(const fvar< T > &x)
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
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 SelDef chns[knumchns]
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
void PrintChiSq2(TString str)
void CenterTitles(TH1 *histo)
Representation of a spectrum in any variable, with associated POT.
float stat_covMatrix[xbins][xbins]
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
const XML_Char const XML_Char * data
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TSpline3 hi("hi", xhi, yhi, 18,"0")
const XML_Char int const XML_Char * value
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void CompareToTrueSignal(TH2F *data, TH2F *data_sig, std::vector< TH2F * > mc, TH1F *sig_weights, TH1F *cc_weights, TH1F *nc_weights, std::string dataname, std::string weightname)
const std::string cv[Ncv]
TMatrixT< double > TMatrixD
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
void PrintChiSq(TString str)
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
void MakePlots(TH2F *data, std::vector< TH2F * > mc, TH1F *sig_weights, TH1F *cc_weights, TH1F *nc_weights, std::vector< TH2F * >shifts_up, std::vector< TH2F * >shifts_down, std::string dataname, std::string weightname)
TGraphAsymmErrors * PlotWithSystErrorBand_mine(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
const Spectrum * Nominal() const
return the nominal universe