56 #include "TGraphErrors.h" 58 #include "TDecompLU.h" 59 #include "TDecompSVD.h" 104 std::stringstream low_X;
105 low_X << std::fixed << std::setprecision(2) <<
106 sample->GetXaxis()->GetBinLowEdge(xbin);
107 std::stringstream hi_X;
108 hi_X << std::fixed << std::setprecision(2) <<
109 sample->GetXaxis()->GetBinUpEdge(xbin);
111 std::stringstream low_Y;
112 low_Y << std::fixed << std::setprecision(2) <<
113 sample->GetYaxis()->GetBinLowEdge(ybin);
114 std::stringstream hi_Y;
115 hi_Y << std::fixed << std::setprecision(2) <<
116 sample->GetYaxis()->GetBinUpEdge(ybin);
119 hi_X.str() +
" " + low_Y.str() +
" #leq E_{e} < " + hi_Y.str();
126 TLatex*
prelim =
new TLatex(0.23, .80, str);
127 prelim->SetTextColor(kBlack);
129 prelim->SetTextSize(0.05);
133 void fcn(Int_t &npar, Double_t *gin, Double_t &
f, Double_t *
par, Int_t iflag){
135 double chisquare = 0;
174 FitVector2[
i][0] = 0;
178 FitVector[0][
i] = ( data - (par[1]*signal + par[0]*numu +
181 FitVector2[
i][0] = ( data -(par[1]*signal + par[0]*numu +
187 float secondterm = Chi2(0,0);
189 chisquare = secondterm;
195 Double_t &
f, Double_t *
par, Int_t iflag){
232 FitVector2[
i][0] = 0;
237 FitVector[0][
i] = ( data - (par[1]*signal + par[0]*numu +
240 FitVector2[
i][0] = ( data -(par[1]*signal + par[0]*numu +
246 float secondterm = Chi2(0,0);
256 std::vector<TH1F*>hMCnom, std::vector<TH1F*>hMCwei,
263 hMC_nom_signal->Add(hMCnom[2],1);
264 hMC_nom_signal->Add(hMCnom[7],1);
266 hMC_nom_bkgd->Add(hMCnom[4],1);
267 hMC_nom_bkgd->Add(hMCnom[5],1);
268 hMC_nom_bkgd->Add(hMCnom[6],1);
271 hMC_wgt_signal->Add(hMCwei[2],1);
272 hMC_wgt_signal->Add(hMCwei[7],1);
274 hMC_wgt_bkgd->Add(hMCwei[2],1);
275 hMC_wgt_bkgd->Add(hMCwei[4],1);
276 hMC_wgt_bkgd->Add(hMCwei[5],1);
277 hMC_wgt_bkgd->Add(hMCwei[6],1);
279 hMC_nom_signal->SetLineColor(
kBlue+3);
280 hMC_nom_bkgd->SetLineColor(
kGreen+3);
283 hMCWeight->Add(hMCwei[1]);
284 hMCWeight->Add(hMCwei[2]);
285 hMCWeight->Add(hMCwei[3]);
286 hMCWeight->Add(hMCwei[4]);
287 hMCWeight->Add(hMCwei[5]);
288 hMCWeight->Add(hMCwei[7]);
289 hMCWeight->SetLineColor(
kRed);
291 hMCnom[0]->SetLineColor(kMagenta-3);
292 hData->SetMarkerStyle(20);
293 hData->SetTitle(caption.c_str());
294 hMCnom[0]->SetTitle(caption.c_str());
296 if(hData->GetMaximum() > hMCnom[0]->GetMaximum())
297 hMCnom[0]->
GetYaxis()->SetRangeUser(0,hData->GetMaximum()*1.3);
301 cCompare->SetLeftMargin(0.10);
302 cCompare->SetRightMargin(0.10);
305 pad1->SetFillStyle(0);
306 pad2->SetFillStyle(0);
307 pad1->SetBottomMargin(0.30);
308 pad2->SetTopMargin(1-0.30);
314 hMCWeight->SetTitle(
"");
315 hMCnom[0]->GetYaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
316 hMCnom[0]->GetXaxis()->SetLabelColor(kWhite);
317 hMCnom[0]->GetYaxis()->SetLabelSize(0.03);
318 hMCnom[0]->GetYaxis()->SetTitleSize(0.035);
319 hMCnom[0]->SetTitle(
"");
320 hMCnom[0]->Draw(
"hist");
323 hMCnom[0]->Draw(
"hist same");
324 hMCWeight->Draw(
"hist same");
326 hMC_wgt_signal->Draw(
"hist same");
327 hMC_nom_signal->Draw(
"hist same");
328 hMC_wgt_bkgd->Draw(
"hist same");
329 hMC_nom_bkgd->Draw(
"hist same");
332 leg->AddEntry(hData,
"Fake Data",
"p");
333 leg->AddEntry(hMCWeight,
"Scaled Total MC",
"l");
334 leg->AddEntry(hMC_wgt_signal,
"Scaled Signal",
"l");
335 leg->AddEntry(hMC_wgt_bkgd,
"Scaled Background",
"l");
336 leg->AddEntry(hMCnom[0],
"Total MC",
"l");
337 leg->AddEntry(hMC_nom_signal,
"Signal",
"l");
338 leg->AddEntry(hMC_nom_bkgd,
"Background",
"l");
347 hratio->GetYaxis()->SetTitle(
"MC/Data");
348 hratio->GetXaxis()->SetTitle(
"Electron Prong CVN Score");
349 hratio->GetXaxis()->CenterTitle();
350 hratio->SetMarkerColor(kBlack);
351 hratio->SetLineColor(kBlack);
352 hratio->SetMarkerStyle(20);
353 hratio->Divide(hData);
354 hratio2->Divide(hData);
355 hratio->GetYaxis()->SetLabelSize(0.03);
356 hratio->GetYaxis()->SetTitleSize(0.035);
357 hratio->GetYaxis()->SetRangeUser(0.4,1.5);
358 hratio->Draw(
"same");
359 hratio2->Draw(
"hist same");
366 sprintf(outname,
"%s%s_%s_postfit_ratio_bin_%i_%i.png",
"Plots/",
368 dataName.c_str(),xbin,ybin);
369 cCompare->SaveAs(outname);
374 std::vector<TH1F*>& ups,
375 std::vector<TH1F*>& dns,
377 float headroom,
bool newaxis)
383 else if(errCol == -1) errCol = col-9;
385 nom->SetLineColor(col);
386 nom->GetXaxis()->CenterTitle();
387 nom->GetYaxis()->CenterTitle();
391 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
393 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
394 const double y = nom->GetBinContent(binIdx);
395 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
397 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
402 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
403 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
404 double lo = y-dns[systIdx]->GetBinContent(binIdx);
407 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
408 if(lo < 0 && hi > 0){
413 if(lo > 0 && hi < 0) {
425 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
428 g->SetFillColor(errCol);
437 std::vector<TH1F*> vPseudoExp;
439 TH2F* covariance_hist =
441 "Covariance Matrix;PID Bin;PID Bin;Covariance",
442 nom->GetXaxis()->GetNbins(),
443 0,nom->GetXaxis()->GetNbins(),
444 nom->GetXaxis()->GetNbins(),
445 0,nom->GetXaxis()->GetNbins());
448 TRandom *r0 =
new TRandom();
449 for(
int i = 0;
i < nPseudo;
i++){
452 for(
uint j = 0;
j < hSysts.size();
j++){
453 double wei = r0->Gaus(0,1);
455 hClone->Add(hSysts[
j],wei);
458 vPseudoExp.push_back(hClone);
463 for(
int j = 1;
j <= hSysts[0]->GetXaxis()->GetNbins();
j++){
464 for(
int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
469 for(
int i = 0;
i < nPseudo;
i++){
470 double val_j = vPseudoExp[
i]->GetBinContent(
j);
471 double val_k = vPseudoExp[
i]->GetBinContent(k);
472 multiple += val_j * val_k;
479 covariance_hist->SetBinContent(
j,k,covariance);
484 sprintf(out,
"Plots/covariance_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
485 TCanvas *
c1 =
new TCanvas(
"c1",
"c1");
486 c1->SetLeftMargin(0.15);
487 c1->SetRightMargin(0.15);
488 covariance_hist->GetXaxis()->CenterTitle();
489 covariance_hist->GetYaxis()->CenterTitle();
490 covariance_hist->GetZaxis()->CenterTitle();
491 covariance_hist->Draw(
"COLZ");
501 std::vector<TH1F*> vPseudoExp;
505 TH2F* covariance_hist =
507 "Covariance Matrix;Systematic Number;Systematic Number;Covariance",
509 0,(
int)hSysts.size(),
511 0,(
int)hSysts.size());
513 TH2F* correlation_hist =
515 "Correlation Matrix;;;Correlation",
517 0,(
int)hSysts.size(),
519 0,(
int)hSysts.size());
522 std::vector<double> vAverages;
523 for(
int j = 0;
j < (
int)hSysts.size();
j++){
525 for(
int i = 0;
i < nom->GetXaxis()->GetNbins();
i++){
526 avg += hSysts[
j]->GetBinContent(
i);
528 avg /= nom->GetXaxis()->GetNbins();
529 vAverages.push_back(avg);
532 for(
int i = 0;
i < (
int)vAverages.size();
i++){
536 const char *
labels[6] = {
"Cal.Shp.",
"Ckv",
"Calib",
"Light",
"XSec",
"PPFX"};
540 for(
int j = 0;
j < (
int) hSysts.size();
j++){
541 for(
int k = 0; k < (
int)hSysts.size(); k++){
542 double summed_term = 0;
544 for(
int i = 0;
i < nom->GetXaxis()->GetNbins();
i++){
545 summed_term += (hSysts[
j]->GetBinContent(
i) - vAverages[
j]) *
546 (hSysts[k]->GetBinContent(
i) - vAverages[k]);
550 covariance_hist->SetBinContent(
j+1,k+1,covariance);
555 for(
int j = 1;
j <= covariance_hist->GetXaxis()->GetNbins();
j++){
556 for(
int k = 1; k <= covariance_hist->GetXaxis()->GetNbins(); k++){
557 double covariance = covariance_hist->GetBinContent(
j,k);
559 covariance/(
sqrt(covariance_hist->GetBinContent(
j,
j)) *
560 sqrt(covariance_hist->GetBinContent(k,k)));
561 correlation_hist->Fill(labels[
j-1],labels[k-1],correlation);
567 sprintf(out,
"Plots/covariance_systs_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
568 TCanvas *
c1 =
new TCanvas(
"c1",
"c1");
569 c1->SetLeftMargin(0.15);
570 c1->SetRightMargin(0.15);
571 covariance_hist->GetXaxis()->CenterTitle();
572 covariance_hist->GetYaxis()->CenterTitle();
573 covariance_hist->GetZaxis()->CenterTitle();
574 covariance_hist->Draw(
"COLZ");
579 correlation_hist->LabelsDeflate(
"X");
580 correlation_hist->LabelsDeflate(
"Y");
581 correlation_hist->LabelsOption(
"v");
582 sprintf(out,
"Plots/correlation_systs_%s_%i_%i.png",
583 pidName.c_str(),xbin,ybin);
584 TCanvas *
c2 =
new TCanvas(
"c2",
"c2");
585 c2->SetLeftMargin(0.15);
586 c2->SetRightMargin(0.15);
587 correlation_hist->GetXaxis()->CenterTitle();
588 correlation_hist->GetYaxis()->CenterTitle();
589 correlation_hist->GetZaxis()->CenterTitle();
590 correlation_hist->Draw(
"COLZ text");
598 std::vector<TH3F*> templates,
599 std::vector<TH3F*> systs_hists,
600 std::vector<TH3F*> systs_hists_up,
601 std::vector<TH3F*> systs_hists_down,
602 std::vector<TH3F*> analysis_templates,
606 std::vector<TH3F*> badentry;
607 if(systs_hists_up.size() != systs_hists_down.size()){
611 std::vector<TH3F*> hResults;
612 for(
uint i = 0;
i < analysis_templates.size();
i++){
613 hResults.push_back((TH3F*)analysis_templates[
i]->
618 for(
int xbin = 0; xbin <= data->GetXaxis()->GetNbins(); xbin++){
619 for(
int ybin = 0; ybin <= data->GetYaxis()->GetNbins(); ybin++){
632 xbin,xbin,ybin,ybin);
634 std::vector<TH1F*> hMCZ;
635 for(
uint numChns = 0; numChns < templates.size(); numChns++){
636 hMCZ.push_back((TH1F*)templates[numChns]->
638 xbin,xbin,ybin,ybin));
643 int pid_bin = hDataZ->GetXaxis()->FindBin(
signalRegion);
644 float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
645 hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
646 float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
648 if(signal_like > 100 &&
652 std::vector<TH1F*> hSystsZ;
655 for(
int systNum = 0; systNum < (
int) systs_hists.size(); systNum++){
658 xbin,xbin,ybin,ybin);
661 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
662 double err = hHolder->GetBinContent(
bin);
663 double cv = hMCZ[0]->GetBinContent(
bin);
664 hHolder->SetBinContent(
bin,
fabs(err-cv));
666 hSystsZ.push_back(hHolder);
673 for(
int systNum = 0; systNum < (
int)systs_hists_up.size(); systNum++){
674 TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
676 TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
681 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
682 double hi = hHolder->GetBinContent(
bin);
683 double lo = hHolder_dwn->GetBinContent(
bin);
684 double cv = hMCZ[0]->GetBinContent(
bin);
685 double val_hi =
fabs(hi - cv);
686 double val_lo =
fabs(lo - cv);
688 double value = (val_hi > val_lo)? val_hi : val_lo;
689 hHolder->SetBinContent(
bin,value);
691 hSystsZ.push_back(hHolder);
702 for(
int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
703 for(
int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
706 hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
707 +hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
718 for(
int x = 1;
x <= hDataZ->GetXaxis()->GetNbins();
x++){
719 float data = hDataZ->GetBinContent(
x);
721 float nue = hMCZ[1]->GetBinContent(
x);
722 float nuebar = hMCZ[2]->GetBinContent(
x);
723 float numu = hMCZ[3]->GetBinContent(
x);
724 float numubar = hMCZ[4]->GetBinContent(
x);
725 float nc = hMCZ[5]->GetBinContent(
x);
726 float other = hMCZ[6]->GetBinContent(
x);
727 float nuenonfiducial = hMCZ[7]->GetBinContent(
x);
730 fSigLike[
x-1] = nue + nuebar + nuenonfiducial;
740 TMinuit *gMinuit =
new TMinuit(3);
741 gMinuit->SetPrintLevel(-1);
742 gMinuit->SetFCN(
fcn);
746 gMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
748 double vstart[3] = {0.8,1.2,0.9};
749 double verror[3] = {0,0,0};
750 double vstep[3] = {0.01,0.01,0.01};
753 gMinuit->mnparm(0,
"Numu-Like Scaling",
754 vstart[0], vstep[0], 0, 0, ierflg);
755 gMinuit->mnparm(1,
"Nue-Like Scaling",
756 vstart[1], vstep[1], 0, 0, ierflg);
757 gMinuit->mnparm(2,
"NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
763 gMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
767 gMinuit->mnexcm(
"SET STR", arglist,1,ierflg);
769 gMinuit->mnexcm(
"SIMPLEX", arglist, 3, ierflg);
770 gMinuit->mnexcm(
"MIGRAD", arglist, 5, ierflg);
771 gMinuit->mnexcm(
"SEEk", arglist, 4, ierflg);
772 gMinuit->mnexcm(
"MINOS", arglist, 4, ierflg);
777 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
778 vstart[0] = fParNewStart;
779 verror[0] = fParNewErr;
780 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
781 vstart[1] = fParNewStart;
782 verror[1] = fParNewErr;
783 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
784 vstart[2] = fParNewStart;
785 verror[2] = fParNewErr;
788 bool troubleFitting =
false;
789 Double_t
fmin = 0,fedm = 0,ferrdef = 0;
790 Int_t npari =0, nparx = 0, istat = 0;
791 gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
793 Double_t corr_mat_test[3][3];
794 gMinuit->mnemat(*corr_mat_test,3);
797 double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
798 corr_mat_test[1][1]);
799 double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
800 corr_mat_test[2][2]);
801 double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
802 corr_mat_test[2][2]);
803 if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
804 fabs(corr_01) > 0.96 ||
fabs(corr_02) > 0.96 ||
fabs(corr_12) > 0.96)
812 TMinuit *jMinuit =
new TMinuit(2);
813 jMinuit->SetPrintLevel(-1);
816 jMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
818 jMinuit->mnparm(0,
"Background", vstart[0],
819 vstep[0], 0, 0, ierflg);
820 jMinuit->mnparm(1,
"Nue CC Scaling", vstart[1],
821 vstep[1], 0, 0, ierflg);
826 jMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
827 jMinuit->mnexcm(
"MIGRAD", arglist, 3, ierflg);
828 jMinuit->mnexcm(
"SEEk", arglist, 3, ierflg);
829 jMinuit->mnexcm(
"MINOS", arglist, 3, ierflg);
832 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
833 vstart[0] = fParNewStart;
834 verror[0] = fParNewErr;
835 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
836 vstart[2] = fParNewStart;
837 verror[2] = fParNewErr;
838 jMinuit->GetParameter(1,fParNewStart,fParNewErr);
839 vstart[1] = fParNewStart;
840 verror[1] = fParNewErr;
843 jMinuit->mnemat(*emat,2);
848 covariance[0][0] = emat[0][0];
849 covariance[0][1] = emat[0][1];
850 covariance[0][2] = 0;
851 covariance[1][0] = emat[1][0];
852 covariance[1][1] = emat[1][1];
853 covariance[1][2] = 0;
854 covariance[2][0] = emat[0][0];
855 covariance[2][1] = emat[0][1];
856 covariance[2][2] = emat[0][0];
859 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
860 vstart[0] = fParNewStart;
861 verror[0] = fParNewErr;
862 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
863 vstart[2] = fParNewStart;
864 verror[2] = fParNewErr;
865 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
866 vstart[1] = fParNewStart;
867 verror[1] = fParNewErr;
871 gMinuit->mnemat(*emat,3);
876 covariance[0][0] = emat[0][0];
877 covariance[0][1] = emat[0][1];
878 covariance[0][2] = emat[0][2];
879 covariance[1][0] = emat[1][0];
880 covariance[1][1] = emat[1][1];
881 covariance[1][2] = emat[1][2];
882 covariance[2][0] = emat[2][0];
883 covariance[2][1] = emat[2][1];
884 covariance[2][2] = emat[2][2];
888 std::cout <<
"******** Skipping Fit For Bin " 893 covariance[0][0] = 0;
894 covariance[0][1] = 0;
895 covariance[0][2] = 0;
896 covariance[1][0] = 0;
897 covariance[1][1] = 0;
898 covariance[1][2] = 0;
899 covariance[2][0] = 0;
900 covariance[2][1] = 0;
901 covariance[2][2] = 0;
905 std::cout <<
"NumuCC Scaling Factor= " << numu_wei.first <<
" +- " <<
907 std::cout <<
"NueCC Scaling Factor= " << nue_wei.first <<
" +- " <<
909 std::cout <<
"NC Scaling Factor= " << nc_wei.first <<
" +- " <<
914 std::vector<TH1F*> hSysts;
915 for(
int z = 0;
z < (
int)systs_hists.size();
z++){
916 hSysts.push_back((TH1F*)systs_hists[
z]->
918 xbin,xbin,ybin,ybin));
920 std::vector<TH1F*> hSystsUp;
921 for(
int z = 0;
z < (
int)systs_hists_up.size();
z++){
922 hSystsUp.push_back((TH1F*)systs_hists_up[
z]->
927 std::vector<TH1F*> hSystsDown;
928 for(
int z = 0;
z < (
int)systs_hists_down.size();
z++){
929 hSystsDown.push_back((TH1F*)systs_hists_down[
z]->
935 std::vector<TH1F*> shifts_up;
936 std::vector<TH1F*> shifts_down;
938 for(
int z = 0;
z < (
int)hSysts.size();
z++){
939 shifts_up.push_back(hSysts[
z]);
943 for(
int z = 0;
z < (
int)hSystsUp.size();
z++){
944 shifts_up.push_back(hSystsUp[
z]);
945 shifts_down.push_back(hSystsDown[z]);
955 std::vector<TH1F*> hMCnom;
956 for(
int z = 0;
z < (
int)hMCZ.size();
z++){
961 hMCZ[1]->Scale(nue_wei.first);
962 hMCZ[2]->Scale(nue_wei.first);
963 hMCZ[3]->Scale(numu_wei.first);
964 hMCZ[4]->Scale(numu_wei.first);
965 hMCZ[5]->Scale(nc_wei.first);
967 hMCZ[7]->Scale(nue_wei.first);
971 PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
974 for(
int numChn = 0; numChn < (
int)hResults.size(); numChn++){
975 for(
int zbin = 0; zbin < hResults[0]->GetZaxis()->GetNbins();zbin++){
976 float nnue = hResults[1]->GetBinContent(xbin,ybin,zbin);
977 float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
978 float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
979 float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
980 float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
981 float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
982 float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
984 float binCont = nnue * nue_wei.first;
985 binCont += nnuebar * nue_wei.first;
986 binCont += nnumu * numu_wei.first;
987 binCont += nnumubar * numu_wei.first;
988 binCont += nnc * nc_wei.first;
990 binCont += nnf * nue_wei.first;
991 hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
994 pow(covariance[0][0],2)*
pow((nnumu+nnumubar),2) +
995 pow(covariance[1][1],2)*
pow((nnue+nnuebar+nnf),2) +
996 pow(covariance[2][2],2)*
pow((nnc),2) +
997 2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) +
998 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
999 2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) +
1000 nnue *
pow(nue_wei.first,2) + nnuebar *
pow(nue_wei.first,2) +
1001 nnf *
pow(nue_wei.first,2) + nnumu *
pow(numu_wei.first,2) +
1002 nnumubar *
pow(numu_wei.first,2) + nnc *
pow(nc_wei.first,2) +
1005 if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1006 covariance[1][2] == 0 && nue_wei.second == 0 &&
1007 numu_wei.second == 0 && nc_wei.second == 0){
1008 binErr2 = nnue + nnuebar + nnf + nnumu + nnumubar +
1015 hResults[0]->SetBinError(xbin,ybin,zbin,
sqrt(binErr2));
1017 std::cout << zbin <<
"\t" << binCont <<
"\t" <<
1021 "\t" <<
pow(covariance[1][1],2)*
pow((nnue+nnuebar+nnf),2) <<
1022 "\t" <<
pow(covariance[2][2],2)*
pow((nnc),2) <<
"\t" <<
1023 2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) <<
1024 "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1025 "\t" << 2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) <<
1026 "\t" << nnue *
pow(nue_wei.first,2) <<
"\t" <<
1027 nnuebar *
pow(nue_wei.first,2) <<
"\t" <<
1028 nnf *
pow(nue_wei.first,2) <<
"\t" <<
1029 nnumu *
pow(numu_wei.first,2) <<
"\t" <<
1030 nnumubar *
pow(numu_wei.first,2) <<
"\t" <<
1031 nnc *
pow(nc_wei.first,2) <<
"\t" << nother <<
std::endl;
1035 float binCont = nnue * nue_wei.first;
1036 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnue,2) +
1037 nnue *
pow(nue_wei.first,2));
1038 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1039 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1042 float binCont = nnuebar * nue_wei.first;
1043 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnuebar,2) +
1044 nnuebar *
pow(nue_wei.first,2) +
1045 pow(nue_wei.second,2)*
pow(nnue,2) +
1046 nnue *
pow(nue_wei.first,2));
1047 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1048 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1051 float binCont = nnumu * numu_wei.first;
1052 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumu,2) +
1053 nnumu *
pow(numu_wei.first,2));
1054 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1055 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1058 float binCont = nnumubar * numu_wei.first;
1059 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumubar,2) +
1060 nnumubar *
pow(numu_wei.first,2));
1061 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1062 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1065 float binCont = nnc * nc_wei.first;
1066 float binErr =
sqrt(
pow(nc_wei.second,2)*
pow(nnc,2) +
1067 nnc *
pow(nc_wei.first,2));
1068 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1069 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1072 float binCont = nnf * nue_wei.first;
1073 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnf,2) +
1074 nnf *
pow(nue_wei.first,2));
1075 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1076 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
T max(const caf::Proxy< T > &a, T b)
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
TSpline3 lo("lo", xlo, ylo, 12,"0")
Cuts and Vars for the 2020 FD DiF Study.
fvar< T > fabs(const fvar< T > &x)
correl_xv GetYaxis() -> SetDecimals()
double stat_covMatrix[pidBins][pidBins]
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
int isnan(const stan::math::var &a)
double covMatrix[pidBins][pidBins]
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
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
std::string getCaption2D(TH3F *sample, int xbin, int ybin)
const std::string cv[Ncv]
void CalculateCovarianceMatrix(TH1F *nom, std::vector< TH1F * > hSysts, int nPseudo, std::string pidName, int xbin, int ybin)
const std::vector< int > kNueCCColorDef
TMatrixT< double > TMatrixD
std::vector< TH3F * > BinByBinTemplateFit(TH3F *data, std::vector< TH3F * > templates, std::vector< TH3F * > systs_hists, std::vector< TH3F * > systs_hists_up, std::vector< TH3F * > systs_hists_down, std::vector< TH3F * > analysis_templates, std::string pidName, std::string varName, std::string dataName)
void PlotFitResults(TH1F *hData, std::vector< TH1F * >hMCnom, std::vector< TH1F * >hMCwei, TGraphAsymmErrors *g, std::string pidName, std::string dataName, std::string caption, int xbin, int ybin)
void PrintCaption(TString str)
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
void CalculateTotalCovariance(TH1F *nom, std::vector< TH1F * > hSysts, std::string pidName, int xbin, int ybin)
std::string UniqueName()
Return a different string each time, for creating histograms.
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)