14 #include "TRatioPlot.h" 20 #include "TLegendEntry.h" 53 -1.0, 0.50, 0.56, 0.62, 0.68,
54 0.74, 0.80, 0.85, 0.88, 0.91,
55 0.94, 0.96, 0.98, 0.99, 1.0
59 0.50, 0.56, 0.62, 0.68,
60 0.74, 0.80, 0.85, 0.88, 0.91,
61 0.94, 0.96, 0.98, 0.99, 1.0
66 0.5, 0.6, 0.7, 0.8, 0.9,
67 1.0, 1.1, 1.2, 1.3, 1.4,
68 1.5, 1.6, 1.7, 1.8, 1.9,
69 2.0, 2.1, 2.2, 2.3, 2.4,
74 0.5, 0.6, 0.7, 0.8, 0.9,
75 1.0, 1.1, 1.2, 1.3, 1.4,
76 1.5, 1.6, 1.7, 1.8, 1.9,
77 2.0, 2.1, 2.2, 2.3, 2.4,
85 1.0, 1.25, 1.50, 1.75,
86 2.0, 2.25, 2.50, 2.75,
87 3.0, 3.25, 3.50, 3.75,
88 4.0, 4.25, 4.50, 4.75,
94 0.05, 0.10, 0.15, 0.20, 0.25,
95 0.30, 0.40, 0.50, 0.60, 0.75,
96 0.90, 1.10, 1.40, 1.80, 2.8
107 string xsec2DTitle =
"#frac{d^{2} #sigma}{d cos(#theta_{#mu}) d T_{#mu}} (cm^{2} / GeV / nucleon #times 10^{39})";
108 string xsecENuTitle =
"#sigma / E_{#nu} (cm^{2} / GeV / nucleon #times 10^{39})";
109 string xsecQ2Title =
"#frac{d #sigma}{d q^{2}} (cm^{2} / GeV^{2} / nucleon #times 10^{39})";
116 int NbinX = histOut->GetXaxis()->GetNbins();
117 int NbinY = histOut->GetYaxis()->GetNbins();
120 for(
int ibinx = 1; ibinx <= NbinX; ++ibinx){
121 for(
int ibiny = 1; ibiny <= NbinY; ++ibiny){
125 histOut->SetBinContent(ibinx, ibiny, histIn->GetBinContent(counter));
126 histOut->SetBinError(ibinx, ibiny, histIn->GetBinError(counter));
130 histOut->GetXaxis()->SetRangeUser(0.5, 1.0);
131 histOut->GetYaxis()->SetRangeUser(0.5, 2.5);
139 for (
int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
140 histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
141 histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
150 for (
int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
151 histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
152 histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
160 TH1I * arr =
new TH1I(
"covTo2DBin",
"", nbins, 1, nbins + 1);
164 for(
int ibinx = 1; ibinx <= NbinX; ++ibinx){
165 for(
int ibiny = 1; ibiny <= NbinY; ++ibiny){
168 arr->SetBinContent(bincounter, hist->GetBin(ibinx, ibiny));
173 cout <<
"Counter=" << bincounter <<
" NBins=" << nbins <<
endl;
174 assert(bincounter == nbins);
181 TH2I * binmap =
new TH2I(
"hist2DBinToCov",
"",
187 for(
int ibinx = 1; ibinx <= NbinX; ++ibinx){
188 for(
int ibiny = 1; ibiny <= NbinY; ++ibiny){
191 binmap->SetBinContent(ibinx, ibiny, bincounter);
196 cout <<
"Counter=" << bincounter <<
endl;
203 assert(covTo2DBin && covTo2DBin->GetNbinsX() > 0);
204 for(
int i = 1;
i <= covTo2DBin->GetNbinsX(); ++
i)
205 if (covTo2DBin->GetBinContent(
i) == globalBin)
213 assert(hist->GetNbinsX() == mat->GetNbinsX());
224 assert(hist->GetNbinsX() == mat.GetNrows());
244 ibinEdge = axis->GetBinLowEdge(
ibin);
246 ibinEdge = axis->GetBinUpEdge(
ibin);
248 throw runtime_error(
"EBinEdge has value " +
to_string(findBinBy) +
", which is not defined.");
249 if (
abs(ibinEdge - binEdge) < epsilon){
255 throw runtime_error(
"Could not find bin matching " +
to_string(binEdge));
260 void OverlayGenerators(map<
string, map<string, TH1* > > generators,
string var,
bool drawSpline =
false, TLegend **
leg = NULL,
double binLowEdge = -1.)
262 int smoothFactor = 5;
263 cout <<
"Starting bin checks for var: " << var <<
endl;
264 for (pair<
string, map<string, TH1* > >
generator : generators){
272 TH1 * histSlice = hist->ProjectionY((hist->GetName() +
string(
"_") +
to_string(generatorBin)).c_str(), generatorBin, generatorBin);
274 histSlice->Rebin(smoothFactor);
275 histSlice->Scale(1.0/
double(smoothFactor));
278 spline =
new TSpline3(histSlice);
283 spline =
new TSpline3(
generator.second[var]);
286 spline->SetLineColor(
generator.second[var]->GetLineColor());
287 spline->SetLineWidth(3);
288 spline->Draw(
"same");
295 TH2 * hist2D = (TH2*)
generator.second[var];
296 hist = hist2D->ProjectionY((hist2D->GetName() +
string(
"_") +
to_string(generatorBin)).c_str(), generatorBin, generatorBin);
300 hist->SetLineColor(
generator.second[var]->GetLineColor());
301 hist->SetLineWidth(3);
302 hist->Draw(
"hist same");
309 TLegend * oldLeg = (TLegend*) (*leg)->Clone();
311 (*leg)->SetFillStyle(0);
312 (*leg)->SetTextSize(0.03);
313 for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
314 TLegendEntry *
entry = (TLegendEntry*) obj;
315 (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
322 TH1 * genClone = (TH1*) gen->Clone();
323 TH1 * genRat = (TH1*) gen->Clone();
324 assert(gen->GetNbinsX() == data->GetNbinsX());
325 genClone->SetTitle(data->GetTitle());
326 genRat->Divide(data);
327 genRat->GetYaxis()->SetTitle(
"Gen / Data");
329 return {genClone, genRat};
332 void PlotGenRatios(map<
string, map<string, TH1* > > generators, TH1* dataFull, TH1* dataStat,
string var, TLegend **
leg = NULL,
double binLowEdge = -1.)
334 TRatioPlot * baseRatio = NULL;
335 vector<pair<TH1*, TH1*> > ratios;
336 for (pair<
string, map<string, TH1* > >
generator : generators){
339 TH2 * hist2D = (TH2*)
generator.second[var];
341 hist = hist2D->ProjectionY((hist->GetName() +
string(
"_") +
to_string(generatorBin)).c_str(), generatorBin, generatorBin);
342 hist->SetLineColor(hist2D->GetLineColor());
345 ratios.push_back(genRatio);
347 baseRatio =
new TRatioPlot(hist, dataFull);
351 (*leg)->AddEntry(hist,
generator.first.c_str());
354 double maxhist = (*std::max_element(ratios.begin(), ratios.end(), [](
auto i,
auto j){
355 return (
i.first->GetMaximum() <
j.first->GetMaximum());
356 })).first->GetMaximum();
357 double maxrat = (*std::max_element(ratios.begin(), ratios.end(), [](
auto i,
auto j){
return i.second->GetMaximum() <
j.second->GetMaximum();})).second->GetMaximum();
358 double minrat = (*std::min_element(ratios.begin(), ratios.end(), [](
auto i,
auto j){
return i.second->GetMinimum() <
j.second->GetMinimum();})).second->GetMinimum();
361 baseRatio->SetSplitFraction(0.4);
362 baseRatio->SetSeparationMargin(0.0);
363 baseRatio->GetLowerRefYaxis()->SetTitle(
"Gen / Data");
364 ((TH1*) baseRatio->GetUpperRefObject())->
SetLineWidth(0.0);
365 ((TH1*) baseRatio->GetUpperRefObject())->
SetTitle(
"Generator Ratios");
366 ((TH1*) baseRatio->GetUpperRefObject())->SetTitleSize(0.035,
"Y");
367 baseRatio->GetLowerRefGraph()->SetMarkerSize(0.0);
368 baseRatio->GetLowerRefGraph()->GetYaxis()->SetTitleSize(0.035);
369 baseRatio->SetGraphDrawOpt(
"X");
370 baseRatio->GetUpperRefYaxis()->SetRangeUser(0.0, maxhist * 1.2);
371 baseRatio->GetLowerRefYaxis()->SetRangeUser(minrat * 0.9, 1.1 * maxrat);
372 baseRatio->GetUpperRefXaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
373 baseRatio->GetLowYaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
374 baseRatio->GetUpperRefYaxis()->SetTitle(dataFull->GetYaxis()->GetTitle());
375 baseRatio->GetUpperRefYaxis()->CenterTitle(
false);
376 baseRatio->GetLowerRefYaxis()->CenterTitle(
true);
377 baseRatio->GetLowerRefYaxis()->SetTitleOffset(1.3);
378 baseRatio->SetGridlines({});
379 baseRatio->Draw(
"noconfint");
380 for(pair<TH1*, TH1*>
ratio : ratios)
382 TH1 * gen =
ratio.first;
383 TH1 * rat =
ratio.second;
384 baseRatio->GetUpperPad()->cd();
385 gen->Draw(
"hist same");
386 baseRatio->GetLowerPad()->cd();
387 rat->Draw(
"hist same");
389 baseRatio->GetUpperPad()->cd();
390 dataFull->Draw(
"e1 same");
391 dataStat->Draw(
"e1 same");
394 TLegend * oldLeg = (TLegend*) (*leg)->Clone();
396 (*leg)->SetFillStyle(0);
397 (*leg)->SetTextSize(0.03);
398 for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
399 TLegendEntry *
entry = (TLegendEntry*) obj;
400 (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
407 void PrintAllSlices(TH2* histSystPlusStat, TH2* histStat,
string histLabel,
string filename,
bool drawSpline, map<
string, map<string, TH1* > > generators = {},
bool printRatios =
false)
409 TCanvas *
canv =
new TCanvas(
"sliceCanv",
"sliceCanv", 1600, 1200);
410 canv->SetLeftMargin(0.14);
411 canv->SetRightMargin(0.05);
412 canv->SetBottomMargin(0.12);
414 for (
int ibinx = 2; ibinx <= histSystPlusStat->GetNbinsX(); ++ibinx){
415 TH1 * proj0 = histSystPlusStat->ProjectionY((histSystPlusStat->GetName() +
string(
"_") +
to_string(ibinx)).c_str(), ibinx, ibinx);
416 proj0->SetTitle(
TString::Format(
"%3.2f < cos(#theta_{#mu}) < %3.2f", histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx), histSystPlusStat->GetXaxis()->GetBinUpEdge(ibinx)));
418 proj0->GetXaxis()->SetRangeUser(0.5, 2.5);
419 proj0->GetYaxis()->SetRangeUser(0.0, proj0->GetMaximum() * 1.25);
420 proj0->SetMarkerStyle(histSystPlusStat->GetMarkerStyle());
423 leg->SetFillStyle(0);
424 leg->AddEntry(proj0, histLabel.c_str());
425 proj0->Draw(
"E1 SAME");
426 OverlayGenerators(generators,
"MuKin", drawSpline, &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
427 TH1 * projStat0 = histStat->ProjectionY((histStat->GetName() +
string(
"_") +
to_string(ibinx)).c_str(), ibinx, ibinx);
428 projStat0->SetLineColor(kGray);
429 projStat0->Draw(
"E1 SAME");
431 canv->SaveAs((filename +
"_bin" +
to_string(ibinx) +
".pdf").c_str());
436 leg->AddEntry(proj0, histLabel.c_str());
437 PlotGenRatios(generators, proj0, projStat0,
"MuKin", &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
439 canv->SaveAs((filename +
"_bin" +
to_string(ibinx) +
"_ratios.pdf").c_str());
451 for (
int ibinx = 1; ibinx <= hist->GetNbinsX(); ++ibinx){
452 hist->SetBinContent(ibinx, hist->GetBinContent(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
453 hist->SetBinError(ibinx, hist->GetBinError(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
461 {
"MuKin",
"T_muCos_theta"},
462 {
"ENu",
"xsecVsEnu_2"},
468 {
"MuKin",
"T_muCos_theta"},
469 {
"ENu",
"xsecVsEnu_2"},
470 {
"Q2",
"xsecVsq2VB_2"}
475 {
"MuKin",
"T_muCos_theta_Soup"},
476 {
"ENu",
"xsecVsEnu_2_Soup"},
477 {
"Q2",
"xsecVsq2_2_Soup"}
482 {
"MuKin",
"T_muCos_theta_Soup"},
483 {
"ENu",
"xsecVsEnu_2_Soup"},
484 {
"Q2",
"xsecVsq2VB_2_Soup"}
492 map<string, TH1*> generatorMap;
503 generatorMap[
"MuKin"]->SetLineColor(col);
504 generatorMap[
"ENu"] ->SetLineColor(col);
505 generatorMap[
"Q2"] ->SetLineColor(col);
507 string appendLabel = histLabel ==
"rebin" || histLabel ==
"genie3_rebin" ?
"rebin" :
"";
508 generatorMap[
"MuKin"]->SetName((label +
"_MuKin" + appendLabel).c_str());
509 generatorMap[
"ENu"] ->SetName((label +
"_ENu" + appendLabel).c_str());
510 generatorMap[
"Q2"] ->SetName((label +
"_Q2" + appendLabel).c_str());
514 generatorMap[
"MuKin"]->Scale(
scaling);
515 generatorMap[
"ENu"] ->Scale(
scaling);
516 generatorMap[
"Q2"] ->Scale(
scaling);
519 generatorMap[
"MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
520 generatorMap[
"MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
521 generatorMap[
"ENu"] ->GetXaxis()->SetRangeUser(0.5, 5.0);
522 generatorMap[
"Q2"] ->GetXaxis()->SetRangeUser(0.0, 2.8);
532 map<string, TH1*> geniePred;
533 geniePred[
"MuKin"] = (TH2*) file->Get((
"geniePrediction/" + weight +
"/MuKin").c_str())->
Clone();
534 geniePred[
"ENu"] = (TH1*) file->Get((
"geniePrediction/" + weight +
"/ENu").c_str())->
Clone();
535 geniePred[
"Q2"] = (TH1*) file->Get((
"geniePrediction/" + weight +
"/Q2").c_str())->
Clone();
536 geniePred[
"MuKin"]->SetLineColor(col);
537 geniePred[
"ENu"] ->SetLineColor(col);
538 geniePred[
"Q2"] ->SetLineColor(col);
539 geniePred[
"MuKin"]->SetName((label +
"_MuKin").c_str());
540 geniePred[
"ENu"] ->SetName((label +
"_ENu" ).c_str());
541 geniePred[
"Q2"] ->SetName((label +
"_Q2" ).c_str());
544 geniePred[
"MuKin"]->Scale(1E39);
545 geniePred[
"ENu"] ->Scale(1E39);
546 geniePred[
"Q2"] ->Scale(1E39);
548 geniePred[
"MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
549 geniePred[
"MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
550 geniePred[
"ENu"] ->GetXaxis()->SetRangeUser(0.5, 5.0);
551 geniePred[
"Q2"] ->GetXaxis()->SetRangeUser(0.0, 2.8);
561 double genBinLow = generator->FindBin(binning[0]);
562 double genBinUp = generator->FindBin(binning[binning.size() - 1]) - 1;
563 double dataBinLow = data->FindBin(binning[0]);
564 double dataBinUp = data->FindBin(binning[binning.size() - 1]) - 1;
565 cout <<
"Scaling " << label <<
": " << data->Integral(dataBinLow, dataBinUp,
"width") / generator->Integral(genBinLow, genBinUp,
"width") <<
endl;
566 generator->Scale(data->Integral(dataBinLow, dataBinUp,
"width") / generator->Integral(genBinLow, genBinUp,
"width"));
572 double genIntegral = 0.;
573 double dataIntegral = 0.;
574 for(
int ibinx = 1; ibinx <= data->GetNbinsX(); ++ibinx)
575 for(
int ibiny = 1; ibiny <= data->GetNbinsY(); ++ibiny)
576 if (data->GetBinContent(ibinx, ibiny) > 0.0){
579 dataIntegral += data->GetBinContent(ibinx, ibiny) * data->GetXaxis()->GetBinWidth(ibinx) * data->GetYaxis()->GetBinWidth(ibiny);
581 for(
int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx)
582 for(
int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
583 int dataBinX = data->GetXaxis()->FindBin(generator->GetXaxis()->GetBinCenter(ibinx));
584 int dataBinY = data->GetYaxis()->FindBin(generator->GetYaxis()->GetBinCenter(ibiny));
585 if (data->GetBinContent(dataBinX, dataBinY) > 0.0){
588 genIntegral += generator->GetBinContent(ibinx, ibiny) * generator->GetXaxis()->GetBinWidth(ibinx) * generator->GetYaxis()->GetBinWidth(ibiny);
591 cout <<
"Scaling: " << label <<
": " << dataIntegral / genIntegral <<
endl;
592 generator->Scale(dataIntegral / genIntegral);
599 for (
int ibinx = binxLow; ibinx <= binxHigh; ++ibinx)
600 widthx += hist->GetXaxis()->GetBinWidth(ibinx);
601 for (
int ibiny = binyLow; ibiny <= binyHigh; ++ibiny)
602 widthy += hist->GetYaxis()->GetBinWidth(ibiny);
603 return widthx * widthy;
606 map<string, map<string, TH1* > >
RebinGenerators(map<
string, map<string, TH1* > > generators)
608 map<string, map<string, TH1* > > rebinnedGenerators;
609 for (pair<
string, map<string, TH1* > > generatorPair : generators){
610 string genLabel = generatorPair.first;
611 map<string, TH1* > &
generator = generatorPair.second;
614 TH2D * mukin_hist =
new TH2D((generator[
"MuKin"]->
GetName() +
string(
"_databins")).c_str(), (genLabel +
" Rebin").c_str(),
616 TH1D * enu_hist =
new TH1D((generator[
"ENu"]->
GetName() +
string(
"_databins")).c_str(), (genLabel +
" Rebin").c_str(),
618 TH1D * q2_hist =
new TH1D((generator[
"Q2"]->
GetName() +
string(
"_databins")).c_str(), (genLabel +
" Rebin").c_str(),
621 mukin_hist->SetLineColor(generator[
"MuKin"]->GetLineColor());
622 enu_hist ->SetLineColor(generator[
"ENu"] ->GetLineColor());
623 q2_hist ->SetLineColor(generator[
"Q2"]->GetLineColor());
638 double binIntegral = ((TH2*) generator[
"MuKin"])->Integral(binxLow, binxHigh, binyLow, binyHigh);
639 int numBins = (binxHigh - binxLow + 1) * (binyHigh - binyLow + 1);
640 mukin_hist->SetBinContent(ibinx + 1, ibiny + 1, binIntegral / numBins);
645 for (
unsigned int ibinx = 0; ibinx <
enuBinEdges.size() - 1; ++ibinx)
651 double binIntegral = generator[
"ENu"]->Integral(binLow, binHigh);
652 enu_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
656 for (
unsigned int ibinx = 0; ibinx <
q2BinEdges.size() - 1; ++ibinx)
662 double binIntegral = generator[
"Q2"]->Integral(binLow, binHigh);
663 q2_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
666 rebinnedGenerators[generatorPair.first][
"MuKin"] = mukin_hist;
667 rebinnedGenerators[generatorPair.first][
"ENu"] = enu_hist;
668 rebinnedGenerators[generatorPair.first][
"Q2"] = q2_hist;
670 return rebinnedGenerators;
678 double naiveChi2 = 0;
679 assert(invCovMat.GetNrows() == invCovMat.GetNcols());
681 assert(invCovMat.GetNrows() == data->GetNbinsX());
682 assert(invCovMat.GetNrows() == generator->GetNbinsX());
683 TH2D * chi2Contribs =
new TH2D((label +
"_enu_chi2_contribs").c_str(), (label +
" E_{#nu} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
685 for (
int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
687 double binChi2 = (generator->GetBinContent(
ibin + 1) - data->GetBinContent(
ibin + 1)) * invCovMat[
ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
688 chi2Contribs->SetBinContent(
ibin + 1, jbin + 1, binChi2);
691 naiveChi2 +=
pow(generator->GetBinContent(
ibin + 1) - data->GetBinContent(
ibin + 1), 2) /
pow(data->GetBinError(
ibin + 1), 2);
693 chi2Contribs->Draw(
"COLZ");
694 cout <<
"\tNaive Chi2: " << naiveChi2 <<
endl;
702 for (
int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
703 chi2 += (generator->GetBinContent(
ibin + 1) - data->GetBinContent(
ibin + 1)) * invCovMat[
ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
710 void fcn_enu(
int &npar,
double *gin,
double &
f,
double *
par,
int iflag)
713 TH1 * scaled_gen = (TH1*) chi2_min_enu_gen->Clone(
"temp");
714 scaled_gen->Scale(par[0]);
717 scaled_gen->Delete();
723 double naiveChi2 = 0;
724 assert(invCovMat.GetNrows() == invCovMat.GetNcols());
726 assert(invCovMat.GetNrows() == data->GetNbinsX());
727 assert(invCovMat.GetNrows() == generator->GetNbinsX());
728 TH2D * chi2Contribs =
new TH2D((label +
"_q2_chi2_contribs").c_str(), (label +
" q^{2} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
730 for (
int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
732 double binChi2 = (generator->GetBinContent(
ibin + 1) - data->GetBinContent(
ibin + 1)) * invCovMat[
ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
733 chi2Contribs->SetBinContent(
ibin + 1, jbin + 1, binChi2);
736 naiveChi2 +=
pow(generator->GetBinContent(
ibin + 1) - data->GetBinContent(
ibin + 1), 2) /
pow(data->GetBinError(
ibin + 1), 2);
738 chi2Contribs->Draw(
"COLZ");
739 cout <<
"\tNaive Chi2: " << naiveChi2 <<
endl;
746 void fcn_q2(
int &npar,
double *gin,
double &
f,
double *
par,
int iflag)
749 TH1 * scaled_gen = (TH1*) chi2_min_q2_gen->Clone(
"temp");
750 scaled_gen->Scale(par[0]);
753 scaled_gen->Delete();
759 double naiveChi2 = 0;
760 assert(invCovMat.GetNrows() == invCovMat.GetNcols());
761 assert(invCovMat.GetNrows() == covToGlobalBinMap->GetNbinsX());
762 assert(data->GetNbinsX() == generator->GetNbinsX());
763 assert(data->GetNbinsY() == generator->GetNbinsY());
764 TH2D * chi2Contribs =
new TH2D((label +
"_mukin_chi2_contribs").c_str(), (label +
"Muon Kinematic #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
766 for (
int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx){
767 for (
int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
768 double xval = generator->GetXaxis()->GetBinLowEdge(ibinx);
769 double yval = generator->GetYaxis()->GetBinLowEdge(ibiny);
772 cerr <<
"\tDid not find bin: (" << generator->GetXaxis()->GetBinLowEdge(ibinx) <<
", " << generator->GetYaxis()->GetBinLowEdge(ibiny) <<
") in list of bins." <<
endl;
773 cerr <<
"\tBinning (2D) does not match! Skipping generator." <<
endl;
779 for (
int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
781 int iGlobalBin = covToGlobalBinMap->GetBinContent(
ibin + 1);
782 int jGlobalBin = covToGlobalBinMap->GetBinContent(jbin + 1);
788 data->GetBinXYZ(iGlobalBin, ibinX, ibinY, temp);
789 data->GetBinXYZ(jGlobalBin, jbinX, jbinY, temp);
790 double binChi2 = (generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY)) * invCovMat[
ibin][jbin] * (generator->GetBinContent(jbinX, jbinY) - data->GetBinContent(jbinX, jbinY));
791 chi2Contribs->SetBinContent(
ibin + 1, jbin + 1, binChi2);
793 naiveChi2 +=
pow(generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY), 2) /
pow(data->GetBinError(ibinX, ibinY), 2);
803 chi2Contribs->Draw(
"COLZ");
804 cout <<
"\tNaive Chi2: " << naiveChi2 <<
endl;
808 void PlotDataGenerators(TH1 * dataSystPlusStat, TH1 * dataStat,
string dataTitle, map<
string, map<string, TH1* > > &generatorMap,
string varName,
bool drawSpline =
false)
810 dataSystPlusStat->Draw(
"E1");
811 TLegend *
leg =
new TLegend(0.6, 0.6, 0.9, 0.85);
812 leg->SetFillStyle(0);
813 leg->AddEntry(dataSystPlusStat, dataTitle.c_str());
815 dataSystPlusStat->Draw(
"E1 SAME");
816 dataStat->SetLineColor(kGray);
817 dataStat->Draw(
"E1 SAME");
823 int nrowcols = hist->GetNbinsX();
824 TMatrixD statErrs(nrowcols, nrowcols);
832 TCanvas * tempCanv =
new TCanvas(
"temp",
"temp", 1600, 1200);
833 tempCanv->SetRightMargin(0.15);
835 TH2D * covHist =
new TH2D(cov);
836 covHist->SetName((label +
"CovHist").c_str());
837 covHist->SetTitle((title +
" Covariance Matrix;;;#times 10^{-78}").c_str());
838 covHist->Draw(
"COLZ");
839 tempCanv->SaveAs((outputFolder +
"/" + label +
"CovMat.pdf").c_str());
841 TH2D * invHist =
new TH2D(inv);
842 invHist->SetName((label +
"InvCovHist").c_str());
843 invHist->SetTitle((title +
" Inv. Cov Matrix;;;#times 10^{78}").c_str());
844 invHist->Draw(
"COLZ");
845 tempCanv->SaveAs((outputFolder +
"/" + label +
"InvCovMat.pdf").c_str());
848 TH2D * identityHist =
new TH2D(identity);
849 identityHist->SetName((label +
"Identity").c_str());
850 identityHist->SetTitle((title +
" Cov #times InvCov").c_str());
851 identityHist->Draw(
"COLZ");
852 tempCanv->SaveAs((outputFolder +
"/" + label +
"Identity.pdf").c_str());
863 gStyle->SetTitleYSize(0.05);
866 TCanvas *
c1 =
new TCanvas(
"c1",
"c1", 1600, 1200);
867 string shapeOnlyFilename =
"CovCurrent/Calculate_covariances_shape_all.root";
868 string totalErrFilename =
"CovCurrent/Calculate_covariances_total_all.root";
870 TFile *
outFile = TFile::Open((outputFolder +
"/numuCCXSecResults.root").c_str(),
"RECREATE");
871 TDirectory * xsecDir = outFile->mkdir(
"results");
872 TDirectory * generatorDir = outFile->mkdir(
"generators");
880 TFile * Stat2DFile = TFile::Open(
"/nova/ana/users/connorj/NumuCC/NERSC_Run_Dec2019/xsecUncertResults.root");
881 TFile * Stat1DFile = TFile::Open(
"/nova/ana/users/connorj/NumuCC/Feb2020_MuonE_Angle_SystFix/xsecUncertResults_1DOnly.root");
882 assert(Stat2DFile && Stat2DFile->IsOpen());
883 assert(Stat1DFile && Stat1DFile->IsOpen());
884 TH2 * statMuKin_originalbins = (TH2*) Stat2DFile->Get(
"cv/xsec_cv_2DProj");
885 TH1 * statENu_originalbins = (TH1*) Stat1DFile->Get(
"cv/xsec_cv_ENu");
886 TH1 * statQ2_originalbins = (TH1*) Stat1DFile->Get(
"cv/xsec_cv_Q2");
887 statMuKin_originalbins->Scale(1E39);
888 statENu_originalbins->Scale(1E39);
889 statQ2_originalbins->Scale(1E39);
909 statMuKin->SetBinContent(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinContent(binxLow, binyLow));
910 statMuKin->SetBinError(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinError(binxLow, binyLow));
915 TDirectory * statErrDir = xsecDir->mkdir(
"statErrs");
921 TMatrixD statMatMuKin(bin2DToCov->GetMaximum(), bin2DToCov->GetMaximum());
922 for(
int ibin = 0;
ibin < statMatMuKin.GetNrows(); ++
ibin){
923 unsigned int binx = 0;
924 unsigned int biny = 0;
925 for(
int ibinx = 1; ibinx <= bin2DToCov->GetNbinsX(); ++ibinx)
926 for(
int ibiny = 1; ibiny <= bin2DToCov->GetNbinsY(); ++ibiny)
927 if (
ibin == bin2DToCov->GetBinContent(ibinx, ibiny))
934 statMatMuKin[
ibin][
ibin] = statMuKin->GetBinError(binx, biny);
944 TFile * totalErrFile = TFile::Open(totalErrFilename.c_str());
945 TH1 * nomMuKinTotal1D = (TH1*) totalErrFile->Get(
"hcv_mukin");
946 TH1 * nomENuTotal_nobininfo = (TH1*) totalErrFile->Get(
"hcv_enu");
947 TH1 * nomQ2Total_nobininfo = (TH1*) totalErrFile->Get(
"hcv_q2");
948 nomMuKinTotal1D ->Scale(1E39);
949 nomENuTotal_nobininfo->Scale(1E39);
950 nomQ2Total_nobininfo ->Scale(1E39);
953 TH2 * covTotalMuKin = (TH2*) totalErrFile->Get(
"cov_mukin");
954 TH2 * covTotalENu = (TH2*) totalErrFile->Get(
"cov_enu");
955 TH2 * covTotalQ2 = (TH2*) totalErrFile->Get(
"cov_q2");
956 covTotalMuKin->Scale(1E78);
957 covTotalENu ->Scale(1E78);
958 covTotalQ2 ->Scale(1E78);
961 TMatrixD covMatTotalMuKin(covTotalMuKin->GetNbinsX(), covTotalMuKin->GetNbinsY());
962 TMatrixD covMatTotalENu (covTotalENu ->GetNbinsX(), covTotalENu ->GetNbinsY());
963 TMatrixD covMatTotalQ2 (covTotalQ2 ->GetNbinsX(), covTotalQ2 ->GetNbinsY());
964 for (
int ibinx = 1; ibinx <= covTotalMuKin->GetNbinsX(); ++ibinx)
965 for (
int ibiny = 1; ibiny <= covTotalMuKin->GetNbinsY(); ++ibiny)
966 covMatTotalMuKin[ibinx - 1][ibiny - 1] = covTotalMuKin->GetBinContent(ibinx, ibiny);
967 for (
int ibinx = 1; ibinx <= covTotalENu->GetNbinsX(); ++ibinx)
968 for (
int ibiny = 1; ibiny <= covTotalENu->GetNbinsY(); ++ibiny)
969 covMatTotalENu[ibinx - 1][ibiny - 1] = covTotalENu->GetBinContent(ibinx, ibiny);
970 for (
int ibinx = 1; ibinx <= covTotalQ2->GetNbinsX(); ++ibinx)
971 for (
int ibiny = 1; ibiny <= covTotalQ2->GetNbinsY(); ++ibiny)
972 covMatTotalQ2[ibinx - 1][ibiny - 1] = covTotalQ2->GetBinContent(ibinx, ibiny);
973 TMatrixD covMatTotalMuKinFull = covMatTotalMuKin + statMatMuKin;
974 TMatrixD covMatTotalENuFull = covMatTotalENu + statMatENu;
975 TMatrixD covMatTotalQ2Full = covMatTotalQ2 + statMatQ2;
976 TMatrixD invCovMatTotalMuKin = covMatTotalMuKinFull;
977 TMatrixD invCovMatTotalENu = covMatTotalENuFull;
978 TMatrixD invCovMatTotalQ2 = covMatTotalQ2Full;
979 invCovMatTotalMuKin.Invert();
980 invCovMatTotalENu .Invert();
981 invCovMatTotalQ2 .Invert();
989 TH2 * nomMuKinTotal =
Construct2D(nomMuKinTotal1D);
990 TH1 * nomENuTotal =
ConstructENu(nomENuTotal_nobininfo);
991 TH1 * nomQ2Total =
ConstructQ2(nomQ2Total_nobininfo);
995 TH1 * nomMuKinTotalStats1D = (TH1*) nomMuKinTotal1D ->
Clone(
"hmukin_cv_2D_stats");
996 TH1 * nomENuTotalStats_nobininfo = (TH1*) nomENuTotal_nobininfo->Clone(
"henu_cv_enuaxis_stats");
997 TH1 * nomQ2TotalStats_nobininfo = (TH1*) nomQ2Total_nobininfo ->
Clone(
"hq2_cv_q2axis_stats");
1001 TH2 * nomMuKinStats =
Construct2D(nomMuKinTotalStats1D);
1002 TH1 * nomENuStats =
ConstructENu(nomENuTotalStats_nobininfo);
1003 TH1 * nomQ2Stats =
ConstructQ2(nomQ2TotalStats_nobininfo);
1010 nomMuKinTotal->SetMarkerStyle(20);
1015 TDirectory * totalUncertResults = xsecDir->mkdir(
"totalUncertResults");
1016 totalUncertResults->cd();
1017 nomMuKinTotal->Write(
"MuKin");
1018 nomENuTotal->Write(
"ENu");
1019 nomQ2Total->Write(
"Q2");
1026 TFile * shapeOnlyFile = TFile::Open(shapeOnlyFilename.c_str());
1027 TH1 * nomMuKinShape1D = (TH1*) shapeOnlyFile->Get(
"hcv_mukin");
1028 TH1 * nomENuShape_nobininfo = (TH1*) shapeOnlyFile->Get(
"hcv_enu");
1029 TH1 * nomQ2Shape_nobininfo = (TH1*) shapeOnlyFile->Get(
"hcv_q2");
1030 nomMuKinShape1D ->Scale(1E39);
1031 nomENuShape_nobininfo->Scale(1E39);
1032 nomQ2Shape_nobininfo ->Scale(1E39);
1035 TH2 * covShapeMuKin = (TH2*) shapeOnlyFile->Get(
"cov_mukin");
1036 TH2 * covShapeENu = (TH2*) shapeOnlyFile->Get(
"cov_enu");
1037 TH2 * covShapeQ2 = (TH2*) shapeOnlyFile->Get(
"cov_q2");
1038 covShapeMuKin->Scale(1E78);
1039 covShapeENu ->Scale(1E78);
1040 covShapeQ2 ->Scale(1E78);
1043 TMatrixD covMatShapeMuKin(covShapeMuKin->GetNbinsX(), covShapeMuKin->GetNbinsY());
1044 TMatrixD covMatShapeENu (covShapeENu ->GetNbinsX(), covShapeENu ->GetNbinsY());
1045 TMatrixD covMatShapeQ2 (covShapeQ2 ->GetNbinsX(), covShapeQ2 ->GetNbinsY());
1046 for (
int ibinx = 1; ibinx <= covShapeMuKin->GetNbinsX(); ++ibinx)
1047 for (
int ibiny = 1; ibiny <= covShapeMuKin->GetNbinsY(); ++ibiny)
1048 covMatShapeMuKin[ibinx - 1][ibiny - 1] = covShapeMuKin->GetBinContent(ibinx, ibiny);
1049 for (
int ibinx = 1; ibinx <= covShapeENu->GetNbinsX(); ++ibinx)
1050 for (
int ibiny = 1; ibiny <= covShapeENu->GetNbinsY(); ++ibiny)
1051 covMatShapeENu[ibinx - 1][ibiny - 1] = covShapeENu->GetBinContent(ibinx, ibiny);
1052 for (
int ibinx = 1; ibinx <= covShapeQ2->GetNbinsX(); ++ibinx)
1053 for (
int ibiny = 1; ibiny <= covShapeQ2->GetNbinsY(); ++ibiny)
1054 covMatShapeQ2[ibinx - 1][ibiny - 1] = covShapeQ2->GetBinContent(ibinx, ibiny);
1055 TMatrixD covMatShapeMuKinFull = covMatShapeMuKin + statMatMuKin;
1056 TMatrixD covMatShapeENuFull = covMatShapeENu + statMatENu;
1057 TMatrixD covMatShapeQ2Full = covMatShapeQ2 + statMatQ2;
1058 TMatrixD invCovMatShapeMuKin = covMatShapeMuKinFull;
1059 TMatrixD invCovMatShapeENu = covMatShapeENuFull;
1060 TMatrixD invCovMatShapeQ2 = covMatShapeQ2Full;
1061 invCovMatShapeMuKin.Invert();
1062 invCovMatShapeENu .Invert();
1063 invCovMatShapeQ2 .Invert();
1067 covarianceDir->cd();
1068 PlotCovMatrices(covMatTotalENuFull, invCovMatTotalENu, outputFolder,
"enuTotal",
"E_{#nu} Total");
1069 PlotCovMatrices(covMatTotalQ2Full, invCovMatTotalQ2, outputFolder,
"q2Total",
"q^{2} Total");
1070 PlotCovMatrices(covMatTotalMuKinFull, invCovMatTotalMuKin, outputFolder,
"mukinTotal",
"Muon Kinematic Total");
1071 PlotCovMatrices(covMatShapeENuFull, invCovMatShapeENu, outputFolder,
"enuShape",
"E_{#nu} Shape-only");
1072 PlotCovMatrices(covMatShapeQ2Full, invCovMatShapeQ2, outputFolder,
"q2Shape",
"q^{2} Shape-only");
1073 PlotCovMatrices(covMatShapeMuKinFull, invCovMatShapeMuKin, outputFolder,
"mukinShape",
"Muon Kinematic Shape-only");
1081 TH2 * nomMuKinShape =
Construct2D(nomMuKinShape1D);
1082 TH1 * nomENuShape =
ConstructENu(nomENuShape_nobininfo);
1083 TH1 * nomQ2Shape =
ConstructQ2(nomQ2Shape_nobininfo);
1102 nomMuKinShape->SetMarkerStyle(20);
1107 TDirectory * shapeUncertResults = xsecDir->mkdir(
"shapeUncertResults");
1108 shapeUncertResults->cd();
1109 nomMuKinShape->Write(
"MuKin");
1110 nomENuShape->Write(
"ENu");
1111 nomQ2Shape->Write(
"Q2");
1116 TH1I * covTo2DBin =
MakeCovTo2DBin(nomMuKinTotal, covMatTotalMuKin.GetNrows());
1118 covTo2DBin->Write(
"CovTo2DGlobalBin");
1119 bin2DToCov->Write(
"Bin2DToCov");
1125 cout <<
"Started loading generators" <<
endl;
1128 TFile * gibuuFile = TFile::Open(
"GiBUU_numuCC_predictions.root");
1129 TFile * neutFile = TFile::Open(
"NEUT54_numuCC_predictions.root");
1130 TFile * nuwroFile = TFile::Open(
"NuWro_numuCC_predictions.root");
1131 TFile * genieFile = TFile::Open(
"GeniePrediction.root");
1132 TFile * genie3File = TFile::Open(
"Genie3Predictions.root");
1133 assert(gibuuFile && gibuuFile->IsOpen());
1134 assert(neutFile && neutFile->IsOpen());
1135 assert(nuwroFile && nuwroFile->IsOpen());
1136 assert(genieFile && genieFile->IsOpen());
1137 assert(genie3File && genie3File->IsOpen());
1138 map<string, map<string, TH1* > > generators;
1143 generators[
"GENIE 2.12.2 - Untuned"] =
LoadGeniePrediction(genieFile,
"novatune", kCyan,
"ppfx");
1147 map<string, map<string, TH1* > > generatorsAreaNorm;
1148 for(pair<
string, map<string, TH1* > >
generator : generators)
1150 generatorsAreaNorm[
generator.first +
" (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() +
string(
"_norm")).c_str());
1151 for (pair<
string, map<string, TH1* > >
generator : generatorsAreaNorm)
1159 map<string, map<string, TH1* > > generatorsVarBin;
1160 generatorsVarBin[
"GiBUU 2019"] =
LoadGenerator(gibuuFile,
"gibuu",
kBlue, 1E39,
"rebin");
1164 generatorsVarBin[
"GENIE 2.12.2 - Untuned"] =
LoadGeniePrediction(genieFile,
"novatune_rebin", kCyan,
"ppfx");
1167 generatorsVarBin[
"GiBUU 2019"][
"ENu"]->Rebin(5);
1168 generatorsVarBin[
"GiBUU 2019"][
"ENu"]->Scale(1./5);
1169 generatorsVarBin[
"GiBUU 2019"][
"ENu"]->GetXaxis()->SetRangeUser(0.0, 5.0);
1170 generatorsVarBin[
"NuWro 2019"][
"ENu"]->Rebin(5);
1171 generatorsVarBin[
"NuWro 2019"][
"ENu"]->Scale(1./5);
1174 generatorsVarBin[
"NEUT 5.4.0"][
"ENu"]->Rebin(10);
1175 generatorsVarBin[
"NEUT 5.4.0"][
"ENu"]->Scale(1./10);
1178 map<string, map<string, TH1* > > generatorsAreaNormVarBin;
1179 for(pair<
string, map<string, TH1* > >
generator : generatorsVarBin)
1180 for(pair<string, TH1*> var :
generator.second)
1181 generatorsAreaNormVarBin[
generator.first +
" (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() +
string(
"_norm")).c_str());
1182 for (pair<
string, map<string, TH1* > >
generator : generatorsAreaNormVarBin)
1188 c1->SetRightMargin(0.15);
1192 map<string, map<string, TH1* > > generatorsRebin =
RebinGenerators(generatorsVarBin);
1193 map<string, map<string, TH1* > > generatorsAreaNormRebin =
RebinGenerators(generatorsAreaNormVarBin);
1194 generatorDir->mkdir(
"standard")->cd();
1197 var.second->Write((
generator.first +
"_" + var.first).c_str());
1198 generatorDir->mkdir(
"normalized")->cd();
1199 for (
auto generator : generatorsAreaNorm)
1201 var.second->Write((
generator.first +
"_" + var.first +
"_norm").c_str());
1202 generatorDir->mkdir(
"databinned")->cd();
1205 var.second->Write((
generator.first +
"_" + var.first +
"_databins").c_str());
1207 generatorDir->mkdir(
"databinned_normalized")->cd();
1208 for (
auto generator : generatorsAreaNormRebin)
1210 var.second->Write((
generator.first +
"_" + var.first +
"_databins_norm").c_str());
1218 map<string, map<string, double > > generatorChi2s;
1219 map<string, map<string, TH1* > > chi2_min_scaled_generators;
1220 map<string, map<string, double > > generatorScaledChi2s;
1221 map<string, string> gen_short_forms
1223 {
"GiBUU 2019",
"GiBUU"},
1224 {
"NEUT 5.4.0",
"NEUT"},
1225 {
"NuWro 2019",
"NuWro"},
1226 {
"GENIE 2.12.2 - Untuned",
"Genie_Untuned"},
1227 {
"GENIE 2.12.2 - NOvA Tune",
"Genie_Tuned"},
1228 {
"GENIE 3.00.06",
"Genie_3"}
1304 string dataTotalTitle =
"Data - Total error";
1305 string dataShapeTitle =
"Data - Shape only errs";
1311 c1->SetRightMargin(0.15);
1312 nomMuKinTotal->Draw(
"COLZ");
1313 c1->SaveAs((outputFolder +
"/xsec_2dProj.pdf").c_str());
1316 c1->SetLeftMargin(0.13);
1317 c1->SetRightMargin(0.05);
1318 c1->SetBottomMargin(0.10);
1321 PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder +
"/muKinSlices",
true, generators);
1322 PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder +
"/muKinSlices_norm",
true, generatorsAreaNorm);
1323 PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder +
"/muKinSlices_rebin",
false, generatorsRebin,
true);
1324 PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder +
"/muKinSlices_norm_rebin",
false, generatorsAreaNormRebin,
true);
1327 nomENuTotal->GetYaxis()->SetRangeUser(0.0, nomENuTotal->GetMaximum() * 1.35);
1328 PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generators,
"ENu",
true);
1329 c1->SaveAs((outputFolder +
"/results_enu.pdf").c_str());
1330 PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNorm,
"ENu",
true);
1331 c1->SaveAs((outputFolder +
"/results_enu_norm.pdf").c_str());
1332 PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generatorsRebin,
"ENu",
false);
1333 c1->SaveAs((outputFolder +
"/results_enu_rebin.pdf").c_str());
1334 PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNormRebin,
"ENu",
false);
1335 c1->SaveAs((outputFolder +
"/results_enu_norm_rebin.pdf").c_str());
1338 nomQ2Total->GetYaxis()->SetRangeUser(0.0, nomQ2Total->GetMaximum() * 1.25);
1339 nomQ2Shape->GetYaxis()->SetRangeUser(0.0, nomQ2Shape->GetMaximum() * 1.25);
1341 c1->SaveAs((outputFolder +
"/results_q2.pdf").c_str());
1342 PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNorm,
"Q2",
true);
1343 c1->SaveAs((outputFolder +
"/results_q2_norm.pdf").c_str());
1344 PlotDataGenerators(nomQ2Total, nomQ2Stats, dataTotalTitle, generatorsRebin,
"Q2",
false);
1345 c1->SaveAs((outputFolder +
"/results_q2_rebin.pdf").c_str());
1346 PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNormRebin,
"Q2",
false);
1347 c1->SaveAs((outputFolder +
"/results_q2_norm_rebin.pdf").c_str());
1352 c1 =
new TCanvas(
"c1",
"c1", 1600, 1320);
1353 leg =
new TLegend();
1354 leg->AddEntry(nomENuTotal, dataTotalTitle.c_str());
1355 PlotGenRatios(generatorsRebin, nomENuTotal, nomENuStats,
"ENu", &leg);
1356 c1->SaveAs((outputFolder +
"/results_enu_ratios.pdf").c_str());
1361 leg =
new TLegend();
1362 leg->AddEntry(nomENuShape, dataShapeTitle.c_str());
1363 PlotGenRatios(generatorsAreaNormRebin, nomENuShape, nomENuStats,
"ENu", &leg);
1364 c1->SaveAs((outputFolder +
"/results_enu_shape_ratios.pdf").c_str());
1369 leg =
new TLegend();
1370 leg->AddEntry(nomQ2Total, dataTotalTitle.c_str());
1371 PlotGenRatios(generatorsRebin, nomQ2Total, nomQ2Stats,
"Q2", &leg);
1372 c1->SaveAs((outputFolder +
"/results_q2_ratios.pdf").c_str());
1377 leg =
new TLegend();
1378 leg->AddEntry(nomQ2Shape, dataShapeTitle.c_str());
1379 PlotGenRatios(generatorsAreaNormRebin, nomQ2Shape, nomQ2Stats,
"Q2", &leg);
1380 c1->SaveAs((outputFolder +
"/results_q2_shape_ratios.pdf").c_str());
void fcn_enu(int &npar, double *gin, double &f, double *par, int iflag)
pair< TH1 *, TH1 * > GeneratorRatio(TH1 *gen, TH1 *data)
TH1 * SetErrorFromCovDiag(TH1 *hist, TH2 *mat)
TMatrixD * chi2_min_enu_invCovMat
map< int, int > MuTperCosBins
TMatrixD * chi2_min_q2_invCovMat
correl_xv GetYaxis() -> SetDecimals()
map< string, TH1 * > LoadGenerator(TFile *file, string label, Color_t col, double scaling=-1, string histLabel="default")
const std::vector< double > angBinEdges_reduced
TH1 * ConstructQ2(TH1 *histIn)
void OverlayGenerators(map< string, map< string, TH1 * > > generators, string var, bool drawSpline=false, TLegend **leg=NULL, double binLowEdge=-1.)
int Get2DBinArea(TH2 *hist, int binxLow, int binxHigh, int binyLow, int binyHigh)
void PlotDataGenerators(TH1 *dataSystPlusStat, TH1 *dataStat, string dataTitle, map< string, map< string, TH1 * > > &generatorMap, string varName, bool drawSpline=false)
map< string, map< string, string > > generator_hist_names
const std::vector< double > mukeBinEdges_reduced
int GetBinByEdge(TAxis *axis, double binEdge, EBinEdge findBinBy=kBinLowEdge)
const std::vector< double > enuBinEdges
map< string, TH1 * > LoadGeniePrediction(TFile *file, string label, Color_t col, string weight)
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
fVtxDx SetMarkerStyle(20)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
int GetCovBin(TH1I *covTo2DBin, int globalBin)
void SigmaDivideENu(TH1 *hist)
void scaling(TH1D *hIn, const double shape_scale)
double CalculateENuChi2(const TMatrixD &invCovMat, TH1 *data, TH1 *generator, string label)
const XML_Char const XML_Char * data
std::string GetName(int i)
void PrintAllSlices(TH2 *histSystPlusStat, TH2 *histStat, string histLabel, string filename, bool drawSpline, map< string, map< string, TH1 * > > generators={}, bool printRatios=false)
TH2 * Construct2D(TH1 *histIn)
TMatrixD StatMatFromHist(TH1 *hist)
correl_xv GetXaxis() -> SetDecimals()
const std::vector< double > q2BinEdges
const std::vector< double > mukeBinEdges
void fcn_q2(int &npar, double *gin, double &f, double *par, int iflag)
const std::vector< double > angBinEdges
TH1 * ConstructENu(TH1 *histIn)
void AreaNormalize(TH1 *generator, TH1 *data, const vector< double > &binning, const string label)
const std::string covarianceDir
void PlotCovariance(string outputFolder=".")
double Calculate1DChi2_simple(const TMatrixD &invCovMat, const TH1 *data, const TH1 *generator)
map< string, map< string, TH1 * > > RebinGenerators(map< string, map< string, TH1 * > > generators)
string cv_name
A script which takes covariance matrices, and constructs our final pretty plots for NumuCC...
double Calculate2DChi2(const TMatrixD &invCovMat, TH2 *data, TH2 *generator, TH1 *covToGlobalBinMap, string label)
TH1I * MakeCovTo2DBin(TH2 *hist, int nbins)
double CalculateQ2Chi2(const TMatrixD &invCovMat, TH1 *data, TH1 *generator, string label)
assert(nhit_max >=nhit_nbins)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
std::string to_string(ModuleType mt)
void PlotCovMatrices(TMatrixD &cov, TMatrixD &inv, string outputFolder, string label, string title)
void PlotGenRatios(map< string, map< string, TH1 * > > generators, TH1 *dataFull, TH1 *dataStat, string var, TLegend **leg=NULL, double binLowEdge=-1.)
void AreaNormalize2D(TH2 *generator, TH2 *data, const string label)
fvar< T > inv(const fvar< T > &x)