10 #include "3FlavorAna/Cuts/NueCutsSecondAna.h" 57 #include "TGraphErrors.h" 59 #include "TDecompLU.h" 60 #include "TDecompSVD.h" 105 std::stringstream low_X;
106 low_X << std::fixed << std::setprecision(2) <<
107 sample->GetXaxis()->GetBinLowEdge(xbin);
108 std::stringstream hi_X;
109 hi_X << std::fixed << std::setprecision(2) <<
110 sample->GetXaxis()->GetBinUpEdge(xbin);
112 std::stringstream low_Y;
113 low_Y << std::fixed << std::setprecision(2) <<
114 sample->GetYaxis()->GetBinLowEdge(ybin);
115 std::stringstream hi_Y;
116 hi_Y << std::fixed << std::setprecision(2) <<
117 sample->GetYaxis()->GetBinUpEdge(ybin);
120 hi_X.str() +
" " + low_Y.str() +
" #leq E_{e} < " + hi_Y.str();
127 TLatex*
prelim =
new TLatex(0.23, .80, str);
128 prelim->SetTextColor(kBlack);
130 prelim->SetTextSize(0.05);
134 void fcn(Int_t &
npar, Double_t *gin, Double_t &
f, Double_t *
par, Int_t iflag){
136 double chisquare = 0;
175 FitVector2[
i][0] = 0;
179 FitVector[0][
i] = ( data - (par[1]*signal + par[0]*numu +
182 FitVector2[
i][0] = ( data -(par[1]*signal + par[0]*numu +
188 float secondterm = Chi2(0,0);
190 chisquare = secondterm;
196 Double_t &
f, Double_t *
par, Int_t iflag){
198 double chisquare = 0;
233 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);
257 float x_hi = 1.0,
float y_lo = 1.0,
260 hWeights->GetXaxis()->SetTitle(xtitle.c_str());
261 hWeights->GetYaxis()->SetTitle(ytitle.c_str());
262 hWeights->SetTitle(title.c_str());
263 hWeights->GetXaxis()->SetRangeUser(x_lo,x_hi);
264 hWeights->GetYaxis()->SetRangeUser(y_lo,y_hi);
265 hWeights->GetZaxis()->SetTitle(ztitle_wei.c_str());
266 hWeights->GetXaxis()->CenterTitle();
267 hWeights->GetYaxis()->CenterTitle();
268 hWeights->GetZaxis()->CenterTitle();
270 hErr->GetXaxis()->SetTitle(xtitle.c_str());
271 hErr->GetYaxis()->SetTitle(ytitle.c_str());
272 hErr->SetTitle(title.c_str());
273 hErr->GetXaxis()->SetRangeUser(x_lo,x_hi);
274 hErr->GetYaxis()->SetRangeUser(y_lo,y_hi);
275 hErr->GetZaxis()->SetTitle(ztitle_err.c_str());
276 hErr->GetXaxis()->CenterTitle();
277 hErr->GetYaxis()->CenterTitle();
278 hErr->GetZaxis()->CenterTitle();
281 hWeights->Draw(
"COLZ");
283 sprintf(out,
"%s%s_%s_%s.png",
"Plots/",pidName.c_str(),
284 dataName.c_str(),outnameWei.c_str());
290 sprintf(out,
"%s%s_%s_%s.png",
"Plots/",pidName.c_str(),
291 dataName.c_str(),outnameErr.c_str());
300 std::vector<TH1F*>hMCnom, std::vector<TH1F*>hMCwei,
307 hMC_nom_signal->Add(hMCnom[2],1);
308 hMC_nom_signal->Add(hMCnom[7],1);
310 hMC_nom_bkgd->Add(hMCnom[4],1);
311 hMC_nom_bkgd->Add(hMCnom[5],1);
312 hMC_nom_bkgd->Add(hMCnom[6],1);
315 hMC_wgt_signal->Add(hMCwei[2],1);
316 hMC_wgt_signal->Add(hMCwei[7],1);
318 hMC_wgt_bkgd->Add(hMCwei[2],1);
319 hMC_wgt_bkgd->Add(hMCwei[4],1);
320 hMC_wgt_bkgd->Add(hMCwei[5],1);
321 hMC_wgt_bkgd->Add(hMCwei[6],1);
323 hMC_nom_signal->SetLineColor(
kBlue+3);
324 hMC_nom_bkgd->SetLineColor(
kGreen+3);
327 hMCWeight->Add(hMCwei[1]);
328 hMCWeight->Add(hMCwei[2]);
329 hMCWeight->Add(hMCwei[3]);
330 hMCWeight->Add(hMCwei[4]);
331 hMCWeight->Add(hMCwei[5]);
332 hMCWeight->Add(hMCwei[7]);
333 hMCWeight->SetLineColor(
kRed);
335 hMCnom[0]->SetLineColor(kMagenta-3);
336 hData->SetMarkerStyle(20);
337 hData->SetTitle(caption.c_str());
338 hMCnom[0]->SetTitle(caption.c_str());
340 if(hData->GetMaximum() > hMCnom[0]->GetMaximum())
341 hMCnom[0]->
GetYaxis()->SetRangeUser(0,hData->GetMaximum()*1.3);
345 cCompare->SetLeftMargin(0.10);
346 cCompare->SetRightMargin(0.10);
349 pad1->SetFillStyle(0);
350 pad2->SetFillStyle(0);
351 pad1->SetBottomMargin(0.30);
352 pad2->SetTopMargin(1-0.30);
358 hMCWeight->SetTitle(
"");
359 hMCnom[0]->GetYaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
360 hMCnom[0]->GetXaxis()->SetLabelColor(kWhite);
361 hMCnom[0]->GetYaxis()->SetLabelSize(0.03);
362 hMCnom[0]->GetYaxis()->SetTitleSize(0.035);
363 hMCnom[0]->SetTitle(
"");
364 hMCnom[0]->Draw(
"hist");
367 hMCnom[0]->Draw(
"hist same");
368 hMCWeight->Draw(
"hist same");
370 hMC_wgt_signal->Draw(
"hist same");
371 hMC_nom_signal->Draw(
"hist same");
372 hMC_wgt_bkgd->Draw(
"hist same");
373 hMC_nom_bkgd->Draw(
"hist same");
376 leg->AddEntry(hData,
"Fake Data",
"p");
377 leg->AddEntry(hMCWeight,
"Scaled Total MC",
"l");
378 leg->AddEntry(hMC_wgt_signal,
"Scaled Signal",
"l");
379 leg->AddEntry(hMC_wgt_bkgd,
"Scaled Background",
"l");
380 leg->AddEntry(hMCnom[0],
"Total MC",
"l");
381 leg->AddEntry(hMC_nom_signal,
"Signal",
"l");
382 leg->AddEntry(hMC_nom_bkgd,
"Background",
"l");
391 hratio->GetYaxis()->SetTitle(
"MC/Data");
392 hratio->GetXaxis()->SetTitle(
"Electron Prong CVN Score");
393 hratio->GetXaxis()->CenterTitle();
394 hratio->SetMarkerColor(kBlack);
395 hratio->SetLineColor(kBlack);
396 hratio->SetMarkerStyle(20);
397 hratio->Divide(hData);
398 hratio2->Divide(hData);
399 hratio->GetYaxis()->SetLabelSize(0.03);
400 hratio->GetYaxis()->SetTitleSize(0.035);
401 hratio->GetYaxis()->SetRangeUser(0.4,1.5);
402 hratio->Draw(
"same");
403 hratio2->Draw(
"hist same");
410 sprintf(outname,
"%s%s_%s_postfit_ratio_bin_%i_%i.png",
"Plots/",
412 dataName.c_str(),xbin,ybin);
413 cCompare->SaveAs(outname);
418 std::vector<TH1F*>& ups,
419 std::vector<TH1F*>& dns,
421 float headroom,
bool newaxis)
427 else if(errCol == -1) errCol = col-9;
429 nom->SetLineColor(col);
430 nom->GetXaxis()->CenterTitle();
431 nom->GetYaxis()->CenterTitle();
433 double yMax = nom->GetBinContent(nom->GetMaximumBin());
435 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
437 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
438 const double y = nom->GetBinContent(binIdx);
439 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
441 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
446 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
447 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
448 double lo = y-dns[systIdx]->GetBinContent(binIdx);
451 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
459 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
462 g->SetFillColor(errCol);
471 std::vector<TH1F*> vPseudoExp;
474 TH2F* covariance_hist =
476 "Covariance Matrix;PID Bin;PID Bin;Covariance",
477 nom->GetXaxis()->GetNbins(),
478 0,nom->GetXaxis()->GetNbins(),
479 nom->GetXaxis()->GetNbins(),
480 0,nom->GetXaxis()->GetNbins());
483 TRandom *r0 =
new TRandom();
484 for(
int i = 0;
i < nPseudo;
i++){
487 for(
uint j = 0;
j < hSysts.size();
j++){
488 double wei = r0->Gaus(0,1);
496 hClone->Add(hSysts[
j],wei);
499 vPseudoExp.push_back(hClone);
504 for(
int j = 1;
j <= hSysts[0]->GetXaxis()->GetNbins();
j++){
505 for(
int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
510 for(
int i = 0;
i < nPseudo;
i++){
511 double val_j = vPseudoExp[
i]->GetBinContent(
j);
512 double val_k = vPseudoExp[
i]->GetBinContent(k);
513 multiple += val_j * val_k;
520 covariance_hist->SetBinContent(
j,k,covariance);
525 sprintf(out,
"Plots/covariance_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
526 TCanvas *
c1 =
new TCanvas(
"c1",
"c1");
528 c1->SetLeftMargin(0.15);
529 c1->SetRightMargin(0.15);
530 covariance_hist->GetXaxis()->CenterTitle();
531 covariance_hist->GetYaxis()->CenterTitle();
532 covariance_hist->GetZaxis()->CenterTitle();
533 covariance_hist->Draw(
"COLZ");
545 TH2F* covariance_hist =
547 "Covariance Matrix;Systematic Number;Systematic Number;Covariance",
549 0,(
int)hSysts.size(),
551 0,(
int)hSysts.size());
553 TH2F* correlation_hist =
555 "Correlation Matrix;;;Correlation",
557 0,(
int)hSysts.size(),
559 0,(
int)hSysts.size());
562 std::vector<double> vAverages;
563 for(
int j = 0;
j < (
int)hSysts.size();
j++){
565 for(
int i = 1;
i <= nominal->GetXaxis()->GetNbins();
i++){
566 avg += hSysts[
j]->GetBinContent(
i);
568 avg /= nominal->GetXaxis()->GetNbins();
569 vAverages.push_back(avg);
572 const char *
labels[6] = {
"Cal.Shp.",
"Ckv",
"Calib",
"Light",
"XSec",
"PPFX"};
576 for(
int j = 0;
j < (
int) hSysts.size();
j++){
577 for(
int k = 0; k < (
int)hSysts.size(); k++){
578 double summed_term = 0;
580 for(
int i = 1;
i <= nominal->GetXaxis()->GetNbins();
i++){
581 summed_term += (hSysts[
j]->GetBinContent(
i) - vAverages[
j]) *
582 (hSysts[k]->GetBinContent(
i) - vAverages[k]);
586 covariance_hist->SetBinContent(
j+1,k+1,covariance);
590 for(
int j = 1;
j <= covariance_hist->GetXaxis()->GetNbins();
j++){
591 for(
int k = 1; k <= covariance_hist->GetXaxis()->GetNbins(); k++){
592 double covariance = covariance_hist->GetBinContent(
j,k);
594 covariance/(
sqrt(covariance_hist->GetBinContent(
j,
j)) *
595 sqrt(covariance_hist->GetBinContent(k,k)));
596 correlation_hist->Fill(labels[
j-1],labels[k-1],correlation);
602 sprintf(out,
"Plots/covariance_systs_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
603 TCanvas *
c1 =
new TCanvas(
"c1",
"c1");
604 c1->SetLeftMargin(0.15);
605 c1->SetRightMargin(0.15);
606 covariance_hist->GetXaxis()->CenterTitle();
607 covariance_hist->GetYaxis()->CenterTitle();
608 covariance_hist->GetZaxis()->CenterTitle();
609 covariance_hist->Draw(
"COLZ");
614 correlation_hist->LabelsDeflate(
"X");
615 correlation_hist->LabelsDeflate(
"Y");
616 correlation_hist->LabelsOption(
"v");
617 sprintf(out,
"Plots/correlation_systs_%s_%i_%i.png",
618 pidName.c_str(),xbin,ybin);
619 TCanvas *
c2 =
new TCanvas(
"c2",
"c2");
620 c2->SetLeftMargin(0.15);
621 c2->SetRightMargin(0.15);
622 correlation_hist->GetXaxis()->CenterTitle();
623 correlation_hist->GetYaxis()->CenterTitle();
624 correlation_hist->GetZaxis()->CenterTitle();
625 correlation_hist->Draw(
"COLZ text");
722 std::vector<TH3F*> templates,
723 std::vector<TH3F*> systs_hists,
724 std::vector<TH3F*> systs_hists_up,
725 std::vector<TH3F*> systs_hists_down,
726 std::vector<TH3F*> analysis_templates,
730 std::vector<TH3F*> badentry;
731 if(systs_hists_up.size() != systs_hists_down.size()){
std::cerr <<
"Unequal up and down shifts. Exiting" <<
std::endl;
734 std::vector<TH3F*> hResults;
735 for(
uint i = 0;
i < analysis_templates.size();
i++){
736 hResults.push_back((TH3F*)analysis_templates[
i]->
741 TH2F* hSigWeights = (TH2F*)data->Project3D(
"yx");
743 TH2F* hSigError = (TH2F*)data->Project3D(
"yx");
745 TH2F* hNumuWeights = (TH2F*)data->Project3D(
"yx");
747 TH2F* hNumuError = (TH2F*)data->Project3D(
"yx");
749 TH2F* hNCWeights = (TH2F*)data->Project3D(
"yx");
751 TH2F* hNCError = (TH2F*)data->Project3D(
"yx");
754 for(
int xbin = 1; xbin <= data->GetXaxis()->GetNbins(); xbin++){
755 for(
int ybin = 1; ybin <= data->GetYaxis()->GetNbins(); ybin++){
768 xbin,xbin,ybin,ybin);
770 std::vector<TH1F*> hMCZ;
771 for(
uint numChns = 0; numChns < templates.size(); numChns++){
772 hMCZ.push_back((TH1F*)templates[numChns]->
774 xbin,xbin,ybin,ybin));
779 int pid_bin = hDataZ->GetXaxis()->FindBin(
signalRegion);
780 float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
781 hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
782 float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
784 if(signal_like > 100 &&
788 std::vector<TH1F*> hSystsZ;
791 for(
int systNum = 0; systNum < (
int) systs_hists.size(); systNum++){
794 xbin,xbin,ybin,ybin);
797 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
798 double err = hHolder->GetBinContent(
bin);
799 double cv = hMCZ[0]->GetBinContent(
bin);
800 hHolder->SetBinContent(
bin,
fabs(err-cv));
802 hSystsZ.push_back(hHolder);
809 for(
int systNum = 0; systNum < (
int)systs_hists_up.size(); systNum++){
810 TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
812 TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
817 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
818 double hi = hHolder->GetBinContent(
bin);
819 double lo = hHolder_dwn->GetBinContent(
bin);
820 double cv = hMCZ[0]->GetBinContent(
bin);
821 double value_hi =
fabs(hi - cv);
822 double value_lo =
fabs(lo - cv);
824 (value_hi >= value_lo)? value = value_hi : value = value_lo;
826 hHolder->SetBinContent(
bin,value);
828 hSystsZ.push_back(hHolder);
841 for(
int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
842 for(
int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
845 hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
846 + hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
857 for(
int x = 1;
x <= hDataZ->GetXaxis()->GetNbins();
x++){
858 float data = hDataZ->GetBinContent(
x);
859 float signal = hMCZ[1]->GetBinContent(
x);
860 float nuebar = hMCZ[2]->GetBinContent(
x);
861 float numu = hMCZ[3]->GetBinContent(
x);
862 float numubar = hMCZ[4]->GetBinContent(
x);
863 float nc = hMCZ[5]->GetBinContent(
x);
864 float other = hMCZ[6]->GetBinContent(
x);
865 float nuenonfiducial = hMCZ[7]->GetBinContent(
x);
867 std::cout <<
"Data: " << data <<
"\tMC: " <<
868 signal + nuebar + numu+numubar+nc+other+nuenonfiducial <<
std::endl;
871 fSigLike[
x-1] = signal + nuebar + nuenonfiducial;
886 TMinuit *gMinuit =
new TMinuit(3);
887 gMinuit->SetPrintLevel(-1);
888 gMinuit->SetFCN(
fcn);
892 gMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
894 double vstart[3] = {0.8,1.2,0.9};
895 double verror[3] = {0,0,0};
896 double vstep[3] = {0.01,0.01,0.01};
899 gMinuit->mnparm(0,
"Numu-Like Scaling",
900 vstart[0], vstep[0], 0, 0, ierflg);
901 gMinuit->mnparm(1,
"Nue-Like Scaling",
902 vstart[1], vstep[1], 0, 0, ierflg);
903 gMinuit->mnparm(2,
"NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
909 gMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
913 gMinuit->mnexcm(
"SET STR", arglist,1,ierflg);
915 gMinuit->mnexcm(
"SIMPLEX", arglist, 3, ierflg);
916 gMinuit->mnexcm(
"MIGRAD", arglist, 5, ierflg);
917 gMinuit->mnexcm(
"SEEk", arglist, 4, ierflg);
918 gMinuit->mnexcm(
"MINOS", arglist, 4, ierflg);
923 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
924 vstart[0] = fParNewStart;
925 verror[0] = fParNewErr;
926 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
927 vstart[1] = fParNewStart;
928 verror[1] = fParNewErr;
929 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
930 vstart[2] = fParNewStart;
931 verror[2] = fParNewErr;
934 bool troubleFitting =
false;
935 Double_t
fmin = 0,fedm = 0,ferrdef = 0;
936 Int_t npari =0, nparx = 0, istat = 0;
937 gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
939 Double_t corr_mat_test[3][3];
940 gMinuit->mnemat(*corr_mat_test,3);
943 double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
944 corr_mat_test[1][1]);
945 double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
946 corr_mat_test[2][2]);
947 double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
948 corr_mat_test[2][2]);
949 if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
950 fabs(corr_01) > 0.96 ||
fabs(corr_02) > 0.96 ||
fabs(corr_12) > 0.96)
958 TMinuit *jMinuit =
new TMinuit(2);
959 jMinuit->SetPrintLevel(-1);
962 jMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
964 jMinuit->mnparm(0,
"Background", vstart[0],
965 vstep[0], 0, 0, ierflg);
966 jMinuit->mnparm(1,
"Nue CC Scaling", vstart[1],
967 vstep[1], 0, 0, ierflg);
972 jMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
973 jMinuit->mnexcm(
"MIGRAD", arglist, 3, ierflg);
974 jMinuit->mnexcm(
"SEEk", arglist, 3, ierflg);
975 jMinuit->mnexcm(
"MINOS", arglist, 3, ierflg);
978 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
979 vstart[0] = fParNewStart;
980 verror[0] = fParNewErr;
981 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
982 vstart[2] = fParNewStart;
983 verror[2] = fParNewErr;
984 jMinuit->GetParameter(1,fParNewStart,fParNewErr);
985 vstart[1] = fParNewStart;
986 verror[1] = fParNewErr;
989 jMinuit->mnemat(*emat,2);
994 covariance[0][0] = emat[0][0];
995 covariance[0][1] = emat[0][1];
996 covariance[0][2] = 0;
997 covariance[1][0] = emat[1][0];
998 covariance[1][1] = emat[1][1];
999 covariance[1][2] = 0;
1000 covariance[2][0] = emat[0][0];
1001 covariance[2][1] = emat[0][1];
1002 covariance[2][2] = emat[0][0];
1005 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1006 vstart[0] = fParNewStart;
1007 verror[0] = fParNewErr;
1008 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1009 vstart[2] = fParNewStart;
1010 verror[2] = fParNewErr;
1011 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1012 vstart[1] = fParNewStart;
1013 verror[1] = fParNewErr;
1016 Double_t emat[3][3];
1017 gMinuit->mnemat(*emat,3);
1022 covariance[0][0] = emat[0][0];
1023 covariance[0][1] = emat[0][1];
1024 covariance[0][2] = emat[0][2];
1025 covariance[1][0] = emat[1][0];
1026 covariance[1][1] = emat[1][1];
1027 covariance[1][2] = emat[1][2];
1028 covariance[2][0] = emat[2][0];
1029 covariance[2][1] = emat[2][1];
1030 covariance[2][2] = emat[2][2];
1034 std::cout <<
"******** Skipping Fit For Bin " 1035 << caption <<
"\t\n\n\n\n" <<
std::endl;
1036 std::cout <<
"******** Skipping Fit For Bin " 1037 << caption <<
"\t" << signal_like <<
"\t" 1038 << signal_like/bkgd_amount <<
"\n\n\n\n" <<
std::endl;
1043 covariance[0][0] = 0;
1044 covariance[0][1] = 0;
1045 covariance[0][2] = 0;
1046 covariance[1][0] = 0;
1047 covariance[1][1] = 0;
1048 covariance[1][2] = 0;
1049 covariance[2][0] = 0;
1050 covariance[2][1] = 0;
1051 covariance[2][2] = 0;
1055 std::cout <<
"NumuCC Scaling Factor= " << numu_wei.first <<
" +- " <<
1057 std::cout <<
"NueCC Scaling Factor= " << nue_wei.first <<
" +- " <<
1059 std::cout <<
"NC Scaling Factor= " << nc_wei.first <<
" +- " <<
1063 hSigWeights->SetBinContent(xbin,ybin,nue_wei.first);
1064 hSigError->SetBinContent(xbin,ybin,nue_wei.second);
1065 hNumuWeights->SetBinContent(xbin,ybin,numu_wei.first);
1066 hNumuError->SetBinContent(xbin,ybin,numu_wei.second);
1067 hNCWeights->SetBinContent(xbin,ybin,nc_wei.first);
1068 hNCError->SetBinContent(xbin,ybin,nc_wei.second);
1072 std::vector<TH1F*> hSysts;
1073 for(
int z = 0;
z < (
int)systs_hists.size();
z++){
1074 hSysts.push_back((TH1F*)systs_hists[
z]->
1076 xbin,xbin,ybin,ybin));
1078 std::vector<TH1F*> hSystsUp;
1079 for(
int z = 0;
z < (
int)systs_hists_up.size();
z++){
1080 hSystsUp.push_back((TH1F*)systs_hists_up[
z]->
1085 std::vector<TH1F*> hSystsDown;
1086 for(
int z = 0;
z < (
int)systs_hists_down.size();
z++){
1087 hSystsDown.push_back((TH1F*)systs_hists_down[
z]->
1093 std::vector<TH1F*> shifts_up;
1094 std::vector<TH1F*> shifts_down;
1096 for(
int z = 0;
z < (
int)hSysts.size();
z++){
1097 shifts_up.push_back(hSysts[
z]);
1101 for(
int z = 0;
z < (
int)hSystsUp.size();
z++){
1102 shifts_up.push_back(hSystsUp[
z]);
1103 shifts_down.push_back(hSystsDown[z]);
1113 std::vector<TH1F*> hMCnom;
1114 for(
int z = 0;
z < (
int)hMCZ.size();
z++){
1119 hMCZ[1]->Scale(nue_wei.first);
1120 hMCZ[2]->Scale(nue_wei.first);
1121 hMCZ[3]->Scale(numu_wei.first);
1122 hMCZ[4]->Scale(numu_wei.first);
1123 hMCZ[5]->Scale(nc_wei.first);
1125 hMCZ[7]->Scale(nue_wei.first);
1129 PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
1132 for(
int numChn = 0; numChn < (
int)hResults.size(); numChn++){
1133 for(
int zbin = 1; zbin <= hResults[0]->GetZaxis()->GetNbins();zbin++){
1134 float nsig = hResults[1]->GetBinContent(xbin,ybin,zbin);
1135 float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
1136 float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
1137 float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
1138 float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
1139 float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
1140 float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
1143 float binCont = nsig * nue_wei.first;
1144 binCont += nnuebar * nue_wei.first;
1145 binCont += nnumu * numu_wei.first;
1146 binCont += nnumubar * numu_wei.first;
1147 binCont += nnc * nc_wei.first;
1149 binCont += nnf * nue_wei.first;
1150 hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
1153 pow(covariance[0][0],2)*
pow((nnumu+nnumubar),2) +
1154 pow(covariance[1][1],2)*
pow((nsig+nnuebar+nnf),2) +
1155 pow(covariance[2][2],2)*
pow((nnc),2) +
1156 2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) +
1157 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
1158 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) +
1159 nsig *
pow(nue_wei.first,2) + nnuebar *
pow(nue_wei.first,2) +
1160 nnf *
pow(nue_wei.first,2) + nnumu *
pow(numu_wei.first,2) +
1161 nnumubar *
pow(numu_wei.first,2) + nnc *
pow(nc_wei.first,2) +
1164 if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1165 covariance[1][2] == 0 && nue_wei.second == 0 &&
1166 numu_wei.second == 0 && nc_wei.second == 0){
1167 binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
1174 hResults[0]->SetBinError(xbin,ybin,zbin,
sqrt(binErr2));
1176 std::cout << zbin <<
"\t" << binCont <<
"\t" <<
1180 "\t" <<
pow(covariance[1][1],2)*
pow((nsig+nnuebar+nnf),2) <<
1181 "\t" <<
pow(covariance[2][2],2)*
pow((nnc),2) <<
"\t" <<
1182 2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) <<
1183 "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1184 "\t" << 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) <<
1185 "\t" << nsig *
pow(nue_wei.first,2) <<
"\t" <<
1186 nnuebar *
pow(nue_wei.first,2) <<
"\t" <<
1187 nnf *
pow(nue_wei.first,2) <<
"\t" <<
1188 nnumu *
pow(numu_wei.first,2) <<
"\t" <<
1189 nnumubar *
pow(numu_wei.first,2) <<
"\t" <<
1190 nnc *
pow(nc_wei.first,2) <<
"\t" << nother <<
std::endl;
1196 float binCont = nsig * nue_wei.first;
1210 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nsig,2) +
1211 nsig *
pow(nue_wei.first,2) +
1212 pow(nue_wei.second,2)*
pow(nnuebar,2) +
1213 nnuebar *
pow(nue_wei.first,2));
1215 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1216 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1219 float binCont = nnuebar * nue_wei.first;
1220 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnuebar,2) +
1221 nnuebar *
pow(nue_wei.first,2));
1222 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1223 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1226 float binCont = nnumu * numu_wei.first;
1227 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumu,2) +
1228 nnumu *
pow(numu_wei.first,2));
1229 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1230 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1233 float binCont = nnumubar * numu_wei.first;
1234 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumubar,2) +
1235 nnumubar *
pow(numu_wei.first,2));
1236 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1237 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1240 float binCont = nnc * nc_wei.first;
1241 float binErr =
sqrt(
pow(nc_wei.second,2)*
pow(nnc,2) +
1242 nnc *
pow(nc_wei.first,2));
1243 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1244 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1247 float binCont = nnf * nue_wei.first;
1248 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnf,2) +
1249 nnf *
pow(nue_wei.first,2));
1250 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1251 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1260 "nue_weights",
"nue_errors",
1261 "Reconstructed Electron Energy, E_{e} (GeV)",
1262 "Reconstructed cos #theta_{e}",
1263 "Signal-Like Template Fit Results",
1264 "Normalization",
"Uncertainty");
1266 "numu_weights",
"numu_errors",
1267 "Reconstructed Electron Energy, E_{e} (GeV)",
1268 "Reconstructed cos #theta_{e}",
1269 "Numu-Like Template Fit Results",
1270 "Normalization",
"Uncertainty");
1272 "nc_weights",
"nc_errors",
1273 "Reconstructed Electron Energy, E_{e} (GeV)",
1274 "Reconstructed cos #theta_{e}",
1275 "NC Template Fit Results",
1276 "Normalization",
"Uncertainty");
1284 std::vector<TH3F*> templates,
1285 std::vector<TH3F*> systs_hists,
1286 std::vector<TH3F*> systs_hists_up,
1287 std::vector<TH3F*> systs_hists_down,
1288 std::vector<TH3F*> analysis_templates,
1292 std::vector<TH3F*> badentry;
1293 if(systs_hists_up.size() != systs_hists_down.size()){
1297 std::vector<TH3F*> hResults;
1298 for(
uint i = 0;
i < templates.size();
i++){
1299 hResults.push_back((TH3F*)templates[
i]->
1303 TH2F* hSigWeights = (TH2F*)data->Project3D(
"yx");
1305 TH2F* hSigError = (TH2F*)data->Project3D(
"yx");
1307 TH2F* hNumuWeights = (TH2F*)data->Project3D(
"yx");
1309 TH2F* hNumuError = (TH2F*)data->Project3D(
"yx");
1311 TH2F* hNCWeights = (TH2F*)data->Project3D(
"yx");
1313 TH2F* hNCError = (TH2F*)data->Project3D(
"yx");
1316 for(
int xbin = 1; xbin <= data->GetXaxis()->GetNbins(); xbin++){
1317 for(
int ybin = 1; ybin <= data->GetYaxis()->GetNbins(); ybin++){
1330 xbin,xbin,ybin,ybin);
1332 std::vector<TH1F*> hMCZ;
1333 for(
uint numChns = 0; numChns < templates.size(); numChns++){
1334 hMCZ.push_back((TH1F*)templates[numChns]->
1336 xbin,xbin,ybin,ybin));
1341 int pid_bin = hDataZ->GetXaxis()->FindBin(
signalRegion);
1342 float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
1343 hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
1344 float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
1346 if(signal_like > 100 &&
1350 std::vector<TH1F*> hSystsZ;
1353 for(
int systNum = 0; systNum < (
int) systs_hists.size(); systNum++){
1356 xbin,xbin,ybin,ybin);
1359 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
1360 double err = hHolder->GetBinContent(
bin);
1361 double cv = hMCZ[0]->GetBinContent(
bin);
1362 hHolder->SetBinContent(
bin,
fabs(err-cv));
1364 hSystsZ.push_back(hHolder);
1392 for(
int systNum = 0; systNum < (
int)systs_hists_up.size(); systNum++){
1393 TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
1395 TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
1400 for(
int bin = 1;
bin <= hHolder->GetXaxis()->GetNbins();
bin++){
1401 double hi = hHolder->GetBinContent(
bin);
1402 double lo = hHolder_dwn->GetBinContent(
bin);
1403 double cv = hMCZ[0]->GetBinContent(
bin);
1404 double value_hi =
fabs(hi - cv);
1405 double value_lo =
fabs(lo - cv);
1407 (value_hi >= value_lo)? value = value_hi : value = value_lo;
1408 hHolder->SetBinContent(
bin,value);
1410 hSystsZ.push_back(hHolder);
1422 for(
int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
1423 for(
int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
1426 hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
1427 +hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
1438 for(
int x = 1;
x <= hDataZ->GetXaxis()->GetNbins();
x++){
1439 float data = hDataZ->GetBinContent(
x);
1440 float signal = hMCZ[1]->GetBinContent(
x);
1441 float nuebar = hMCZ[2]->GetBinContent(
x);
1442 float numu = hMCZ[3]->GetBinContent(
x);
1443 float numubar = hMCZ[4]->GetBinContent(
x);
1444 float nc = hMCZ[5]->GetBinContent(
x);
1445 float other = hMCZ[6]->GetBinContent(
x);
1446 float nuenonfiducial = hMCZ[7]->GetBinContent(
x);
1449 fSigLike[
x-1] = signal + nuebar + nuenonfiducial;
1459 TMinuit *gMinuit =
new TMinuit(3);
1460 gMinuit->SetPrintLevel(-1);
1461 gMinuit->SetFCN(
fcn);
1462 Double_t arglist[4];
1465 gMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
1467 double vstart[3] = {0.8,1.2,0.9};
1468 double verror[3] = {0,0,0};
1469 double vstep[3] = {0.01,0.01,0.01};
1472 gMinuit->mnparm(0,
"Numu-Like Scaling",
1473 vstart[0], vstep[0], 0, 0, ierflg);
1474 gMinuit->mnparm(1,
"Nue-Like Scaling",
1475 vstart[1], vstep[1], 0, 0, ierflg);
1476 gMinuit->mnparm(2,
"NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1482 gMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
1486 gMinuit->mnexcm(
"SET STR", arglist,1,ierflg);
1488 gMinuit->mnexcm(
"SIMPLEX", arglist, 3, ierflg);
1489 gMinuit->mnexcm(
"MIGRAD", arglist, 5, ierflg);
1490 gMinuit->mnexcm(
"SEEk", arglist, 4, ierflg);
1491 gMinuit->mnexcm(
"MINOS", arglist, 4, ierflg);
1494 double fParNewStart;
1496 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1497 vstart[0] = fParNewStart;
1498 verror[0] = fParNewErr;
1499 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1500 vstart[1] = fParNewStart;
1501 verror[1] = fParNewErr;
1502 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1503 vstart[2] = fParNewStart;
1504 verror[2] = fParNewErr;
1507 bool troubleFitting =
false;
1508 Double_t
fmin = 0,fedm = 0,ferrdef = 0;
1509 Int_t npari =0, nparx = 0, istat = 0;
1510 gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1512 Double_t corr_mat_test[3][3];
1513 gMinuit->mnemat(*corr_mat_test,3);
1516 double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
1517 corr_mat_test[1][1]);
1518 double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
1519 corr_mat_test[2][2]);
1520 double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
1521 corr_mat_test[2][2]);
1522 if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1523 fabs(corr_01) > 0.96 ||
fabs(corr_02) > 0.96 ||
fabs(corr_12) > 0.96)
1524 troubleFitting=
true;
1531 TMinuit *jMinuit =
new TMinuit(2);
1532 jMinuit->SetPrintLevel(-1);
1535 jMinuit->mnexcm(
"SET ERR", arglist, 1, ierflg);
1537 jMinuit->mnparm(0,
"Background", vstart[0],
1538 vstep[0], 0, 0, ierflg);
1539 jMinuit->mnparm(1,
"Nue CC Scaling", vstart[1],
1540 vstep[1], 0, 0, ierflg);
1545 jMinuit->mnexcm(
"SIMPLEX", arglist, 2, ierflg);
1546 jMinuit->mnexcm(
"MIGRAD", arglist, 3, ierflg);
1547 jMinuit->mnexcm(
"SEEk", arglist, 3, ierflg);
1548 jMinuit->mnexcm(
"MINOS", arglist, 3, ierflg);
1551 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1552 vstart[0] = fParNewStart;
1553 verror[0] = fParNewErr;
1554 jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1555 vstart[2] = fParNewStart;
1556 verror[2] = fParNewErr;
1557 jMinuit->GetParameter(1,fParNewStart,fParNewErr);
1558 vstart[1] = fParNewStart;
1559 verror[1] = fParNewErr;
1561 Double_t emat[2][2];
1562 jMinuit->mnemat(*emat,2);
1567 covariance[0][0] = emat[0][0];
1568 covariance[0][1] = emat[0][1];
1569 covariance[0][2] = 0;
1570 covariance[1][0] = emat[1][0];
1571 covariance[1][1] = emat[1][1];
1572 covariance[1][2] = 0;
1573 covariance[2][0] = emat[0][0];
1574 covariance[2][1] = emat[0][1];
1575 covariance[2][2] = emat[0][0];
1578 gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1579 vstart[0] = fParNewStart;
1580 verror[0] = fParNewErr;
1581 gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1582 vstart[2] = fParNewStart;
1583 verror[2] = fParNewErr;
1584 gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1585 vstart[1] = fParNewStart;
1586 verror[1] = fParNewErr;
1589 Double_t emat[3][3];
1590 gMinuit->mnemat(*emat,3);
1595 covariance[0][0] = emat[0][0];
1596 covariance[0][1] = emat[0][1];
1597 covariance[0][2] = emat[0][2];
1598 covariance[1][0] = emat[1][0];
1599 covariance[1][1] = emat[1][1];
1600 covariance[1][2] = emat[1][2];
1601 covariance[2][0] = emat[2][0];
1602 covariance[2][1] = emat[2][1];
1603 covariance[2][2] = emat[2][2];
1607 std::cout <<
"******** Skipping Fit For Bin " 1608 << caption <<
"\t" << signal_like <<
"\t" 1609 << signal_like/bkgd_amount <<
"\n\n\n\n" <<
std::endl;
1614 covariance[0][0] = 0;
1615 covariance[0][1] = 0;
1616 covariance[0][2] = 0;
1617 covariance[1][0] = 0;
1618 covariance[1][1] = 0;
1619 covariance[1][2] = 0;
1620 covariance[2][0] = 0;
1621 covariance[2][1] = 0;
1622 covariance[2][2] = 0;
1626 std::cout <<
"NumuCC Scaling Factor= " << numu_wei.first <<
" +- " <<
1628 std::cout <<
"NueCC Scaling Factor= " << nue_wei.first <<
" +- " <<
1630 std::cout <<
"NC Scaling Factor= " << nc_wei.first <<
" +- " <<
1634 hSigWeights->SetBinContent(xbin,ybin,nue_wei.first);
1635 hSigError->SetBinContent(xbin,ybin,nue_wei.second);
1636 hNumuWeights->SetBinContent(xbin,ybin,numu_wei.first);
1637 hNumuError->SetBinContent(xbin,ybin,numu_wei.second);
1638 hNCWeights->SetBinContent(xbin,ybin,nc_wei.first);
1639 hNCError->SetBinContent(xbin,ybin,nc_wei.second);
1643 std::vector<TH1F*> hSysts;
1644 for(
int z = 0;
z < (
int)systs_hists.size();
z++){
1645 hSysts.push_back((TH1F*)systs_hists[
z]->
1647 xbin,xbin,ybin,ybin));
1649 std::vector<TH1F*> hSystsUp;
1650 for(
int z = 0;
z < (
int)systs_hists_up.size();
z++){
1651 hSystsUp.push_back((TH1F*)systs_hists_up[
z]->
1656 std::vector<TH1F*> hSystsDown;
1657 for(
int z = 0;
z < (
int)systs_hists_down.size();
z++){
1658 hSystsDown.push_back((TH1F*)systs_hists_down[
z]->
1664 std::vector<TH1F*> shifts_up;
1665 std::vector<TH1F*> shifts_down;
1667 for(
int z = 0;
z < (
int)hSysts.size();
z++){
1668 shifts_up.push_back(hSysts[
z]);
1672 for(
int z = 0;
z < (
int)hSystsUp.size();
z++){
1673 shifts_up.push_back(hSystsUp[
z]);
1674 shifts_down.push_back(hSystsDown[z]);
1684 std::vector<TH1F*> hMCnom;
1685 for(
int z = 0;
z < (
int)hMCZ.size();
z++){
1690 hMCZ[1]->Scale(nue_wei.first);
1691 hMCZ[2]->Scale(nue_wei.first);
1692 hMCZ[3]->Scale(numu_wei.first);
1693 hMCZ[4]->Scale(numu_wei.first);
1694 hMCZ[5]->Scale(nc_wei.first);
1696 hMCZ[7]->Scale(nue_wei.first);
1700 PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
1703 for(
int numChn = 0; numChn < (
int)hResults.size(); numChn++){
1704 for(
int zbin = 1; zbin <= hResults[0]->GetZaxis()->GetNbins();zbin++){
1705 float nsig = hResults[1]->GetBinContent(xbin,ybin,zbin);
1706 float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
1707 float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
1708 float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
1709 float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
1710 float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
1711 float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
1713 float binCont = nsig * nue_wei.first;
1714 binCont += nnuebar * nue_wei.first;
1715 binCont += nnumu * numu_wei.first;
1716 binCont += nnumubar * numu_wei.first;
1717 binCont += nnc * nc_wei.first;
1719 binCont += nnf * nue_wei.first;
1720 hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
1723 pow(covariance[0][0],2)*
pow((nnumu+nnumubar),2) +
1724 pow(covariance[1][1],2)*
pow((nsig+nnuebar+nnf),2) +
1725 pow(covariance[2][2],2)*
pow((nnc),2) +
1726 2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) +
1727 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
1728 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) +
1729 nsig *
pow(nue_wei.first,2) + nnuebar *
pow(nue_wei.first,2) +
1730 nnf *
pow(nue_wei.first,2) + nnumu *
pow(numu_wei.first,2) +
1731 nnumubar *
pow(numu_wei.first,2) + nnc *
pow(nc_wei.first,2) +
1734 if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1735 covariance[1][2] == 0 && nue_wei.second == 0 &&
1736 numu_wei.second == 0 && nc_wei.second == 0){
1737 binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
1744 hResults[0]->SetBinError(xbin,ybin,zbin,
sqrt(binErr2));
1746 std::cout << zbin <<
"\t" << binCont <<
"\t" <<
1750 "\t" <<
pow(covariance[1][1],2)*
pow((nsig+nnuebar+nnf),2) <<
1751 "\t" <<
pow(covariance[2][2],2)*
pow((nnc),2) <<
"\t" <<
1752 2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) <<
1753 "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1754 "\t" << 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) <<
1755 "\t" << nsig *
pow(nue_wei.first,2) <<
"\t" <<
1756 nnuebar *
pow(nue_wei.first,2) <<
"\t" <<
1757 nnf *
pow(nue_wei.first,2) <<
"\t" <<
1758 nnumu *
pow(numu_wei.first,2) <<
"\t" <<
1759 nnumubar *
pow(numu_wei.first,2) <<
"\t" <<
1760 nnc *
pow(nc_wei.first,2) <<
"\t" << nother <<
std::endl;
1764 float binCont = nsig * nue_wei.first;
1765 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nsig,2) +
1766 nsig *
pow(nue_wei.first,2));
1767 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1768 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1771 float binCont = nnuebar * nue_wei.first;
1772 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnuebar,2) +
1773 nnuebar *
pow(nue_wei.first,2));
1774 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1775 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1778 float binCont = nnumu * numu_wei.first;
1779 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumu,2) +
1780 nnumu *
pow(numu_wei.first,2));
1781 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1782 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1785 float binCont = nnumubar * numu_wei.first;
1786 float binErr =
sqrt(
pow(numu_wei.second,2)*
pow(nnumubar,2) +
1787 nnumubar *
pow(numu_wei.first,2));
1788 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1789 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1792 float binCont = nnc * nc_wei.first;
1793 float binErr =
sqrt(
pow(nc_wei.second,2)*
pow(nnc,2) +
1794 nnc *
pow(nc_wei.first,2));
1795 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1796 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1799 float binCont = nnf * nue_wei.first;
1800 float binErr =
sqrt(
pow(nue_wei.second,2)*
pow(nnf,2) +
1801 nnf *
pow(nue_wei.first,2));
1802 hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1803 hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1812 "nue_weights",
"nue_errors",
1813 "Reconstructed Electron Energy, E_{e} (GeV)",
1814 "Reconstructed cos #theta_{e}",
1815 "Signal-Like Template Fit Results",
1816 "Normalization",
"Uncertainty");
1818 "numu_weights",
"numu_errors",
1819 "Reconstructed Electron Energy, E_{e} (GeV)",
1820 "Reconstructed cos #theta_{e}",
1821 "Numu-Like Template Fit Results",
1822 "Normalization",
"Uncertainty");
1824 "nc_weights",
"nc_errors",
1825 "Reconstructed Electron Energy, E_{e} (GeV)",
1826 "Reconstructed cos #theta_{e}",
1827 "NC Template Fit Results",
1828 "Normalization",
"Uncertainty");
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)
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
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)
correl_xv GetYaxis() -> SetDecimals()
void MakeWeightAndErrorPlot(TH2F *hWeights, TH2F *hErr, std::string pidName, std::string dataName, std::string outnameWei, std::string outnameErr, std::string ytitle, std::string xtitle, std::string title, std::string ztitle_wei, std::string ztitle_err, float x_lo=0.75, float x_hi=1.0, float y_lo=1.0, float y_hi=10.0)
std::vector< TH3F * > BinByBinTemplateFit_TemplateResults(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 CalculateCovarianceMatrix(TH1F *nom, std::vector< TH1F * > hSysts, int nPseudo, std::string pidName, int xbin, int ybin)
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)
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static const unsigned int npar
int isnan(const stan::math::var &a)
const std::vector< int > kNueCCColorDef
const XML_Char const XML_Char * data
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double stat_covMatrix[pidBins][pidBins]
TSpline3 hi("hi", xhi, yhi, 18,"0")
const XML_Char int const XML_Char * value
const std::string cv[Ncv]
TMatrixT< double > TMatrixD
std::string getCaption2D(TH3F *sample, int xbin, int ybin)
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
double covMatrix[pidBins][pidBins]
void PrintCaption(TString str)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
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 CalculateTotalCovariance(TH1F *nominal, std::vector< TH1F * > hSysts, std::string pidName, int xbin, int ybin)
std::string UniqueName()
Return a different string each time, for creating histograms.