7#include "TGraphAsymmErrors.h"
10#pragma GCC diagnostic ignored "-Wuseless-cast"
18 random = std::make_unique<TRandom3>(Get<int>(
fitMan->
raw()[
"General"][
"Seed"], __FILE__, __LINE__));
25 clock = std::make_unique<TStopwatch>();
26 stepClock = std::make_unique<TStopwatch>();
29 debug = GetFromManager<bool>(
fitMan->
raw()[
"General"][
"Debug"],
false, __FILE__ , __LINE__);
32 auto outfile = Get<std::string>(
fitMan->
raw()[
"General"][
"OutputFile"], __FILE__ , __LINE__);
35 auto_save = Get<int>(
fitMan->
raw()[
"General"][
"MCMC"][
"AutoSave"], __FILE__ , __LINE__);
46 outTree =
new TTree(
"posteriors",
"Posterior_Distributions");
62 if (debug) debugFile.open((outfile+
".log").c_str());
66 fTestLikelihood = GetFromManager<bool>(
fitMan->
raw()[
"General"][
"Fitter"][
"FitTestLikelihood"],
false, __FILE__ , __LINE__);
88 TDirectory* MaCh3Version =
outputFile->mkdir(
"MaCh3Engine");
91 if (std::getenv(
"MaCh3_ROOT") == NULL) {
97 if (std::getenv(
"MACH3") == NULL) {
102 std::string header_path = std::string(std::getenv(
"MACH3"));
103 header_path +=
"/version.h";
104 FILE* file = fopen(header_path.c_str(),
"r");
107 header_path = std::string(std::getenv(
"MaCh3_ROOT"));
108 header_path +=
"/version.h";
114 TMacro versionHeader(
"version_header",
"version_header");
115 versionHeader.ReadFile(header_path.c_str());
116 versionHeader.Write();
119 Engine.Write(
GetName().c_str());
121 MaCh3Version->Write();
133 for(
unsigned int i = 0; i <
systematics.size(); ++i)
136 for(
unsigned int i = 0; i <
samples.size(); ++i)
161 bool SaveProposal = GetFromManager<bool>(
fitMan->
raw()[
"General"][
"SaveProposal"],
false, __FILE__ , __LINE__);
165 cov->SetBranches(*
outTree, SaveProposal);
178 for (
size_t i = 0; i <
samples.size(); ++i) {
179 std::stringstream oss, oss2;
180 oss <<
"LogL_sample_" << i;
181 oss2 << oss.str() <<
"/D";
186 std::stringstream oss, oss2;
187 oss <<
"LogL_systematic_" <<
systematics[i]->GetName();
188 oss2 << oss.str() <<
"/D";
200 MACH3LOG_INFO(
"-------------------- Starting MCMC --------------------");
203 debugFile <<
"----- Starting MCMC -----" << std::endl;
215 for (
size_t i = 0; i <
samples.size(); ++i) {
216 samples[i]->CleanMemoryBeforeFit();
235 debugFile <<
"\n\n" <<
step <<
" steps took " <<
clock->RealTime() <<
" seconds to complete. (" <<
clock->RealTime() /
step <<
"s / step).\n" <<
accCount<<
" steps were accepted." << std::endl;
249 for (
const auto &s :
samples) {
250 if (s->GetTitle() == sample->
GetTitle()) {
280 CorrMatrix->Write((cov->
GetName() + std::string(
"_Corr")).c_str());
299 TFile *infile =
new TFile(FitName.c_str(),
"READ");
300 TTree *posts = infile->Get<TTree>(
"posteriors");
303 posts->SetBranchAddress(
"step",&step_val);
304 posts->SetBranchAddress(
"LogL",&log_val);
308 TDirectory* CovarianceFolder = infile->Get<TDirectory>(
"CovarianceFolder");
310 std::string ConfigName =
"Config_" +
systematics[s]->GetName();
311 TMacro *ConfigCov = CovarianceFolder->Get<TMacro>(ConfigName.c_str());
313 if (ConfigCov !=
nullptr) {
317 YAML::Node ConfigCurrent =
systematics[s]->GetConfig();
321 MACH3LOG_ERROR(
"Yaml configs in previous chain and current one are different", FitName);
328 CovarianceFolder->Close();
329 delete CovarianceFolder;
332 for (
int i = 0; i <
systematics[s]->GetNumParams(); ++i) {
333 posts->SetBranchAddress(
systematics[s]->GetParName(i).c_str(), &branch_vals[i]);
335 posts->GetEntry(posts->GetEntries()-1);
337 for (
int i = 0; i <
systematics[s]->GetNumParams(); ++i) {
340 MACH3LOG_ERROR(
"Parameter {} is unvitalised with value {}", i, branch_vals[i]);
341 MACH3LOG_ERROR(
"Please check more precisely chain you passed {}", FitName);
352 for (
int i = 0; i <
systematics[s]->GetNumParams(); ++i) {
353 posts->SetBranchAddress(
systematics[s]->GetParName(i).c_str(),
nullptr);
373 if (
fitMan ==
nullptr)
return;
376 if (GetFromManager<bool>(
fitMan->
raw()[
"General"][
"ProcessMCMC"],
false, __FILE__ , __LINE__)){
382 TVectorD *Central =
nullptr;
383 TVectorD *Errors =
nullptr;
384 TVectorD *Central_Gauss =
nullptr;
385 TVectorD *Errors_Gauss =
nullptr;
386 TVectorD *Peaks =
nullptr;
389 Processor.
GetPostfit(Central, Errors, Central_Gauss, Errors_Gauss, Peaks);
393 TMatrixDSym *Covariance =
nullptr;
394 TMatrixDSym *Correlation =
nullptr;
404 MACH3LOG_INFO(
"Opening output again to update with means..");
405 outputFile =
new TFile(Get<std::string>(
fitMan->
raw()[
"General"][
"Output"][
"Filename"], __FILE__, __LINE__).c_str(),
"UPDATE");
407 Central->Write(
"PDF_Means");
408 Errors->Write(
"PDF_Errors");
409 Central_Gauss->Write(
"Gauss_Means");
410 Errors_Gauss->Write(
"Errors_Gauss");
411 Covariance->Write(
"Covariance");
412 Correlation->Write(
"Correlation");
424 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs)
426 TStopwatch clockRace;
428 for(
int Lap = 0; Lap < NLaps; ++Lap) {
432 MACH3LOG_INFO(
"It took {:.4f} s to reweights {} times sample: {}", clockRace.RealTime(), NLaps,
samples[ivs]->GetTitle());
433 MACH3LOG_INFO(
"On average {:.6f}", clockRace.RealTime()/NLaps);
436 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs)
438 TStopwatch clockRace;
440 for(
int Lap = 0; Lap < NLaps; ++Lap) {
444 MACH3LOG_INFO(
"It took {:.4f} s to calculate GetLikelihood {} times sample: {}", clockRace.RealTime(), NLaps,
samples[ivs]->GetTitle());
445 MACH3LOG_INFO(
"On average {:.6f}", clockRace.RealTime()/NLaps);
448 std::vector<std::vector<double>> StepsValuesBefore(
systematics.size());
450 StepsValuesBefore[s] =
systematics[s]->GetProposed();
453 TStopwatch clockRace;
455 for(
int Lap = 0; Lap < NLaps; ++Lap) {
459 MACH3LOG_INFO(
"It took {:.4f} s to propose step {} times cov: {}", clockRace.RealTime(), NLaps,
systematics[s]->GetName());
460 MACH3LOG_INFO(
"On average {:.6f}", clockRace.RealTime()/NLaps);
463 systematics[s]->SetParameters(StepsValuesBefore[s]);
467 TStopwatch clockRace;
469 for(
int Lap = 0; Lap < NLaps; ++Lap) {
473 MACH3LOG_INFO(
"It took {:.4f} s to calculate get likelihood {} times cov: {}", clockRace.RealTime(), NLaps,
systematics[s]->GetName());
474 MACH3LOG_INFO(
"On average {:.6f}", clockRace.RealTime()/NLaps);
482 bool isScanRanges =
false;
484 if(
fitMan->
raw()[
"LLHScan"][
"ScanRanges"]){
485 YAML::Node scanRangesList =
fitMan->
raw()[
"LLHScan"][
"ScanRanges"];
486 for (
auto it = scanRangesList.begin(); it != scanRangesList.end(); ++it) {
487 std::string itname = it->first.as<std::string>();
488 std::vector<double> itrange = it->second.as<std::vector<double>>();
490 scanRanges[itname] = itrange;
494 MACH3LOG_INFO(
"There are no user-defined parameter ranges, so I'll use default param bounds for LLH Scans");
503 for(
unsigned int is = 0; is < SkipVector.size(); ++is)
505 if(ParamName.substr(0, SkipVector[is].length()) == SkipVector[is])
524 bool PlotLLHScanBySample = GetFromManager<bool>(
fitMan->
raw()[
"LLHScan"][
"LLHScanBySample"],
false, __FILE__ , __LINE__);;
525 auto SkipVector = GetFromManager<std::vector<std::string>>(
fitMan->
raw()[
"LLHScan"][
"LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
529 std::vector<TDirectory *> Cov_LLH(
systematics.size());
530 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc )
532 std::string NameTemp =
systematics[ivc]->GetName();
533 NameTemp = NameTemp.substr(0, NameTemp.find(
"_cov")) +
"_LLH";
534 Cov_LLH[ivc] =
outputFile->mkdir(NameTemp.c_str());
537 std::vector<TDirectory *> SampleClass_LLH(
samples.size());
538 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
540 std::string NameTemp =
samples[ivs]->GetTitle();
541 SampleClass_LLH[ivs] =
outputFile->mkdir(NameTemp.c_str());
544 TDirectory *Sample_LLH =
outputFile->mkdir(
"Sample_LLH");
545 TDirectory *Total_LLH =
outputFile->mkdir(
"Total_LLH");
547 std::vector<TDirectory *>SampleSplit_LLH;
548 if(PlotLLHScanBySample)
551 int SampleIterator = 0;
552 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
554 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is )
556 SampleSplit_LLH[SampleIterator] =
outputFile->mkdir((
samples[ivs]->GetSampleName(is)+
"_LLH").c_str());
562 const int n_points = GetFromManager<int>(
fitMan->
raw()[
"LLHScan"][
"LLHScanPoints"], 100, __FILE__ , __LINE__);
563 double nSigma = GetFromManager<int>(
fitMan->
raw()[
"LLHScan"][
"LLHScanSigma"], 1., __FILE__, __LINE__);
566 const int countwidth = int(
double(n_points)/
double(5));
569 std::map<std::string, std::vector<double>> scanRanges;
577 int npars = cov->GetNumParams();
578 bool IsPCA = cov->IsPCA();
579 if (IsPCA) npars = cov->GetNParameters();
580 for (
int i = 0; i < npars; ++i)
583 std::string name = cov->GetParFancyName(i);
584 if (IsPCA) name +=
"_PCA";
589 double prior = cov->GetParInit(i);
590 if (IsPCA) prior = cov->GetPCAHandler()->GetParCurrPCA(i);
595 double lower = prior - nSigma*cov->GetDiagonalError(i);
596 double upper = prior + nSigma*cov->GetDiagonalError(i);
599 lower = prior - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
600 upper = prior + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
601 MACH3LOG_INFO(
"eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
611 auto it = scanRanges.find(name);
612 if (it != scanRanges.end() && it->second.size() == 2) {
613 lower = it->second[0];
614 upper = it->second[1];
615 MACH3LOG_INFO(
"Found matching param name for setting specified range for {}", name);
616 MACH3LOG_INFO(
"Range for {} = [{:.2f}, {:.2f}]", name, lower, upper);
622 lower = std::max(lower, cov->GetLowerBound(i));
623 upper = std::min(upper, cov->GetUpperBound(i));
624 MACH3LOG_INFO(
"Scanning {} with {} steps, from [{:.2f} , {:.2f}], prior = {:.2f}", name, n_points, lower, upper, prior);
627 auto hScan = std::make_unique<TH1D>((name +
"_full").c_str(), (name +
"_full").c_str(), n_points, lower, upper);
628 hScan->SetTitle((std::string(
"2LLH_full, ") + name +
";" + name +
"; -2(ln L_{sample} + ln L_{xsec+flux} + ln L_{det})").c_str());
629 hScan->SetDirectory(
nullptr);
631 auto hScanSam = std::make_unique<TH1D>((name +
"_sam").c_str(), (name +
"_sam").c_str(), n_points, lower, upper);
632 hScanSam->SetTitle((std::string(
"2LLH_sam, ") + name +
";" + name +
"; -2(ln L_{sample})").c_str());
633 hScanSam->SetDirectory(
nullptr);
635 std::vector<std::unique_ptr<TH1D>> hScanSample(
samples.size());
636 std::vector<double> nSamLLH(
samples.size());
637 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
639 std::string NameTemp =
samples[ivs]->GetTitle();
640 hScanSample[ivs] = std::make_unique<TH1D>((name+
"_"+NameTemp).c_str(), (name+
"_" + NameTemp).c_str(), n_points, lower, upper);
641 hScanSample[ivs]->SetDirectory(
nullptr);
642 hScanSample[ivs]->SetTitle((
"2LLH_" + NameTemp +
", " + name +
";" + name +
"; -2(ln L_{" + NameTemp +
"})").c_str());
646 std::vector<std::unique_ptr<TH1D>> hScanCov(
systematics.size());
648 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc )
650 std::string NameTemp =
systematics[ivc]->GetName();
651 NameTemp = NameTemp.substr(0, NameTemp.find(
"_cov"));
652 hScanCov[ivc] = std::make_unique<TH1D>((name +
"_" + NameTemp).c_str(), (name +
"_" + NameTemp).c_str(), n_points, lower, upper);
653 hScanCov[ivc]->SetDirectory(
nullptr);
654 hScanCov[ivc]->SetTitle((
"2LLH_" + NameTemp +
", " + name +
";" + name +
"; -2(ln L_{" + NameTemp +
"})").c_str());
658 std::vector<TH1D*> hScanSamSplit;
659 std::vector<double> sampleSplitllh;
660 if(PlotLLHScanBySample)
662 int SampleIterator = 0;
663 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
667 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is )
669 hScanSamSplit[SampleIterator] =
new TH1D((name+
samples[ivs]->GetSampleName(is)).c_str(), (name+
samples[ivs]->GetSampleName(is)).c_str(), n_points, lower, upper);
670 hScanSamSplit[SampleIterator]->SetTitle((std::string(
"2LLH_sam, ") + name +
";" + name +
"; -2(ln L_{sample})").c_str());
677 for (
int j = 0; j < n_points; ++j)
679 if (j % countwidth == 0)
684 cov->GetPCAHandler()->SetParPropPCA(i, hScan->GetBinCenter(j+1));
687 cov->SetParProp(i, hScan->GetBinCenter(j+1));
691 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
699 double samplellh = 0;
701 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
703 nSamLLH[ivs] =
samples[ivs]->GetLikelihood();
704 samplellh += nSamLLH[ivs];
707 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc )
710 totalllh += nCovLLH[ivc];
713 totalllh += samplellh;
715 if(PlotLLHScanBySample)
717 int SampleIterator = 0;
718 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
720 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is)
722 sampleSplitllh[SampleIterator] =
samples[ivs]->GetSampleLikelihood(is);
728 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs ) {
729 hScanSample[ivs]->SetBinContent(j+1, 2*nSamLLH[ivs]);
731 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc ) {
732 hScanCov[ivc]->SetBinContent(j+1, 2*nCovLLH[ivc]);
735 hScanSam->SetBinContent(j+1, 2*samplellh);
736 hScan->SetBinContent(j+1, 2*totalllh);
738 if(PlotLLHScanBySample)
740 int SampleIterator = 0;
741 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
743 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is)
745 hScanSamSplit[SampleIterator]->SetBinContent(j+1, 2*sampleSplitllh[SampleIterator]);
751 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc )
754 hScanCov[ivc]->Write();
757 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
759 SampleClass_LLH[ivs]->cd();
760 hScanSample[ivs]->Write();
767 if(PlotLLHScanBySample)
769 int SampleIterator = 0;
770 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
772 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is)
774 SampleSplit_LLH[SampleIterator]->cd();
775 hScanSamSplit[SampleIterator]->Write();
776 delete hScanSamSplit[SampleIterator];
784 cov->GetPCAHandler()->SetParPropPCA(i, prior);
786 cov->SetParProp(i, prior);
791 for(
unsigned int ivc = 0; ivc <
systematics.size(); ++ivc )
793 Cov_LLH[ivc]->Write();
797 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
799 SampleClass_LLH[ivs]->Write();
800 delete SampleClass_LLH[ivs];
809 if(PlotLLHScanBySample)
811 int SampleIterator = 0;
812 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs )
814 for(
int is = 0; is <
samples[ivs]->GetNsamples(); ++is )
816 SampleSplit_LLH[SampleIterator]->Write();
817 delete SampleSplit_LLH[SampleIterator];
828 TDirectory *Sample_LLH =
outputFile->Get<TDirectory>(
"Sample_LLH");
831 if(Sample_LLH->IsZombie())
833 MACH3LOG_WARN(
"Couldn't find Sample_LLH, it looks like LLH scan wasn't run, will do this now");
839 const int npars = cov->GetNumParams();
840 std::vector<double> StepScale(npars);
841 for (
int i = 0; i < npars; ++i)
843 std::string name = cov->GetParFancyName(i);
845 StepScale[i] = cov->GetIndivStepScale(i);
846 TH1D* LLHScan = Sample_LLH->Get<TH1D>((name+
"_sam").c_str());
847 if(LLHScan ==
nullptr)
849 MACH3LOG_WARN(
"Couldn't find LLH scan, for {}, skipping", name);
852 const double LLH_val = std::max(LLHScan->GetBinContent(1), LLHScan->GetBinContent(LLHScan->GetNbinsX()));
854 if(LLH_val < 0.001)
continue;
859 const double Var = 1.;
860 const double approxSigma = TMath::Abs(Var)/std::sqrt(LLH_val);
863 const double NewStepScale = approxSigma * 2.38/std::sqrt(npars);
864 StepScale[i] = NewStepScale;
868 cov->SetIndivStepScale(StepScale);
869 cov->SaveUpdatedMatrixConfig();
882 TDirectory *Sample_2DLLH =
outputFile->mkdir(
"Sample_2DLLH");
883 auto SkipVector = GetFromManager<std::vector<std::string>>(
fitMan->
raw()[
"LLHScan"][
"LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
886 const int n_points = GetFromManager<int>(
fitMan->
raw()[
"LLHScan"][
"2DLLHScanPoints"], 20, __FILE__ , __LINE__);
888 const int countwidth = int(
double(n_points)/
double(5));
890 std::map<std::string, std::vector<double>> scanRanges;
893 const double nSigma = GetFromManager<int>(
fitMan->
raw()[
"LLHScan"][
"LLHScanSigma"], 1., __FILE__, __LINE__);
900 int npars = cov->GetNumParams();
901 bool IsPCA = cov->IsPCA();
902 if (IsPCA) npars = cov->GetNParameters();
904 for (
int i = 0; i < npars; ++i)
906 std::string name_x = cov->GetParFancyName(i);
907 if (IsPCA) name_x +=
"_PCA";
910 double prior_x = cov->GetParInit(i);
911 if (IsPCA) prior_x = cov->GetPCAHandler()->GetParCurrPCA(i);
915 double lower_x = prior_x - nSigma*cov->GetDiagonalError(i);
916 double upper_x = prior_x + nSigma*cov->GetDiagonalError(i);
919 lower_x = prior_x - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
920 upper_x = prior_x + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
921 MACH3LOG_INFO(
"eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
930 auto it = scanRanges.find(name_x);
931 if (it != scanRanges.end() && it->second.size() == 2) {
932 lower_x = it->second[0];
933 upper_x = it->second[1];
934 MACH3LOG_INFO(
"Found matching param name for setting specified range for {}", name_x);
935 MACH3LOG_INFO(
"Range for {} = [{:.2f}, {:.2f}]", name_x, lower_x, upper_x);
940 lower_x = std::max(lower_x, cov->GetLowerBound(i));
941 upper_x = std::min(upper_x, cov->GetUpperBound(i));
945 for (
int j = 0; j < i; ++j)
947 std::string name_y = cov->GetParFancyName(j);
948 if (IsPCA) name_y +=
"_PCA";
953 double prior_y = cov->GetParInit(j);
954 if (IsPCA) prior_y = cov->GetPCAHandler()->GetParCurrPCA(j);
957 double lower_y = prior_y - nSigma*cov->GetDiagonalError(j);
958 double upper_y = prior_y + nSigma*cov->GetDiagonalError(j);
961 lower_y = prior_y - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
962 upper_y = prior_y + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
963 MACH3LOG_INFO(
"eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(j));
972 auto it = scanRanges.find(name_y);
973 if (it != scanRanges.end() && it->second.size() == 2) {
974 lower_y = it->second[0];
975 upper_y = it->second[1];
976 MACH3LOG_INFO(
"Found matching param name for setting specified range for {}", name_y);
977 MACH3LOG_INFO(
"Range for {} = [{:.2f}, {:.2f}]", name_y, lower_y, upper_y);
982 lower_y = std::max(lower_y, cov->GetLowerBound(j));
983 upper_y = std::min(upper_y, cov->GetUpperBound(j));
984 MACH3LOG_INFO(
"Scanning X {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_x, n_points, lower_x, upper_x, prior_x);
985 MACH3LOG_INFO(
"Scanning Y {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_y, n_points, lower_y, upper_y, prior_y);
987 auto hScanSam = std::make_unique<TH2D>((name_x +
"_" + name_y +
"_sam").c_str(), (name_x +
"_" + name_y +
"_sam").c_str(),
988 n_points, lower_x, upper_x, n_points, lower_y, upper_y);
989 hScanSam->SetDirectory(
nullptr);
990 hScanSam->GetXaxis()->SetTitle(name_x.c_str());
991 hScanSam->GetYaxis()->SetTitle(name_y.c_str());
992 hScanSam->GetZaxis()->SetTitle(
"2LLH_sam");
995 for (
int x = 0; x < n_points; ++x)
997 if (x % countwidth == 0)
1000 for (
int y = 0; y < n_points; ++y)
1004 cov->GetPCAHandler()->SetParPropPCA(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
1005 cov->GetPCAHandler()->SetParPropPCA(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
1008 cov->SetParProp(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
1009 cov->SetParProp(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
1012 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs) {
1017 double samplellh = 0;
1018 for(
unsigned int ivs = 0; ivs <
samples.size(); ++ivs) {
1019 samplellh +=
samples[ivs]->GetLikelihood();
1021 hScanSam->SetBinContent(x+1, y+1, 2*samplellh);
1029 cov->GetPCAHandler()->SetParPropPCA(i, prior_x);
1030 cov->GetPCAHandler()->SetParPropPCA(j, prior_y);
1032 cov->SetParProp(i, prior_x);
1033 cov->SetParProp(j, prior_y);
1038 Sample_2DLLH->Write();
1039 delete Sample_2DLLH;
1051 constexpr int numVar = 5;
1053 constexpr int sigmaArray[numVar] = {-3, -1, 0, 1, 3};
1059 constexpr int nRelevantModes = 2;
1061 bool DoByMode = GetFromManager<int>(
fitMan->
raw()[
"General"][
"DoByMode"],
false, __FILE__ , __LINE__);
1064 bool PlotLLHperBin =
false;
1065 auto SkipVector = GetFromManager<std::vector<std::string>>(
fitMan->
raw()[
"LLHScan"][
"LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
1069 TMatrixDSym *Cov = cov->GetCovMatrix();
1073 MACH3LOG_ERROR(
"Using PCAed matrix not implemented within sigma var code, I am sorry :(");
1078 for (
int i = 0; i < cov->GetNumParams(); ++i)
1081 std::string name = cov->GetParFancyName(i);
1086 TDirectory* dirArryDial =
outputFile->mkdir(name.c_str());
1089 int SampleIterator = 0;
1091 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++ )
1093 for(
int k = 0; k <
samples[ivs]->GetNsamples(); k++ )
1095 std::string title = std::string(
samples[ivs]->GetPDF(k)->
GetName());
1097 dirArrySample[SampleIterator] = dirArryDial->mkdir(title.c_str());
1103 double init = cov->GetParInit(i);
1105 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x(numVar);
1106 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y(numVar);
1107 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x_norm(numVar);
1108 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y_norm(numVar);
1111 TH1D ****sigmaArray_mode_x =
nullptr;
1112 TH1D ****sigmaArray_mode_y =
nullptr;
1115 sigmaArray_mode_x =
new TH1D***[numVar]();
1116 sigmaArray_mode_y =
new TH1D***[numVar]();
1122 for (
int j = 0; j < numVar; ++j)
1125 double paramVal = cov->GetParInit(i)+sigmaArray[j]*std::sqrt((*Cov)(i,i));
1128 paramVal = std::max(cov->GetLowerBound(i), std::min(paramVal, cov->GetUpperBound(i)));
1131 cov->SetParProp(i, paramVal);
1133 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++) {
1150 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++ )
1152 for (
int k = 0; k <
samples[ivs]->GetNsamples(); ++k)
1155 std::ostringstream ss;
1157 std::string parVarTitle = name +
"_" + ss.str();
1159 auto currSamp = M3::Clone<TH2Poly>(
static_cast<TH2Poly*
>(
samples[ivs]->GetPDF(k)));
1161 std::string title_long = std::string(currSamp->GetName())+
"_"+parVarTitle;
1165 std::vector<double> xbins;
1166 std::vector<double> ybins;
1174 MaCh3Modes_t RelevantModes[nRelevantModes] = {
samples[ivs]->GetMaCh3Modes()->GetMode(
"CCQE"),
samples[ivs]->GetMaCh3Modes()->GetMode(
"2p2h")};
1176 sigmaArray_mode_x[j][SampleIterator] =
new TH1D*[nRelevantModes]();
1177 sigmaArray_mode_y[j][SampleIterator] =
new TH1D*[nRelevantModes]();
1179 std::string mode_title_long;
1181 for(
int ir = 0; ir < nRelevantModes; ir++)
1183 auto currSampMode = M3::Clone<TH2Poly>(
static_cast<TH2Poly*
>(
samples[ivs]->GetPDFMode(k, RelevantModes[ir])));
1185 mode_title_long = title_long +
"_" +
samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir]);
1186 currSampMode->SetNameTitle(mode_title_long.c_str(), mode_title_long.c_str());
1187 dirArrySample[SampleIterator]->cd();
1188 currSampMode->Write();
1190 sigmaArray_mode_x[j][SampleIterator][ir] =
PolyProjectionX(currSampMode.get(), (mode_title_long+
"_xProj").c_str(), xbins);
1191 sigmaArray_mode_x[j][SampleIterator][ir]->SetDirectory(
nullptr);
1192 sigmaArray_mode_y[j][SampleIterator][ir] =
PolyProjectionY(currSampMode.get(), (mode_title_long+
"_yProj").c_str(), ybins);
1193 sigmaArray_mode_y[j][SampleIterator][ir]->SetDirectory(
nullptr);
1200 auto currLLHSamp = M3::Clone<TH2Poly>(
static_cast<TH2Poly*
>(
samples[ivs]->GetPDF(k)));
1201 currLLHSamp->Reset(
"");
1202 currLLHSamp->Fill(0.0, 0.0, 0.0);
1204 TH2Poly* MCpdf =
static_cast<TH2Poly*
>(
samples[ivs]->GetPDF(k));
1205 TH2Poly* Datapdf =
static_cast<TH2Poly*
>(
samples[ivs]->GetData(k));
1206 TH2Poly* W2pdf =
samples[ivs]->GetW2(k);
1208 for(
int bin = 1; bin < currLLHSamp->GetNumberOfBins()+1; bin++)
1210 const double mc = MCpdf->GetBinContent(bin);
1211 const double dat = Datapdf->GetBinContent(bin);
1212 const double w2 = W2pdf->GetBinContent(bin);
1213 currLLHSamp->SetBinContent(bin,
samples[ivs]->GetTestStatLLH(dat, mc, w2));
1215 currLLHSamp->SetNameTitle((title_long+
"_LLH").c_str() ,(title_long+
"_LLH").c_str());
1216 dirArrySample[SampleIterator]->cd();
1217 currLLHSamp->Write();
1221 sigmaArray_x[j][SampleIterator] = std::unique_ptr<TH1D>(
PolyProjectionX(currSamp.get(), (title_long+
"_xProj").c_str(), xbins));
1222 sigmaArray_x[j][SampleIterator]->SetDirectory(
nullptr);
1223 sigmaArray_x[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetXaxis()->GetTitle());
1224 sigmaArray_y[j][SampleIterator] = std::unique_ptr<TH1D>(
PolyProjectionY(currSamp.get(), (title_long+
"_yProj").c_str(), ybins));
1225 sigmaArray_y[j][SampleIterator]->SetDirectory(
nullptr);
1226 sigmaArray_y[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetYaxis()->GetTitle());
1228 sigmaArray_x_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_x[j][SampleIterator].get());
1229 sigmaArray_x_norm[j][SampleIterator]->Scale(1.,
"width");
1230 sigmaArray_y_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_y[j][SampleIterator].get());
1231 sigmaArray_y_norm[j][SampleIterator]->SetDirectory(
nullptr);
1232 sigmaArray_y_norm[j][SampleIterator]->Scale(1.,
"width");
1234 currSamp->SetNameTitle(title_long.c_str(), title_long.c_str());
1235 dirArrySample[k]->cd();
1238 sigmaArray_x[j][k]->Write();
1239 sigmaArray_y[j][k]->Write();
1246 cov->SetParProp(i, init);
1250 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++ )
1252 for (
int k = 0; k <
samples[ivs]->GetNsamples(); ++k)
1254 std::string title = std::string(
samples[ivs]->GetPDF(k)->
GetName()) +
"_" + name;
1255 auto var_x =
MakeAsymGraph(sigmaArray_x[1][SampleIterator].get(), sigmaArray_x[2][SampleIterator].get(), sigmaArray_x[3][SampleIterator].get(), (title+
"_X").c_str());
1256 auto var_y =
MakeAsymGraph(sigmaArray_y[1][SampleIterator].get(), sigmaArray_y[2][SampleIterator].get(), sigmaArray_y[3][SampleIterator].get(), (title+
"_Y").c_str());
1258 auto var_x_norm =
MakeAsymGraph(sigmaArray_x_norm[1][SampleIterator].get(), sigmaArray_x_norm[2][SampleIterator].get(), sigmaArray_x_norm[3][SampleIterator].get(), (title+
"_X_norm").c_str());
1259 auto var_y_norm =
MakeAsymGraph(sigmaArray_y_norm[1][SampleIterator].get(), sigmaArray_y_norm[2][SampleIterator].get(), sigmaArray_y_norm[3][SampleIterator].get(), (title+
"_Y_norm").c_str());
1261 dirArrySample[SampleIterator]->cd();
1264 var_x_norm->Write();
1265 var_y_norm->Write();
1272 MaCh3Modes_t RelevantModes[nRelevantModes] = {
samples[ivs]->GetMaCh3Modes()->GetMode(
"CCQE"),
samples[ivs]->GetMaCh3Modes()->GetMode(
"2p2h")};
1274 for(
int ir = 0; ir < nRelevantModes;ir++)
1276 auto var_mode_x =
MakeAsymGraph(sigmaArray_mode_x[1][SampleIterator][ir], sigmaArray_mode_x[2][SampleIterator][ir], sigmaArray_mode_x[3][SampleIterator][ir], (title+
"_"+
samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+
"_X").c_str());
1277 auto var_mode_y =
MakeAsymGraph(sigmaArray_mode_y[1][SampleIterator][ir], sigmaArray_mode_y[2][SampleIterator][ir], sigmaArray_mode_y[3][SampleIterator][ir], (title+
"_"+
samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+
"_Y").c_str());
1279 dirArrySample[SampleIterator]->cd();
1280 var_mode_x->Write();
1281 var_mode_y->Write();
1289 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++ )
1291 for (
int k = 0; k <
samples[ivs]->GetNsamples(); ++k)
1293 dirArrySample[SampleIterator]->Close();
1294 delete dirArrySample[SampleIterator];
1299 dirArryDial->Close();
1304 for (
int j = 0; j < numVar; ++j)
1307 for(
unsigned int ivs = 0; ivs <
samples.size(); ivs++ )
1309 for (
int k = 0; k <
samples[ivs]->GetNsamples(); ++k)
1311 for(
int ir = 0; ir < nRelevantModes;ir++)
1313 delete sigmaArray_mode_x[j][SampleIterator][ir];
1314 delete sigmaArray_mode_y[j][SampleIterator][ir];
1316 delete[] sigmaArray_mode_x[j][SampleIterator];
1317 delete[] sigmaArray_mode_y[j][SampleIterator];
1322 delete[] sigmaArray_mode_x;
1323 delete[] sigmaArray_mode_y;
#define _MaCh3_Safe_Include_Start_
KS: Avoiding warning checking for headers.
#define _MaCh3_Safe_Include_End_
KS: Restore warning checking after including external headers.
std::vector< TString > BranchNames
MaCh3Plotting::PlottingManager * man
TH1D * PolyProjectionX(TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
WP: Poly Projectors.
std::unique_ptr< TGraphAsymmErrors > MakeAsymGraph(TH1 *sigmaArrayLeft, TH1 *sigmaArrayCentr, TH1 *sigmaArrayRight, const std::string &title)
Used by sigma variation, check how 1 sigma changes spectra.
TH1D * PolyProjectionY(TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
WP: Poly Projectors.
#define MACH3LOG_CRITICAL
int MaCh3Modes_t
Enumerator of MaCh3Mode.
TMacro YAMLtoTMacro(const YAML::Node &yaml_node, const std::string &name)
Convert a YAML node to a ROOT TMacro object.
YAML::Node TMacroToYAML(const TMacro ¯o)
KS: Convert a ROOT TMacro object to a YAML node.
bool compareYAMLNodes(const YAML::Node &node1, const YAML::Node &node2)
Compare if yaml nodes are identical.
void RunLLHScan()
Perform a 1D likelihood scan.
void AddSystObj(ParameterHandlerBase *cov)
This function adds a Covariance object to the analysis framework. The Covariance object will be utili...
void AddSampleHandler(SampleHandlerBase *sample)
This function adds a sample PDF object to the analysis framework. The sample PDF object will be utili...
std::unique_ptr< TRandom3 > random
Random number.
std::vector< SampleHandlerBase * > samples
Sample holder.
int stepStart
step start if restarting
bool CheckSkipParameter(const std::vector< std::string > &SkipVector, const std::string &ParamName) const
KS: Check whether we want to skip parameter using skip vector.
void ProcessMCMC()
Process MCMC output.
int accCount
counts accepted steps
virtual std::string GetName() const
Get name of class.
bool OutputPrepared
Checks if output prepared not repeat some operations.
void SaveOutput()
Save output and close files.
TFile * outputFile
Output.
void SaveSettings()
Save the settings that the MCMC was run with.
manager * fitMan
The manager.
unsigned int step
current state
void PrepareOutput()
Prepare the output file.
bool SettingsSaved
Checks if setting saved not repeat some operations.
double accProb
current acceptance prob
void GetStepScaleBasedOnLLHScan()
LLH scan is good first estimate of step scale.
virtual void StartFromPreviousFit(const std::string &FitName)
Allow to start from previous fit/chain.
bool FileSaved
Checks if file saved not repeat some operations.
std::vector< double > sample_llh
store the llh breakdowns
void RunSigmaVar()
Perform a 2D and 1D sigma var for all samples.
double stepTime
Time of single step.
std::unique_ptr< TStopwatch > clock
tells global time how long fit took
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
TDirectory * CovFolder
Output cov folder.
TDirectory * SampleFolder
Output sample folder.
void DragRace(const int NLaps=100)
Calculates the required time for each sample or covariance object in a drag race simulation....
void Run2DLLHScan()
Perform a 2D likelihood scan.
unsigned int TotalNSamples
Total number of samples used.
double logLCurr
current likelihood
std::vector< double > syst_llh
systematic llh breakdowns
int auto_save
auto save every N steps
FitterBase(manager *const fitMan)
Constructor.
bool fTestLikelihood
Necessary for some fitting algorithms like PSO.
virtual ~FitterBase()
Destructor for the FitterBase class.
TTree * outTree
Output tree with posteriors.
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
bool GetScaneRange(std::map< std::string, std::vector< double > > &scanRanges)
YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D....
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
void Initialise()
Scan chain, what parameters we have and load information from covariance matrices.
const std::vector< TString > & GetBranchNames() const
Get the vector of branch names from root file.
void GetPostfit(TVectorD *&Central, TVectorD *&Errors, TVectorD *&Central_Gauss, TVectorD *&Errors_Gauss, TVectorD *&Peaks)
Get the post-fit results (arithmetic and Gaussian)
void DrawCovariance()
Draw the post-fit covariances.
void DrawPostfit()
Draw the post-fit comparisons.
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
Custom exception class for MaCh3 errors.
Base class responsible for handling of systematic error parameters. Capable of using PCA or using ada...
Class responsible for handling implementation of samples used in analysis, reweighting and returning ...
virtual void SaveAdditionalInfo(TDirectory *Dir)
Store additional info in a chan.
The manager class is responsible for managing configurations and settings.
void SaveSettings(TFile *const OutputFile)
Add manager useful information's to TFile, in most cases to Fitter.
YAML::Node const & raw()
Return config.
int GetNumParams() const
Get total number of parameters.
TH2D * GetCorrelationMatrix()
KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plott...
TMatrixDSym * GetCovMatrix() const
Return covariance matrix.
std::string GetName() const
Get name of covariance.
double GetParInit(const int i) const
Get prior parameter value.
YAML::Node GetConfig() const
Getter to return a copy of the YAML node.
virtual std::string GetTitle() const
virtual M3::int_t GetNsamples()
static constexpr const double _BAD_DOUBLE_
Default value used for double initialisation.
TFile * Open(const std::string &Name, const std::string &Type, const std::string &File, const int Line)
Opens a ROOT file with the given name and mode.
static constexpr const double _LARGE_LOGL_
Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such s...
int GetNThreads()
number of threads which we need for example for TRandom3
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
int getValue(const std::string &Type)
CW: Get info like RAM.