4#pragma GCC diagnostic ignored "-Wfloat-conversion"
5#pragma GCC diagnostic ignored "-Wuseless-cast"
13 MACH3LOG_DEBUG(
"With OpenMP and {} threads", omp_get_max_threads());
19 else MACH3LOG_INFO(
"Using alternative method of statistical fluctuation, which is much slower");
44 rnd = std::make_unique<TRandom3>();
77 lnLHist = std::make_unique<TH1D>(
"lnLHist_predpredfluc",
"lnLHist_predpredfluc", 100, 1, -1);
79 lnLHist->GetXaxis()->SetTitle(
"-2LLH (Pred Fluc, Pred)");
80 lnLHist->GetYaxis()->SetTitle(
"Counts");
82 lnLHist_drawdata = std::make_unique<TH1D>(
"lnLHist_drawdata",
"lnLHist_drawdata", 100, 1, -1);
87 lnLHist_drawfluc = std::make_unique<TH1D>(
"lnLHist_drawpredfluc",
"lnLHist_drawpredfluc", 100, 1, -1);
92 lnLHist_drawflucdraw = std::make_unique<TH1D>(
"lnLHist_drawflucdraw",
"lnLHist_drawflucdraw", 100, 1, -1);
97 lnLDrawHist = std::make_unique<TH2D>(
"lnLDrawHist",
"lnLDrawHist", 50, 1, -1, 50, 1, -1);
99 lnLDrawHist->GetXaxis()->SetTitle(
"-2LLH_{Pred Fluc, Draw}");
100 lnLDrawHist->GetYaxis()->SetTitle(
"-2LLH_{Data, Draw}");
102 lnLFlucHist = std::make_unique<TH2D>(
"lnLFlucHist",
"lnLFlucHist", 50, 1, -1, 50, 1, -1);
104 lnLFlucHist->GetXaxis()->SetTitle(
"-2LLH_{Draw Fluc, Draw}");
105 lnLFlucHist->GetYaxis()->SetTitle(
"-2LLH_{Data, Draw}");
107 lnLDrawHistRate = std::make_unique<TH2D>(
"lnLDrawHistRate",
"lnLDrawHistRate", 50, 1, -1, 50, 1, -1);
113 lnLFlucHist_ProjectX = std::make_unique<TH2D>(
"lnLFlucHist_ProjectX",
"lnLFlucHist_ProjectX", 50, 1, -1, 50, 1, -1);
125 std::stringstream ss;
126 ss <<
"Draws/" << binwidth;
127 RandomHist->GetYaxis()->SetTitle(ss.str().c_str());
175 if(
DataHist[i] ==
nullptr)
continue;
178 for (
int k = 1; k <=
maxBins[i]; ++k)
192 if(
DataHist[i] ==
nullptr)
continue;
221 MACH3LOG_ERROR(
"Size of SampleVector input != number of defined samples");
224 MACH3LOG_ERROR(
"Something has gone wrong with making the Samples");
235 const int Length = int(Data.size());
238 for (
int i = 0; i < Length; ++i) {
239 if (Data[i] ==
nullptr) {
245 std::string classname = std::string(
DataHist[i]->Class_Name());
246 if(classname ==
"TH2Poly")
248 DataHist[i] =
static_cast<TH2Poly*
>(Data[i]->Clone());
255 MACH3LOG_ERROR(
"Right now I only support TH2Poly but I am ambitious piece of code and surely will have more support in the future");
266 const int Length = int(Nominal.size());
270 for (
int i = 0; i < Length; ++i)
272 if (Nominal[i] ==
nullptr) {
288 NominalHist[i] =
static_cast<TH2Poly*
>(Nominal[i]->Clone());
290 W2NomHist[i] =
static_cast<TH2Poly*
>(NomW2[i]->Clone());
307 for (
int i = 0; i < Length; ++i) {
309 if (Nominal[i] !=
nullptr)
311 std::string name = std::string(
NominalHist[i]->GetName());
312 name = name.substr(0, name.find(
"_nom"));
319 for (
int j = 0; j <=
maxBins[i]; ++j)
323 lnLHist_Mean[i]->SetNameTitle((name+
"_MeanlnL").c_str(), (name+
"_MeanlnL").c_str());
325 lnLHist_Mean[i]->GetZaxis()->SetTitle(
"-2lnL_{sample} #times sign(MC-data)");
327 lnLHist_Mode[i]->SetNameTitle((name+
"_ModelnL").c_str(), (name+
"_ModelnL").c_str());
329 lnLHist_Mode[i]->GetZaxis()->SetTitle(
"-2lnL_{sample} #times sign(MC-data)");
331 lnLHist_Mean_ProjectX[i]->SetNameTitle((name+
"_MeanlnL_ProjectX").c_str(), (name+
"_MeanlnL_ProjectX").c_str());
335 MeanHist[i]->SetNameTitle((name+
"_mean").c_str(), (name+
"_mean").c_str());
337 MeanHist[i]->GetZaxis()->SetTitle(
"Mean");
341 MeanHistCorrected[i]->SetNameTitle((name+
"_mean_corrected").c_str(), (name+
"_mean_corrected").c_str());
345 std::vector<double> xbins;
346 std::vector<double> ybins;
352 ViolinHists_ProjectX[i] =
new TH2D((name+
"_Violin_ProjectX").c_str(), (name+
"_Violin_ProjectX").c_str(),
int(xbins.size()-1), &xbins[0] , 400, 0, MaxBinning);
357 ViolinHists_ProjectY[i] =
new TH2D((name+
"_Violin_ProjectY").c_str(), (name+
"_Violin_ProjectY").c_str(),
int(ybins.size()-1), &ybins[0] , 400, 0, MaxBinning);
362 ModeHist[i]->SetNameTitle((name+
"_mode").c_str(), (name+
"_mode").c_str());
364 ModeHist[i]->GetZaxis()->SetTitle(
"Mode");
366 W2MeanHist[i]->SetNameTitle((name+
"_w2mean").c_str(), (name+
"_w2mean").c_str());
368 W2MeanHist[i]->GetZaxis()->SetTitle(
"W2Mean");
370 W2ModeHist[i]->SetNameTitle((name+
"_w2mode").c_str(), (name+
"_w2mode").c_str());
372 W2ModeHist[i]->GetZaxis()->SetTitle(
"W2Mode");
375 lnLHist_Mean1D[i] =
new TH1D((name+
"_MeanlnL1D").c_str(),(name+
"_MeanlnL1D").c_str(), 50, 1, -1);
379 lnLHist_Mode1D[i] =
new TH1D((name+
"_ModelnL1D").c_str(),(name+
"_ModelnL1D").c_str(), 50, 1, -1);
397 for (
int i = 0; i < Length; ++i)
408void SampleSummary::AddThrow(std::vector<TH2Poly*> &SampleVector, std::vector<TH2Poly*> &W2Vec,
const double LLHPenalty,
const double Weight,
const int DrawNumber) {
414 const int size = int(SampleVector.size());
426 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
428 const int nXBins = 500;
430 std::string name = std::string(SampleVector[SampleNum]->GetName());
431 for (
int i = 1; i <=
maxBins[SampleNum]; ++i)
434 TH2PolyBin* bin =
static_cast<TH2PolyBin*
>(SampleVector[SampleNum]->GetBins()->At(i-1));
437 std::stringstream ss2;
439 ss2 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
440 ss2 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
442 PosteriorHist[SampleNum][i] = std::make_unique<TH1D>(ss2.str().c_str(), ss2.str().c_str(),nXBins, 1, -1);
444 w2Hist[SampleNum][i] = std::make_unique<TH1D>((
"w2_"+ss2.str()).c_str(), (
"w2_"+ss2.str()).c_str(),nXBins, 1, -1);
445 w2Hist[SampleNum][i]->SetDirectory(
nullptr);
448 std::string betaName =
"#beta_param_";
449 BetaHist[SampleNum][i] = std::make_unique<TH1D>((betaName + ss2.str()).c_str(), (betaName + ss2.str()).c_str(), 70, 1, -1);
450 BetaHist[SampleNum][i]->SetDirectory(
nullptr);
451 BetaHist[SampleNum][i]->GetXaxis()->SetTitle(
"#beta parameter value");
452 BetaHist[SampleNum][i]->GetYaxis()->SetTitle(
"Counts");
461 #pragma omp parallel for
463 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
465 if (SampleVector[SampleNum] ==
nullptr)
continue;
468 for (
int i = 1; i <=
maxBins[SampleNum]; ++i) {
469 const double Content = SampleVector[SampleNum]->GetBinContent(i);
471 const double w2 = W2Vec[SampleNum]->GetBinContent(i);
472 w2Hist[SampleNum][i]->Fill(w2, Weight);
475 const double data =
DataHist[SampleNum]->GetBinContent(i);
477 BetaHist[SampleNum][i]->Fill(BetaParam, Weight);
496 for (
int SampleNum = 0; SampleNum <
nSamples; SampleNum++)
498 if (
DataHist[SampleNum] ==
nullptr)
continue;
505 constexpr int nXBins = 500;
507 std::string name = std::string(
NominalHist[SampleNum]->GetName());
508 name = name.substr(0, name.find(
"_nom"));
511 for (
int i = 1; i <=
maxBins[SampleNum]; i++)
514 TH2PolyBin* bin =
static_cast<TH2PolyBin*
>(
NominalHist[SampleNum]->GetBins()->At(i-1));
517 std::stringstream ss2;
519 ss2 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
520 ss2 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
523 PosteriorHist_ByMode[SampleNum][j][i] =
new TH1D((name+ss2.str()).c_str(),(name+ss2.str()).c_str(),nXBins, 1, -1);
526 MeanHist_ByMode[SampleNum][j]->SetNameTitle((name+
"_mean").c_str(), (name+
"_mean").c_str());
535 #pragma omp parallel for
537 for (
int SampleNum = 0; SampleNum <
nSamples; SampleNum++)
539 if (
DataHist[SampleNum] ==
nullptr)
continue;
545 for (
int i = 1; i <=
maxBins[SampleNum]; ++i)
547 const double Content = SampleVector_ByMode[SampleNum][j]->GetBinContent(i);
561 TempString.replace(TempString.find(
".root"), 5, std::string(
"_procsW2.root"));
562 Outputfile =
new TFile(TempString.c_str(),
"RECREATE");
598 OutputTree =
new TTree(
"LLH_draws",
"LLH_draws");
606 while (SampleName.find(
" ") != std::string::npos) {
607 SampleName.replace(SampleName.find(
" "), 1, std::string(
"_"));
611 while (SampleName.find(
"-") != std::string::npos) {
612 SampleName.replace(SampleName.find(
"-"), 1, std::string(
"_"));
680 MACH3LOG_INFO(
"Made Prior/Posterior Predictive, it took {:.2f}s, now writing...", timer.RealTime());
710 RatioHistMean->GetZaxis()->SetTitle(
"Data/Mean");
712 RatioHistMode->GetZaxis()->SetTitle(
"Data/Mode");
714 RatioHistNom->GetZaxis()->SetTitle(
"Data/Nom");
726 TH1D *MeanHistCorrectedProjectX =
nullptr;
728 TH1D *MeanHistCorrectedProjectY =
nullptr;
759 while (SampleName.find(
"-") != std::string::npos) {
760 SampleName.replace(SampleName.find(
"-"), 1, std::string(
"_"));
762 OutputTree->Draw((SampleName+
"_data_draw:"+SampleName+
"_drawfluc_draw>>htemp").c_str());
763 TH2D *TempHistogram =
static_cast<TH2D*
>(gDirectory->Get(
"htemp")->Clone());
764 TempHistogram->GetXaxis()->SetTitle(
"-2LLH(Draw Fluc, Draw)");
765 TempHistogram->GetYaxis()->SetTitle(
"-2LLH(Data, Draw)");
766 TempHistogram->SetNameTitle((
SampleNames[i]+
"_drawfluc_draw").c_str(), (
SampleNames[i]+
"_drawfluc_draw").c_str());
768 TempHistogram->Write();
769 delete TempHistogram;
772 OutputTree->Draw((SampleName+
"_data_draw:"+SampleName+
"_predfluc_draw>>htemp2").c_str());
773 TH2D *TempHistogram2 =
static_cast<TH2D*
>(gDirectory->Get(
"htemp2")->Clone());
774 TempHistogram2->GetXaxis()->SetTitle(
"-2LLH(Pred Fluc, Draw)");
775 TempHistogram2->GetYaxis()->SetTitle(
"-2LLH(Data, Draw)");
776 TempHistogram2->SetNameTitle((
SampleNames[i]+
"_predfluc_draw").c_str(), (
SampleNames[i]+
"_predfluc_draw").c_str());
778 TempHistogram2->Write();
779 delete TempHistogram2;
782 OutputTree->Draw((SampleName+
"_rate_data_draw:"+SampleName+
"_rate_predfluc_draw>>htemp3").c_str());
783 TH2D *TempHistogram3 =
static_cast<TH2D*
>(gDirectory->Get(
"htemp3")->Clone());
784 TempHistogram3->GetXaxis()->SetTitle(
"-2LLH(Pred Fluc, Draw)");
785 TempHistogram3->GetYaxis()->SetTitle(
"-2LLH(Data, Draw)");
786 TempHistogram3->SetNameTitle((
SampleNames[i]+
"_rate_predfluc_draw").c_str(), (
SampleNames[i]+
"_rate_predfluc_draw").c_str());
788 TempHistogram3->Write();
789 delete TempHistogram3;
792 OutputTree->Draw((SampleName+
"_data_draw_ProjectX:"+SampleName+
"_drawfluc_draw_ProjectX>>htemp4").c_str());
793 TH2D *TempHistogram4 =
static_cast<TH2D*
>(gDirectory->Get(
"htemp4")->Clone());
796 TempHistogram4->SetNameTitle((
SampleNames[i]+
"_drawfluc_draw_ProjectX").c_str(), (
SampleNames[i]+
"_drawfluc_draw_ProjectX").c_str());
798 TempHistogram4->Write();
799 delete TempHistogram4;
806 RatioHistMean->Write();
807 RatioHistMode->Write();
808 RatioHistNom->Write();
815 DataNormHist->Write();
816 NomNormHist->Write();
817 MeanNormHist->Write();
818 ModeNormHist->Write();
821 NomProjectX->Write();
822 MeanProjectX->Write();
823 ModeProjectX->Write();
824 if(
DoBetaParam) MeanHistCorrectedProjectX->Write();
828 NomProjectY->Write();
829 MeanProjectY->Write();
830 ModeProjectY->Write();
831 if(
DoBetaParam) MeanHistCorrectedProjectY->Write();
834 W2NomProjectX->Write();
835 W2MeanProjectX->Write();
836 W2ModeProjectX->Write();
838 W2NomProjectY->Write();
839 W2MeanProjectY->Write();
840 W2ModeProjectY->Write();
845 TDirectory* DebugDir =
Dir[i]->mkdir(
"Debug");
847 for (
int b = 1; b <=
maxBins[i]; ++b)
854 TempLine->SetLineColor(kRed);
855 TempLine->SetLineWidth(2);
857 auto TempLineData = std::make_unique<TLine>(
DataHist[i]->GetBinContent(b),
PosteriorHist[i][b]->GetMinimum(),
859 TempLineData->SetLineColor(kGreen);
860 TempLineData->SetLineWidth(2);
865 Fitter->SetLineColor(kRed-5);
867 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
868 Legend->SetFillColor(0);
869 Legend->SetFillStyle(0);
870 Legend->SetLineWidth(0);
871 Legend->SetLineColor(0);
872 Legend->AddEntry(TempLineData.get(), Form(
"Data #mu=%.2f",
DataHist[i]->GetBinContent(b)),
"l");
873 Legend->AddEntry(TempLine.get(), Form(
"Prior #mu=%.2f",
NominalHist[i]->GetBinContent(b)),
"l");
875 Legend->AddEntry(Fitter, Form(
"Gauss, #mu=%.2f#pm%.2f", Fitter->GetParameter(1), Fitter->GetParameter(2)),
"l");
876 std::string TempTitle = std::string(
PosteriorHist[i][b]->GetName());
878 TempTitle +=
"_canv";
879 TCanvas *TempCanvas =
new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
880 TempCanvas->SetGridx();
881 TempCanvas->SetGridy();
882 TempCanvas->SetRightMargin(0.03);
883 TempCanvas->SetBottomMargin(0.08);
884 TempCanvas->SetLeftMargin(0.10);
885 TempCanvas->SetTopMargin(0.06);
888 TempLine->Draw(
"same");
889 TempLineData->Draw(
"same");
890 Fitter->Draw(
"same");
891 Legend->Draw(
"same");
925 MeanProjectX_ByMode->Write();
926 MeanProjectY_ByMode->Write();
930 for (
int b = 1; b <=
maxBins[i]; ++b)
935 delete MeanProjectX_ByMode;
936 delete MeanProjectY_ByMode;
940 delete RatioHistMean;
941 delete RatioHistMode;
961 delete W2NomProjectX;
962 delete W2MeanProjectX;
963 delete W2ModeProjectX;
965 delete W2NomProjectY;
966 delete W2MeanProjectY;
967 delete W2ModeProjectY;
981 double llh_total_temp = 0.0;
985 #pragma omp parallel for reduction(+:llh_total_temp)
987 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
993 double negLogL_Mean = 0.0;
994 double negLogL_Mode = 0.0;
997 for (
int j = 1; j <
maxBins[SampleNum]+1; ++j)
1000 TH1D *W2Projection =
w2Hist[SampleNum][j].get();
1003 const double nData =
DataHist[SampleNum]->GetBinContent(j);
1007 const double nMean = Projection->GetMean();
1008 const double nMeanError = Projection->GetRMS();
1009 const double nMode = Projection->GetBinCenter(Projection->GetMaximumBin());
1012 const double nW2Mean = W2Projection->GetMean();
1013 const double nW2Mode = W2Projection->GetBinCenter(W2Projection->GetMaximumBin());
1015 double TempLLH_Mean = 0.0;
1016 double TempLLH_Mode = 0.0;
1024 negLogL_Mean += 2*TempLLH_Mean;
1025 negLogL_Mode += 2*TempLLH_Mode;
1028 MeanHist[SampleNum]->SetBinContent(j,
MeanHist[SampleNum]->GetBinContent(j)+nMean);
1030 MeanHist[SampleNum]->SetBinError(j, nMeanError);
1034 TH1D *BetaTemp =
BetaHist[SampleNum][j].get();
1035 const double nBetaMean = BetaTemp->GetMean();
1036 const double nBetaMeanError = BetaTemp->GetRMS();
1040 const double ErrorTemp = std::sqrt( (nBetaMean*nMeanError) * (nBetaMean*nMeanError) + (nMean*nBetaMeanError) * (nMean*nBetaMeanError));
1045 ModeHist[SampleNum]->SetBinContent(j,
ModeHist[SampleNum]->GetBinContent(j)+nMode);
1047 ModeHist[SampleNum]->SetBinError(j, nModeError);
1055 lnLHist_Mean[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
1056 lnLHist_Mode[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mode);
1066 for (
int i = 1; i <
maxBins[SampleNum]+1; ++i)
1074 const double nMean = Projection->GetMean();
1075 const double nMeanError = Projection->GetRMS();
1084 llh_total_temp += negLogL_Mean;
1088 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1101 const double nMean = MeanProjectX->GetBinContent(j);
1102 const double nW2Mean = W2MeanProjectX->GetBinContent(j);
1104 double TempLLH_Mean = 0.0;
1111 delete MeanProjectX;
1112 delete W2MeanProjectX;
1118 std::stringstream ss;
1120 lnLHist->SetTitle((std::string(
lnLHist->GetTitle())+
"_"+ss.str()).c_str());
1140 double AveragePenalty = 0;
1143 std::vector<double> LLH_PredFluc_V(
nThrows);
1144 std::vector<double> LLH_DataDraw_V(
nThrows);
1145 std::vector<double> LLH_DrawFlucDraw_V(
nThrows);
1150 for (
unsigned int i = 0; i <
nThrows; ++i)
1157 double total_llh_data_draw_temp = 0.0;
1158 double total_llh_drawfluc_draw_temp = 0.0;
1159 double total_llh_predfluc_draw_temp = 0.0;
1161 double total_llh_rate_data_draw_temp = 0.0;
1162 double total_llh_rate_predfluc_draw_temp = 0.0;
1164 double total_llh_data_drawfluc_temp = 0.0;
1165 double total_llh_data_predfluc_temp = 0.0;
1166 double total_llh_draw_pred_temp = 0.0;
1167 double total_llh_drawfluc_pred_temp = 0.0;
1168 double total_llh_drawfluc_predfluc_temp = 0.0;
1169 double total_llh_predfluc_pred_temp = 0.0;
1170 double total_llh_datafluc_draw_temp = 0.0;
1172 double total_llh_data_draw_ProjectX_temp = 0.0;
1173 double total_llh_drawfluc_draw_ProjectX_temp = 0.0;
1180 std::vector<TH2Poly*> FluctHist(
nSamples);
1182 std::vector<TH2Poly*> FluctDrawHist(
nSamples);
1184 std::vector<TH2Poly*> DataFlucHist(
nSamples);
1187 std::vector<TH1D*> FluctDrawHistProjectX(
nSamples);
1188 std::vector<TH1D*> DrawHistProjectX(
nSamples);
1189 std::vector<TH1D*> DrawHistProjectY(
nSamples);
1190 std::vector<TH1D*> DrawW2HistProjectX(
nSamples);
1193 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1195 FluctHist[SampleNum] =
static_cast<TH2Poly*
>(
MeanHist[SampleNum]->Clone());
1196 FluctDrawHist[SampleNum] =
static_cast<TH2Poly*
>(
MeanHist[SampleNum]->Clone());
1197 DataFlucHist[SampleNum] =
static_cast<TH2Poly*
>(
MeanHist[SampleNum]->Clone());
1199 FluctDrawHistProjectX[SampleNum] =
static_cast<TH1D*
>(
DataHist_ProjectX[SampleNum]->Clone());
1202 TH2Poly *DrawHist =
MCVector[i][SampleNum];
1203 TH2Poly *DrawW2Hist =
W2MCVector[i][SampleNum];
1206 DrawHistProjectX[SampleNum] =
ProjectPoly(DrawHist,
true, SampleNum);
1207 DrawW2HistProjectX[SampleNum] =
ProjectPoly(DrawW2Hist,
true, SampleNum);
1208 DrawHistProjectY[SampleNum] =
ProjectPoly(DrawHist,
false, SampleNum);
1212 #pragma omp parallel for reduction(+:total_llh_data_draw_temp, total_llh_drawfluc_draw_temp, total_llh_predfluc_draw_temp, total_llh_rate_data_draw_temp, total_llh_rate_predfluc_draw_temp, total_llh_data_drawfluc_temp, total_llh_data_predfluc_temp, total_llh_draw_pred_temp, total_llh_drawfluc_pred_temp, total_llh_drawfluc_predfluc_temp, total_llh_predfluc_pred_temp, total_llh_datafluc_draw_temp, total_llh_data_draw_ProjectX_temp, total_llh_drawfluc_draw_ProjectX_temp)
1215 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1218 TH2Poly *DrawHist =
MCVector[i][SampleNum];
1219 TH2Poly *DrawW2Hist =
W2MCVector[i][SampleNum];
1221 if (DrawHist ==
nullptr)
continue;
1265 const double DataDrawLLH =
GetLLH(
DataHist[SampleNum], DrawHist, DrawW2Hist);
1267 total_llh_data_draw_temp += DataDrawLLH;
1270 const double DrawFlucDrawLLH =
GetLLH(FluctDrawHist[SampleNum], DrawHist, DrawW2Hist);
1272 total_llh_drawfluc_draw_temp += DrawFlucDrawLLH;
1275 const double PredFlucDrawLLH =
GetLLH(FluctHist[SampleNum], DrawHist, DrawW2Hist);
1277 total_llh_predfluc_draw_temp += PredFlucDrawLLH;
1283 total_llh_rate_data_draw_temp += RateDataDrawLLH;
1288 total_llh_rate_predfluc_draw_temp += RatePredFlucDrawLLH;
1292 const double DataDrawFlucLLH =
GetLLH(
DataHist[SampleNum], FluctDrawHist[SampleNum], DrawW2Hist);
1294 total_llh_data_drawfluc_temp += DataDrawFlucLLH;
1299 total_llh_data_predfluc_temp += DataPredFlucLLH;
1304 total_llh_draw_pred_temp += DrawPredLLH;
1309 total_llh_drawfluc_pred_temp += DrawFlucPredLLH;
1312 const double DrawFlucPredFlucLLH =
GetLLH(FluctDrawHist[SampleNum], FluctHist[SampleNum],
W2MeanHist[SampleNum]);
1314 total_llh_drawfluc_predfluc_temp += DrawFlucPredFlucLLH;
1319 total_llh_predfluc_pred_temp += PredFlucPredLLH;
1322 const double DataFlucDrawLLH =
GetLLH(DataFlucHist[SampleNum], DrawHist, DrawW2Hist);
1324 total_llh_datafluc_draw_temp += DataFlucDrawLLH;
1334 const double DataDrawLLH_ProjectX =
GetLLH(
DataHist_ProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1336 total_llh_data_draw_ProjectX_temp += DataDrawLLH_ProjectX;
1338 const double DrawFlucDrawLLH_ProjectX =
GetLLH(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1340 total_llh_drawfluc_draw_ProjectX_temp += DrawFlucDrawLLH_ProjectX;
1348 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1350 delete FluctHist[SampleNum];
1351 delete FluctDrawHist[SampleNum];
1352 delete DataFlucHist[SampleNum];
1353 delete FluctDrawHistProjectX[SampleNum];
1354 delete DrawHistProjectX[SampleNum];
1355 delete DrawHistProjectY[SampleNum];
1356 delete DrawW2HistProjectX[SampleNum];
1414 AveragePenalty = AveragePenalty/double(
nThrows);
1415 MACH3LOG_INFO(
"Average LLH penalty over toys is {:.2f}", AveragePenalty);
1417 unsigned int Accept_PredFluc = 0;
1418 unsigned int Accept_DrawFluc = 0;
1419 for (
unsigned int i = 0; i <
nThrows; ++i)
1421 if (LLH_DataDraw_V[i] > LLH_DrawFlucDraw_V[i]) Accept_DrawFluc++;
1422 if (LLH_DataDraw_V[i] > LLH_PredFluc_V[i]) Accept_PredFluc++;
1424 const double pvalue_DrawFluc = double(Accept_DrawFluc)/double(
nThrows);
1425 const double pvalue_PredFluc = double(Accept_PredFluc)/double(
nThrows);
1427 MACH3LOG_INFO(
"Calculated exact p-value using Fluctuation of Draw: {:.2f}", pvalue_DrawFluc);
1428 MACH3LOG_INFO(
"Calculated exact p-value using Fluctuation of Prediction: {:.2f}", pvalue_PredFluc);
1451 const double TotalIntegral = Histogram->Integral();
1454 double llh_reference = 0.0;
1456 llh_reference = llh_ref;
1460 for (
int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) {
1461 const double xvalue = Histogram->GetBinCenter(i+1);
1462 if (xvalue >= llh_reference) {
1463 Above += Histogram->GetBinContent(i+1);
1466 const double pvalue = Above/TotalIntegral;
1467 std::stringstream ss;
1468 ss << int(Above) <<
"/" << int(TotalIntegral) <<
"=" << pvalue;
1469 Histogram->SetTitle((std::string(Histogram->GetTitle())+
"_"+ss.str()).c_str());
1472 auto TempLine = std::make_unique<TLine>(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum());
1473 TempLine->SetLineColor(kBlack);
1474 TempLine->SetLineWidth(2);
1477 TH1D *TempHistogram =
static_cast<TH1D*
>(Histogram->Clone());
1478 TempHistogram->SetFillStyle(1001);
1479 TempHistogram->SetFillColor(kRed);
1480 for (
int i = 0; i < TempHistogram->GetNbinsX(); ++i) {
1481 if (TempHistogram->GetBinCenter(i+1) < llh_reference) {
1482 TempHistogram->SetBinContent(i+1, 0.0);
1486 auto Legend = std::make_unique<TLegend>(0.6, 0.6, 0.9, 0.9);
1487 Legend->SetFillColor(0);
1488 Legend->SetFillStyle(0);
1489 Legend->SetLineWidth(0);
1490 Legend->SetLineColor(0);
1491 Legend->AddEntry(TempLine.get(), Form(
"Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue),
"l");
1492 Legend->AddEntry(Histogram, Form(
"LLH, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()),
"l");
1493 std::string Title = Histogram->GetName();
1495 TCanvas *TempCanvas =
new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024);
1496 TempCanvas->SetGridx();
1497 TempCanvas->SetGridy();
1499 TempHistogram->Draw(
"same");
1500 TempLine->Draw(
"same");
1501 Legend->Draw(
"same");
1503 TempCanvas->Write();
1505 delete TempHistogram;
1514 auto TempLine = std::make_unique<TLine>(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum());
1515 TempLine->SetLineColor(kRed);
1516 TempLine->SetLineWidth(2);
1518 TF1 *Fitter =
new TF1(
"Fit",
"gaus", Histogram->GetBinLowEdge(1), Histogram->GetBinLowEdge(Histogram->GetNbinsX()+1));
1519 Histogram->Fit(Fitter,
"RQ");
1520 Fitter->SetLineColor(kRed-5);
1523 for (
int z = 0; z < Histogram->GetNbinsX(); ++z) {
1524 const double xvalue = Histogram->GetBinCenter(z+1);
1525 if (xvalue >= DataRate) {
1526 Above += Histogram->GetBinContent(z+1);
1529 const double pvalue = Above/Histogram->Integral();
1530 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1531 Legend->SetFillColor(0);
1532 Legend->SetFillStyle(0);
1533 Legend->SetLineWidth(0);
1534 Legend->SetLineColor(0);
1535 Legend->AddEntry(TempLine.get(), Form(
"Data, %.0f, p-value=%.2f", DataRate, pvalue),
"l");
1536 Legend->AddEntry(Histogram, Form(
"MC, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()),
"l");
1537 Legend->AddEntry(Fitter, Form(
"Gauss, #mu=%.1f#pm%.1f", Fitter->GetParameter(1), Fitter->GetParameter(2)),
"l");
1538 std::string TempTitle = std::string(Histogram->GetName());
1539 TempTitle +=
"_canv";
1540 TCanvas *TempCanvas =
new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1541 TempCanvas->SetGridx();
1542 TempCanvas->SetGridy();
1543 TempCanvas->SetRightMargin(0.03);
1544 TempCanvas->SetBottomMargin(0.08);
1545 TempCanvas->SetLeftMargin(0.10);
1546 TempCanvas->SetTopMargin(0.06);
1549 TempLine->Draw(
"same");
1550 Fitter->Draw(
"same");
1551 Legend->Draw(
"same");
1552 TempCanvas->Write();
1563 const double llh =
GetLLH(DatHist, MCHist, W2Hist);
1564 std::stringstream ss;
1565 ss <<
"_2LLH=" << llh;
1566 MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
1567 MACH3LOG_INFO(
"{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), DatHist->Integral(), MCHist->Integral(), llh);
1574 const double llh =
GetLLH(DatHist, MCHist, W2Hist);
1575 std::stringstream ss;
1576 ss <<
"_2LLH=" << llh;
1577 MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
1585 for (
int i = 1; i < DatHist->GetNumberOfBins()+1; ++i)
1587 const double data = DatHist->GetBinContent(i);
1588 const double mc = MCHist->GetBinContent(i);
1589 const double w2 = W2Hist->GetBinContent(i);
1600 for (
int i = 1; i <= DatHist->GetXaxis()->GetNbins(); ++i)
1602 const double data = DatHist->GetBinContent(i);
1603 const double mc = MCHist->GetBinContent(i);
1604 const double w2 = W2Hist->GetBinContent(i);
1615 TDirectory *BetaDir =
Outputfile->mkdir(
"BetaParameters");
1618 int originalErrorLevel = gErrorIgnoreLevel;
1621 gErrorIgnoreLevel = kFatal;
1624 std::vector<TDirectory *> DirBeta(
nSamples);
1628 DirBeta[i] = BetaDir->mkdir((
SampleNames[i]).c_str());
1632 for (
int j = 1; j <
maxBins[i]+1; ++j)
1634 const double data =
DataHist[i]->GetBinContent(j);
1635 const double mc =
NominalHist[i]->GetBinContent(j);
1636 const double w2 =
W2NomHist[i]->GetBinContent(j);
1640 auto TempLine = std::unique_ptr<TLine>(
new TLine(BetaPrior,
BetaHist[i][j]->GetMinimum(), BetaPrior,
BetaHist[i][j]->GetMaximum()));
1641 TempLine->SetLineColor(kRed);
1642 TempLine->SetLineWidth(2);
1645 TF1 *Fitter =
new TF1(
"Fit",
"gaus",
BetaHist[i][j]->GetBinLowEdge(1),
BetaHist[i][j]->GetBinLowEdge(
BetaHist[i][j]->GetNbinsX()+1));
1647 Fitter->SetLineColor(kRed-5);
1649 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1650 Legend->SetFillColor(0);
1651 Legend->SetFillStyle(0);
1652 Legend->SetLineWidth(0);
1653 Legend->SetLineColor(0);
1654 Legend->AddEntry(TempLine.get(), Form(
"Prior #mu=%.4f, N_{data}=%.0f", BetaPrior, data),
"l");
1655 Legend->AddEntry(
BetaHist[i][j].get(), Form(
"Post, #mu=%.4f#pm%.4f",
BetaHist[i][j]->GetMean(),
BetaHist[i][j]->GetRMS()),
"l");
1656 Legend->AddEntry(Fitter, Form(
"Gauss, #mu=%.4f#pm%.4f", Fitter->GetParameter(1), Fitter->GetParameter(2)),
"l");
1657 std::string TempTitle = std::string(
BetaHist[i][j]->GetName());
1659 TempTitle +=
"_canv";
1660 TCanvas *TempCanvas =
new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1661 TempCanvas->SetGridx();
1662 TempCanvas->SetGridy();
1663 TempCanvas->SetRightMargin(0.03);
1664 TempCanvas->SetBottomMargin(0.08);
1665 TempCanvas->SetLeftMargin(0.10);
1666 TempCanvas->SetTopMargin(0.06);
1669 TempLine->Draw(
"same");
1670 Fitter->Draw(
"same");
1671 Legend->Draw(
"same");
1672 TempCanvas->Write();
1678 DirBeta[i]->Write();
1684 gErrorIgnoreLevel = originalErrorLevel;
1697 std::vector<double> NEvents_Sample(
nSamples);
1698 double event_rate = 0.;
1701 TTree* Event_Rate_Tree =
new TTree(
"Event_Rate_draws",
"Event_Rate_draws");
1702 Event_Rate_Tree->Branch(
"Event_Rate", &event_rate);
1709 while (SampleName.find(
"-") != std::string::npos) {
1710 SampleName.replace(SampleName.find(
"-"), 1, std::string(
"_"));
1712 Event_Rate_Tree->Branch((SampleName+
"_Event_Rate").c_str(), &NEvents_Sample[i]);
1716 auto EventHist = std::make_unique<TH1D>(
"EventHist",
"Total Event Rate", 100, 1, -1);
1717 EventHist->SetDirectory(
nullptr);
1718 EventHist->GetXaxis()->SetTitle(
"Total event rate");
1719 EventHist->GetYaxis()->SetTitle(
"Counts");
1720 EventHist->SetLineWidth(2);
1723 std::vector<std::unique_ptr<TH1D>> SumHist(
nSamples);
1726 std::string name = std::string(
NominalHist[i]->GetName());
1727 name = name.substr(0, name.find(
"_nom"));
1729 SumHist[i] = std::make_unique<TH1D>((name+
"_sum").c_str(),(name+
"_sum").c_str(), 100, 1, -1);
1730 SumHist[i]->GetXaxis()->SetTitle(
"N_{events}");
1731 SumHist[i]->GetYaxis()->SetTitle(
"Counts");
1733 std::stringstream ss;
1735 SumHist[i]->SetTitle((std::string(SumHist[i]->GetTitle())+
"_"+ss.str()).c_str());
1738 for (
unsigned int it = 0; it <
nThrows; ++it)
1740 double event_rate_temp = 0.;
1743 #pragma omp parallel for reduction(+:event_rate_temp)
1745 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1749 SumHist[SampleNum]->Fill(NEvents_Sample[SampleNum],
WeightVector[it]);
1751 event_rate_temp += NEvents_Sample[SampleNum];
1753 event_rate = event_rate_temp;
1754 EventHist->Fill(event_rate);
1755 Event_Rate_Tree->Fill();
1757 Event_Rate_Tree->Write();
1758 delete Event_Rate_Tree;
1760 double DataRate = 0.0;
1762 #pragma omp parallel for reduction(+:DataRate)
1770 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1772 Dir[SampleNum]->cd();
1778 TDirectory *CorrDir =
Outputfile->mkdir(
"Correlations");
1781 TMatrixDSym* SampleCorrelation =
new TMatrixDSym(
nSamples);
1782 std::vector<std::vector<std::unique_ptr<TH2D>>> SamCorr(
nSamples);
1787 (*SampleCorrelation)(i,i) = 1.0;
1788 const double Min_i = SumHist[i]->GetXaxis()->GetBinLowEdge(1);
1789 const double Max_i = SumHist[i]->GetXaxis()->GetBinUpEdge(SumHist[i]->GetNbinsX()+1);
1792 const double Min_j = SumHist[j]->GetXaxis()->GetBinLowEdge(1);
1793 const double Max_j = SumHist[j]->GetXaxis()->GetBinUpEdge(SumHist[j]->GetNbinsX()+1);
1796 SamCorr[i][j] = std::make_unique<TH2D>(Form(
"SamCorr_%i_%i", i, j), Form(
"SamCorr_%i_%i", i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
1797 SamCorr[i][j]->SetDirectory(
nullptr);
1798 SamCorr[i][j]->SetMinimum(0);
1799 SamCorr[i][j]->GetXaxis()->SetTitle(
SampleNames[i].c_str());
1800 SamCorr[i][j]->GetYaxis()->SetTitle(
SampleNames[j].c_str());
1801 SamCorr[i][j]->GetZaxis()->SetTitle(
"Events");
1807 #pragma omp parallel for
1811 for (
int j = 0; j <= i; ++j)
1814 if (j == i)
continue;
1816 for (
unsigned int it = 0; it <
nThrows; ++it)
1820 SamCorr[i][j]->Smooth();
1823 (*SampleCorrelation)(i,j) = SamCorr[i][j]->GetCorrelationFactor();
1824 (*SampleCorrelation)(j,i) = (*SampleCorrelation)(i,j);
1829 hSamCorr->SetDirectory(
nullptr);
1830 hSamCorr->GetZaxis()->SetTitle(
"Correlation");
1831 hSamCorr->SetMinimum(-1);
1832 hSamCorr->SetMaximum(1);
1833 hSamCorr->GetXaxis()->SetLabelSize(0.015);
1834 hSamCorr->GetYaxis()->SetLabelSize(0.015);
1839 hSamCorr->GetXaxis()->SetBinLabel(i+1,
SampleNames[i].c_str());
1843 hSamCorr->GetYaxis()->SetBinLabel(j+1,
SampleNames[j].c_str());
1845 const double corr = (*SampleCorrelation)(i,j);
1846 hSamCorr->SetBinContent(i+1, j+1, corr);
1849 hSamCorr->Draw(
"colz");
1850 hSamCorr->Write(
"Sample_Corr");
1852 SampleCorrelation->Write(
"Sample_Correlation");
1853 delete SampleCorrelation;
1857 for (
int j = 0; j <= i; ++j)
1860 if (j == i)
continue;
1861 SamCorr[i][j]->Write();
1866 bool DoPerKinemBin =
false;
1870 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
1872 TMatrixDSym* KinCorrelation =
new TMatrixDSym(
maxBins[SampleNum]);
1873 std::vector<std::vector<std::unique_ptr<TH2D>>> KinCorr(
maxBins[SampleNum]);
1874 for (
int i = 0; i <
maxBins[SampleNum]; ++i)
1876 KinCorr[i].resize(
maxBins[SampleNum]);
1877 (*KinCorrelation)(i,i) = 1.0;
1879 const double Min_i =
PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinLowEdge(1);
1883 TH2PolyBin* bin =
static_cast<TH2PolyBin*
>(
NominalHist[SampleNum]->GetBins()->At(i));
1885 std::stringstream ss2;
1886 ss2 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
1887 ss2 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
1889 for (
int j = 0; j <
maxBins[SampleNum]; ++j)
1891 const double Min_j =
PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinLowEdge(1);
1895 KinCorr[i][j] = std::make_unique<TH2D>( Form(
"Kin_%i_%i_%i", SampleNum, i, j),
1896 Form(
"Kin_%i_%i_%i", SampleNum, i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
1897 KinCorr[i][j]->SetDirectory(
nullptr);
1898 KinCorr[i][j]->SetMinimum(0);
1900 KinCorr[i][j]->GetXaxis()->SetTitle(ss2.str().c_str());
1902 bin =
static_cast<TH2PolyBin*
>(
NominalHist[SampleNum]->GetBins()->At(j));
1904 std::stringstream ss3;
1905 ss3 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
1906 ss3 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
1907 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1908 KinCorr[i][j]->GetZaxis()->SetTitle(
"Events");
1913 #pragma omp parallel for
1915 for (
int i = 0; i <
maxBins[SampleNum]; ++i)
1917 for (
int j = 0; j <= i; ++j)
1920 if (j == i)
continue;
1922 for (
unsigned int it = 0; it <
nThrows; ++it)
1924 KinCorr[i][j]->Fill(
MCVector[it][SampleNum]->GetBinContent(i+1),
MCVector[it][SampleNum]->GetBinContent(j+1));
1926 KinCorr[i][j]->Smooth();
1929 (*KinCorrelation)(i,j) = KinCorr[i][j]->GetCorrelationFactor();
1930 (*KinCorrelation)(j,i) = (*KinCorrelation)(i,j);
1936 hKinCorr->SetDirectory(
nullptr);
1937 hKinCorr->GetZaxis()->SetTitle(
"Correlation");
1938 hKinCorr->SetMinimum(-1);
1939 hKinCorr->SetMaximum(1);
1940 hKinCorr->GetXaxis()->SetLabelSize(0.015);
1941 hKinCorr->GetYaxis()->SetLabelSize(0.015);
1944 for (
int i = 0; i <
maxBins[SampleNum]; ++i)
1947 TH2PolyBin* bin =
static_cast<TH2PolyBin*
>(
NominalHist[SampleNum]->GetBins()->At(i));
1949 std::stringstream ss2;
1950 ss2 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
1951 ss2 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
1952 hKinCorr->GetXaxis()->SetBinLabel(i+1, ss2.str().c_str());
1954 for (
int j = 0; j <
maxBins[SampleNum]; ++j)
1956 bin =
static_cast<TH2PolyBin*
>(
NominalHist[SampleNum]->GetBins()->At(j));
1958 std::stringstream ss3;
1959 ss3 <<
"p_{#mu} (" << bin->GetXMin() <<
"-" << bin->GetXMax() <<
")";
1960 ss3 <<
" cos#theta_{#mu} (" << bin->GetYMin() <<
"-" << bin->GetYMax() <<
")";
1961 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1963 hKinCorr->GetYaxis()->SetBinLabel(j+1, ss3.str().c_str());
1965 const double corr = (*KinCorrelation)(i,j);
1966 hKinCorr->SetBinContent(i+1, j+1, corr);
1969 hKinCorr->Draw(
"colz");
1970 hKinCorr->Write((
SampleNames[SampleNum] +
"_Corr").c_str());
1972 KinCorrelation->Write((
SampleNames[SampleNum] +
"_Correlation").c_str());
1973 delete KinCorrelation;
1990 MACH3LOG_INFO(
"Not calculating correlations per each kinematic bin");
2000 EventHist_ByMode[j] =
new TH1D(Form(
"EventHist_%s", ModeName.c_str()), Form(
"Total Event Rate %s", ModeName.c_str()), 100, 1, -1);
2001 EventHist_ByMode[j]->GetXaxis()->SetTitle(
"Total event rate");
2002 EventHist_ByMode[j]->GetYaxis()->SetTitle(
"Counts");
2003 EventHist_ByMode[j]->SetLineWidth(2);
2007 for (
unsigned int it = 0; it <
nThrows; ++it)
2011 double event_rate_temp = 0.;
2013 #pragma omp parallel for reduction(+:event_rate_temp)
2015 for (
int SampleNum = 0; SampleNum <
nSamples; SampleNum++)
2019 EventHist_ByMode[j]->Fill(event_rate_temp);
2028 TMatrixDSym* ModeCorrelation =
new TMatrixDSym(
Modes->
GetNModes()+1);
2035 (*ModeCorrelation)(i,i) = 1.0;
2037 const double Min_i = EventHist_ByMode[i]->GetXaxis()->GetBinLowEdge(1);
2038 const double Max_i = EventHist_ByMode[i]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[i]->GetNbinsX()+1);
2041 const double Min_j = EventHist_ByMode[j]->GetXaxis()->GetBinLowEdge(1);
2042 const double Max_j = EventHist_ByMode[j]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[j]->GetNbinsX()+1);
2045 ModeCorr[i][j] =
new TH2D(Form(
"ModeCorr_%i_%i",i,j), Form(
"ModeCorr_%i_%i",i,j), 70, Min_i, Max_i, 70, Min_j, Max_j);
2046 ModeCorr[i][j]->SetDirectory(
nullptr);
2047 ModeCorr[i][j]->SetMinimum(0);
2050 ModeCorr[i][j]->GetZaxis()->SetTitle(
"Events");
2056 #pragma omp parallel for
2060 for (
int j = 0; j <= i; ++j)
2063 if (j == i)
continue;
2065 for (
unsigned int it = 0; it <
nThrows; ++it)
2067 double Integral_X = 0.;
2068 double Integral_Y = 0.;
2069 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
2074 ModeCorr[i][j]->Fill(Integral_X, Integral_Y);
2076 ModeCorr[i][j]->Smooth();
2079 (*ModeCorrelation)(i,j) = ModeCorr[i][j]->GetCorrelationFactor();
2080 (*ModeCorrelation)(j,i) = (*ModeCorrelation)(i,j);
2085 hModeCorr->SetDirectory(
nullptr);
2086 hModeCorr->GetZaxis()->SetTitle(
"Correlation");
2087 hModeCorr->SetMinimum(-1);
2088 hModeCorr->SetMaximum(1);
2089 hModeCorr->GetXaxis()->SetLabelSize(0.015);
2090 hModeCorr->GetYaxis()->SetLabelSize(0.015);
2101 const double corr = (*ModeCorrelation)(i,j);
2102 hModeCorr->SetBinContent(i+1, j+1, corr);
2105 hModeCorr->Draw(
"colz");
2106 hModeCorr->Write(
"Mode_Corr");
2111 for (
int j = 0; j <= i; ++j)
2114 if (j == i)
continue;
2115 ModeCorr[i][j]->Write();
2123 delete ModeCorr[i][j];
2125 delete[] ModeCorr[i];
2128 ModeCorrelation->Write(
"Mode_Correlation");
2129 delete ModeCorrelation;
2133 delete EventHist_ByMode[j];
2141 MACH3LOG_INFO(
"Calculating correlations took {:.2f}s", timer.RealTime());
2149 TH1D* Projection =
nullptr;
2152 name = std::string(Histogram->GetName()) +
"_x";
2153 Projection = Histogram->ProjectionX(name.c_str(), 1, Histogram->GetYaxis()->GetNbins(),
"e");
2155 name = std::string(Histogram->GetName()) +
"_y";
2156 Projection = Histogram->ProjectionY(name.c_str(), 1, Histogram->GetXaxis()->GetNbins(),
"e");
2165 std::vector<double> xbins;
2166 std::vector<double> ybins;
2169 TH1D* Projection =
nullptr;
2172 name = std::string(Histogram->GetName()) +
"_x";
2173 Projection =
PolyProjectionX(Histogram, name.c_str(), xbins, MakeErrorHist);
2175 name = std::string(Histogram->GetName()) +
"_y";
2176 Projection =
PolyProjectionY(Histogram, name.c_str(), ybins, MakeErrorHist);
2230 double DataRate = 0.0;
2231 double BinsRate = 0.0;
2233 #pragma omp parallel for reduction(+:DataRate, BinsRate)
2237 if (
DataHist[i] ==
nullptr)
continue;
2244 MACH3LOG_INFO(
"Calculated Bayesian Information Criterion using global number of events: {:.2f}", EventRateBIC);
2245 MACH3LOG_INFO(
"Calculated Bayesian Information Criterion using global number of bins: {:.2f}", BinBasedBIC);
2257 #pragma omp parallel for reduction(+:Dbar)
2259 for (
unsigned int i = 0; i <
nThrows; ++i)
2261 double LLH_temp = 0.;
2262 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum)
2275 const double p_D = std::fabs(Dbar - Dhat);
2278 const double DIC_stat = Dhat + 2 * p_D;
2279 MACH3LOG_INFO(
"Effective number of parameters following DIC formalism is equal to: {:.2f}", p_D);
2293 #pragma omp parallel for reduction(+:lppd, p_WAIC)
2295 for (
int SampleNum = 0; SampleNum <
nSamples; ++SampleNum) {
2297 for (
int i = 1; i <=
nBins; ++i) {
2298 double mean_llh = 0.;
2299 double sum_exp_llh = 0;
2300 double mean_llh_squared = 0.;
2302 for (
unsigned int s = 0; s <
nThrows; ++s) {
2303 const double data =
DataHist[SampleNum]->GetBinContent(i);
2304 const double mc =
MCVector[s][SampleNum]->GetBinContent(i);
2305 const double w2 =
W2MCVector[s][SampleNum]->GetBinContent(i);
2311 double LLH_temp = -neg_LLH_temp;
2313 mean_llh += LLH_temp;
2314 mean_llh_squared += LLH_temp * LLH_temp;
2315 sum_exp_llh += std::exp(LLH_temp);
2322 sum_exp_llh = std::log(sum_exp_llh);
2325 lppd += sum_exp_llh;
2328 p_WAIC += mean_llh_squared - (mean_llh * mean_llh);
2333 double WAIC = -2 * (lppd - p_WAIC);
2334 MACH3LOG_INFO(
"Effective number of parameters following WAIC formalism is equal to: {:.2f}", p_WAIC);
TH1D * PolyProjectionX(TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
WP: Poly Projectors.
void MakeFluctuatedHistogramAlternative(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check.
void NormaliseTH2Poly(TH2Poly *Histogram)
Helper to Normalise histograms.
void MakeFluctuatedHistogramStandard(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using default fast method.
void FastViolinFill(TH2D *violin, TH1D *hist_1d)
KS: Fill Violin histogram with entry from a toy.
TH1D * PolyProjectionY(TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
WP: Poly Projectors.
double NoOverflowIntegral(TH2Poly *poly)
WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow.
TH2Poly * RatioPolys(TH2Poly *NumHist, TH2Poly *DenomHist)
Helper to make ratio of TH2Polys.
TH2Poly * NormalisePoly(TH2Poly *Histogram)
WP: Helper to Normalise histograms.
double GetBIC(const double llh, const int data, const int nPars)
Get the Bayesian Information Criterion (BIC) or Schwarz information criterion (also SIC,...
void Get2DBayesianpValue(TH2D *Histogram)
Calculates the 2D Bayesian p-value and generates a visualization.
double GetModeError(TH1D *hpost)
Get the mode error from a TH1D.
double GetBetaParameter(const double data, const double mc, const double w2, TestStatistic TestStat)
KS: Calculate Beta parameter which will be different based on specified test statistic.
Custom exception class for MaCh3 errors.
int GetNModes() const
KS: Get number of modes, keep in mind actual number is +1 greater due to unknown category.
std::string GetMaCh3ModeName(const int Index)
KS: Get normal name of mode, if mode not known you will get UNKNOWN_BAD.
Class responsible for handling implementation of samples used in analysis, reweighting and returning ...
virtual void SetupBinning(const M3::int_t Selection, std::vector< double > &BinningX, std::vector< double > &BinningY)
virtual std::string GetKinVarLabel(const int sample, const int Dimension)
void AddThrowByMode(std::vector< std::vector< TH2Poly * > > &SampleVector_ByMode)
KS: Add histograms for each mode.
void StudyInformationCriterion(M3::kInfCrit Criterion)
Information Criterion.
std::vector< TH2Poly * > lnLHist_Mode
The LLH distribution in pmu cosmu for using the mode in each bin.
std::vector< double > LLHPenaltyVector
Vector to hold the penalty term.
std::vector< double > llh_data_draw_ProjectX
Projection X (most likely muon momentum) of LLH.
void MakeCutEventRate(TH1D *Histogram, const double DataRate)
Make the 1D Event Rate Hist.
SampleHandlerBase * SampleHandler
Pointer to SampleHandler object, mostly used to get sample names, binning etc.
~SampleSummary()
Destructor.
std::vector< double > llh_predfluc_draw
Fluctuated Predictive vs Draw.
void AddThrow(std::vector< TH2Poly * > &MCHist, std::vector< TH2Poly * > &W2Hist, const double LLHPenalty=0.0, const double Weight=1.0, const int DrawNumber=0)
KS: Add histograms with throws.
std::vector< std::vector< TH2Poly * > > MCVector
Vector of vectors which holds the loaded MC histograms.
std::vector< TH2Poly * > lnLHist_Mean
The LLH distribution in pmu cosmu for using the mean in each bin.
void CalcLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly, will modify MC hist tittle to include LLH.
double total_llh_predfluc_pred
Fluctuated Predictive vs Predictive.
double total_llh_drawfluc_draw_ProjectX
Fluctuated Draw vs Draw for projection X (most likely muon momentum)
std::vector< std::vector< TH2Poly * > > W2MCVector
Vector of vectors which holds the loaded W2 histograms.
std::vector< TH2Poly * > W2MeanHist
Pointer to the w2 histograms (for mean values).
double total_llh_data_draw
Data vs Draw.
std::unique_ptr< TH1D > RandomHist
Holds the history of which entries have been drawn in the MCMC file.
double total_llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using Rate.
void MakeChi2Hists()
Make the fluctuated histograms (2D and 1D) for the chi2s Essentially taking the MCMC draws and calcul...
std::vector< double > llh_data_predfluc
Data vs Fluctuated Predictive.
void StudyBIC()
Study Bayesian Information Criterion (BIC) .
std::unique_ptr< TRandom3 > rnd
Random number generator.
std::vector< TH1D * > lnLHist_Mode1D
Holds the bin-by-bin LLH for the mode posterior predictive vs the data.
std::string OutputName
Output filename.
void MakePredictive()
Finalise the distributions from the thrown samples.
std::vector< double > llh_predfluc_pred
Fluctuated Predictive vs Predictive.
std::vector< std::vector< std::unique_ptr< TH1D > > > BetaHist
Distribution of beta parameters in Barlow Beeston formalisms.
int nModelParams
Number of parameters.
std::vector< TH2D * > ViolinHists_ProjectX
Posterior predictive but for X projection but as a violin plot.
bool first_pass
KS: Hacky flag to let us know if this is first toy.
std::vector< TH2Poly * > NominalHist
The nominal histogram for the selection.
std::vector< std::vector< std::unique_ptr< TH1D > > > PosteriorHist
The posterior predictive for the whole selection: this gets built after adding in the toys....
double total_llh_draw_pred
Draw vs Predictive.
std::vector< double > llh_drawfluc_draw
Fluctuated Draw vs Draw.
double total_llh_data_draw_ProjectX
Data vs Draw for projection X (most likely muon momentum)
std::vector< TH2Poly * > MeanHistCorrected
The posterior predictive distribution in pmu cosmu using the mean after applying Barlow-Beeston Corre...
std::vector< double > llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
MaCh3Modes * Modes
MaCh3 Modes.
std::vector< double > llh_data_draw
Data vs Draw.
double total_llh_data_drawfluc
Data vs Fluctuated Draw.
void AddData(std::vector< TH2Poly * > &DataHist)
KS: Add data histograms.
TH1D * ProjectPoly(TH2Poly *Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist=false)
Helper to project TH2Poly onto axis.
std::unique_ptr< TH1D > lnLHist
The histogram containing the lnL for each throw.
std::vector< TH1D * > lnLHist_Mean1D
Holds the bin-by-bin LLH for the mean posterior predictive vs the data.
TTree * OutputTree
TTree which we save useful information to.
bool DoByModePlots
By mode variables.
std::vector< TH1D * > lnLHist_Mean_ProjectX
The LLH distribution in pmu using the mean in each bin.
void PrepareOutput()
KS: Prepare output tree and necessary variables.
SampleSummary(const int n_Samples, const std::string &Filename, SampleHandlerBase *const sample, const int nSteps)
Constructor.
double total_llh_drawfluc_draw
Fluctuated Draw vs Draw.
std::vector< TH2Poly * > DataHist
The data histogram for the selection.
std::vector< double > llh_draw_pred
Draw vs Predictive.
double total_llh_predfluc_draw
Fluctuated Predictive vs Draw.
int Debug
Tells Debug level to save additional histograms.
std::vector< double > llh_drawfluc_draw_ProjectX
bool CheckSamples(const int Length)
Check the length of samples agrees.
std::unique_ptr< TH2D > lnLDrawHist
The 2D lnLhist, showing (draw vs data) and (draw vs fluct), anything above y=x axis is the p-value.
void AddNominal(std::vector< TH2Poly * > &NominalHist, std::vector< TH2Poly * > &W2Nom)
KS: Add prior histograms.
void StudyDIC()
KS: Get the Deviance Information Criterion (DIC) .
std::vector< TH2D * > ViolinHists_ProjectY
Posterior predictive but for Y projection but as a violin plot.
std::vector< double > llh_data_drawfluc
Data vs Fluctuated Draw.
std::vector< TH2Poly * > ModeHist
The posterior predictive distribution in pmu cosmu using the mode.
std::vector< std::string > SampleNames
name for each sample
std::vector< TH2Poly * > W2NomHist
Pointer to the w2 histograms (for nominal values).
double total_llh_drawfluc_pred
Fluctuated Draw vs Predictive.
void StudyWAIC()
KS: Get the Watanabe-Akaike information criterion (WAIC) .
std::vector< TDirectory * > Dir
Directory for each sample.
std::vector< double > llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using rate only.
unsigned int nThrows
Number of throws.
std::vector< TH1D * > lnLHist_Sample_DrawData
The histogram containing the lnL (draw vs data) for each throw for each sample.
std::vector< std::vector< std::unique_ptr< TH1D > > > w2Hist
The posterior predictive for the whole selection: this gets built after adding in the toys....
TH1D **** PosteriorHist_ByMode
Histogram which corresponds to each bin in the sample's th2poly.
std::unique_ptr< TH1D > lnLHist_drawflucdraw
The lnLhist for the draw vs draw fluctuated.
TFile * Outputfile
Output filename.
double llh_total
Total LLH for the posterior predictive distribution.
void Write()
KS: Write results into root file.
std::unique_ptr< TH2D > lnLDrawHistRate
The 2D lnLhist, showing (draw vs data) and (draw vs fluct), using rate, anything above y=x axis is th...
double total_llh_rate_data_draw
Rate Data vs Draw.
std::unique_ptr< TH1D > lnLHist_drawdata
The lnLhist for the draw vs data.
TestStatistic likelihood
Type of likelihood for example Poisson, Barlow-Beeston or Ice Cube.
int nSamples
Number of samples.
void MakeCutLLH()
Make the cut LLH histogram.
std::vector< std::vector< std::vector< TH2Poly * > > > MCVectorByMode
Vector of vectors which holds the loaded MC histograms for each mode.
std::vector< double > llh_drawfluc_pred
Fluctuated Draw vs Predictive.
std::vector< std::vector< TH2Poly * > > MeanHist_ByMode
The posterior predictive distribution in pmu cosmu using the mean.
std::vector< double > llh_datafluc_draw
Fluctuated Data vs Draw.
std::vector< TH1D * > DataHist_ProjectX
The data histogram for the selection X projection.
void MakeCutLLH1D(TH1D *Histogram, double llh_ref=-999)
std::vector< TH1D * > DataHist_ProjectY
The data histogram for the selection Y projection.
void PlotBetaParameters()
KS: In Barlow Beeston we have Beta Parameters which scale generated MC.
double GetLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly.
std::vector< double > WeightVector
Vector holding weight.
void MakeFluctuatedHistogram(TH1D *FluctHist, TH1D *PolyHist)
Make Poisson fluctuation of TH1D hist.
std::vector< double > llh_rate_data_draw
Data vs Draw using rate only.
bool isPriorPredictive
bool whether we have Prior or Posterior Predictive
bool DoBetaParam
Are we making Beta Histograms.
std::vector< TH2Poly * > W2ModeHist
Pointer to the w2 histograms (for mode values).
std::unique_ptr< TH2D > lnLFlucHist
The 2D lnLHist, showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-val...
bool StandardFluctuation
KS: We have two methods for Poissonian fluctuation.
std::unique_ptr< TH2D > lnLFlucHist_ProjectX
The 2D lnLHist but for ProjectionX histogram (pmu), showing (draw vs data) and (draw vs draw fluct),...
std::unique_ptr< TH1D > lnLHist_drawfluc
The lnLhist for the draw vs MC fluctuated.
bool doShapeOnly
bool whether to normalise each toy to have shape based p-value and pos pred distribution
std::vector< TH1D * > lnLHist_Sample_PredflucDraw
The histogram containing the lnL (draw vs pred fluct) for each throw for each sample.
double total_llh_data_predfluc
Data vs Fluctuated Predictive.
double llh_penalty
LLH penalty for each throw.
std::vector< TH2Poly * > MeanHist
The posterior predictive distribution in pmu cosmu using the mean.
unsigned int nChainSteps
Number of throws by user.
std::vector< int > maxBins
Max Number of Bins per each sample.
TH1D * ProjectHist(TH2D *Histogram, const bool ProjectX)
Helper to project TH2D onto axis.
double total_llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
double total_llh_datafluc_draw
Fluctuated Data vs Draw.
void StudyKinematicCorrelations()
KS: Study how correlated are sample or kinematic bins.
std::vector< TH1D * > lnLHist_Sample_DrawflucDraw
The histogram containing the lnL (draw vs draw fluct) for each throw for each sample.
virtual std::string GetSampleName(int Sample) const =0
MaCh3Modes * GetMaCh3Modes() const
Return pointer to MaCh3 modes.
double GetTestStatLLH(double data, double mc) const
Calculate test statistic for a single bin using Poisson.
kInfCrit
KS: Different Information Criterion tests mostly based Gelman paper.
@ kWAIC
Watanabe-Akaike information criterion.
@ kInfCrits
This only enumerates.
@ kBIC
Bayesian Information Criterion.
@ kDIC
Deviance Information Criterion.
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.