From b13888901938390a336dc0d40476177214bf8e65 Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Tue, 19 May 2026 09:49:43 +0000 Subject: [PATCH 1/4] Implement posterior to Chi2 conversion --- Diagnostics/ProcessMCMC.cpp | 1 + Fitters/MCMCProcessor.cpp | 32 ++++++++++++++++++++++++++++++++ Fitters/MCMCProcessor.h | 3 +++ Fitters/StatisticalUtils.cpp | 24 ++++++++++++++++++++++++ Fitters/StatisticalUtils.h | 11 +++++++++++ 5 files changed, 71 insertions(+) diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 1e59671b0..2ab45670f 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -205,6 +205,7 @@ void ProcessMCMC(const std::string& inputFile) //KS: When creating covariance matrix longest time is spend on caching every step, since we already cached we can run some fancy covariance stability diagnostic if(GetFromManager(Settings["DiagnoseCovarianceMatrix"], false)) DiagnoseCovarianceMatrix(Processor.get(), inputFile); } + Processor->ProduceChi2(GetFromManager(Settings["Chi2Group"], "Osc")); if(GetFromManager(Settings["JarlskogAnalysis"], true)) Processor->PerformJarlskogAnalysis(); if(GetFromManager(Settings["ProducePMNSElements"], true)) Processor->ProducePMNSElements(); if(GetFromManager(Settings["ProduceUnitarityTriangles"], true)) Processor->ProduceUnitarityTriangles(); diff --git a/Fitters/MCMCProcessor.cpp b/Fitters/MCMCProcessor.cpp index 1a225a94f..ea1461181 100644 --- a/Fitters/MCMCProcessor.cpp +++ b/Fitters/MCMCProcessor.cpp @@ -1764,6 +1764,38 @@ void MCMCProcessor::DrawCorrelations1D() { gStyle->SetOptTitle(OptTitle); } + +// ********************* +// Convert posterior likelihood to Delta Chi2 used for comparison with frequentists fitter +void MCMCProcessor::ProduceChi2(const std::string& GroupName) const { +// ********************* + if(GroupName == "") return; + MACH3LOG_INFO("Starting {}", __func__); + TDirectory* Chi2Folder = OutputFile->mkdir("DeltaChi2"); + + Chi2Folder->cd(); + for (int iPar = 0; iPar < nDraw; iPar++) + { + std::string GroupNameCurr; + if(ParamType[iPar] == kXSecPar){ + const int InternalNumeration = iPar - ParamTypeStartPos[kXSecPar]; + GroupNameCurr = ParameterGroup[InternalNumeration]; + } else { + GroupNameCurr = "Other"; // Use Other for all legacy params + } + if (IamVaried[iPar] == false) continue; + if (GroupName != "All" && GroupNameCurr != GroupName) continue; + + auto Chi2 = GetDeltaChi2(hpost[iPar]); + RemoveFitter(Chi2.get(), "Gauss"); + + Chi2->Write(); + } + Chi2Folder->Close(); + delete Chi2Folder; + OutputFile->cd(); +} + // ********************* // Make fancy Credible Intervals plots void MCMCProcessor::MakeCredibleRegions(const std::vector& CredibleRegions, diff --git a/Fitters/MCMCProcessor.h b/Fitters/MCMCProcessor.h index eaef3284f..ef2bcef85 100644 --- a/Fitters/MCMCProcessor.h +++ b/Fitters/MCMCProcessor.h @@ -179,6 +179,9 @@ class MCMCProcessor { const std::string& Title, const double EvaluationPoint) const ; + /// @brief Convert posterior likelihood to Delta Chi2 used for comparison with frequentists fitter + void ProduceChi2(const std::string& GroupName) const; + /// @brief Reweight Prior by giving new central value and new error /// @param Names Parameter names for which we do reweighting /// @param NewCentral New central value for which we reweight diff --git a/Fitters/StatisticalUtils.cpp b/Fitters/StatisticalUtils.cpp index a24b97a29..e3ff03d63 100644 --- a/Fitters/StatisticalUtils.cpp +++ b/Fitters/StatisticalUtils.cpp @@ -672,7 +672,31 @@ void Get2DBayesianpValue(TH2D *Histogram) { TempCanvas->Write(); } +// **************** +// Converts posterior likelihood to dchi2 +std::unique_ptr GetDeltaChi2(TH1D* posterior_probability_hist) { +// **************** + auto delta_chi2 = M3::Clone(posterior_probability_hist); + TString title = delta_chi2->GetTitle(); + title+=";#Delta#chi^{2}"; + + delta_chi2->GetYaxis()->SetTitle("#Delta#chi^{2}"); + + int max_bin = delta_chi2->GetMaximumBin(); + double max_content = delta_chi2->GetBinContent(max_bin); + double NewMaximum = M3::_BAD_DOUBLE_; + for(int iBin = 1; iBin < delta_chi2->GetNbinsX()+1; iBin++) { + double bin_content = delta_chi2->GetBinContent(iBin); + if(bin_content == 0) bin_content = 1.0 ; + + double chi2_likelihood = -2*std::log(bin_content/max_content); + delta_chi2->SetBinContent(iBin, chi2_likelihood); + NewMaximum = std::max(NewMaximum, chi2_likelihood); + } + delta_chi2->SetMaximum(NewMaximum*1.1); + return delta_chi2; +} // **************** void PassErrorToRatioPlot(TH1D* RatioHist, TH1D* Hist1, TH1D* DataHist) { diff --git a/Fitters/StatisticalUtils.h b/Fitters/StatisticalUtils.h index c14a1f166..3efadadf3 100644 --- a/Fitters/StatisticalUtils.h +++ b/Fitters/StatisticalUtils.h @@ -173,6 +173,17 @@ double ComputeKLDivergence(TH2Poly* DataPoly, TH2Poly* PolyMC); /// @return The combined p-value, representing the overall significance. double FisherCombinedPValue(const std::vector& pvalues); +/// @brief Convert a posterior probability histogram into a \f$\Delta\chi^2\f$ distribution. Using +/// the likelihood-ratio definition: +/// +/// \f[ +/// \Delta\chi^2 = -2 \ln\left(\frac{L}{L_{\max}}\right) +/// \f] +/// +/// @param posterior_probability_hist Pointer to a TH1D histogram containing posterior probabilities +/// @note based on CompareMaCh3PThetaDeltaChi2.C +std::unique_ptr GetDeltaChi2(TH1D* posterior_probability_hist); + /// @brief Thin MCMC Chain, to save space and maintain low autocorrelations. /// /// @param FilePath Path to MCMC chain you want to thin From 92414b341d4b81cd8e9ec5bfa02713a3ce8f8cbf Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Tue, 19 May 2026 15:52:17 +0000 Subject: [PATCH 2/4] adress comment by Henry --- Fitters/StatisticalUtils.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Fitters/StatisticalUtils.cpp b/Fitters/StatisticalUtils.cpp index e3ff03d63..4619b0218 100644 --- a/Fitters/StatisticalUtils.cpp +++ b/Fitters/StatisticalUtils.cpp @@ -677,13 +677,15 @@ void Get2DBayesianpValue(TH2D *Histogram) { std::unique_ptr GetDeltaChi2(TH1D* posterior_probability_hist) { // **************** auto delta_chi2 = M3::Clone(posterior_probability_hist); - TString title = delta_chi2->GetTitle(); - title+=";#Delta#chi^{2}"; - delta_chi2->GetYaxis()->SetTitle("#Delta#chi^{2}"); int max_bin = delta_chi2->GetMaximumBin(); double max_content = delta_chi2->GetBinContent(max_bin); + if (max_content == 0) { + MACH3LOG_ERROR("Histogram {}, has larges bin with 0", delta_chi2->GetTitle()); + MACH3LOG_ERROR("This suggest you skewed binning for posterior probability or something else"); + throw MaCh3Exception(__FILE__, __LINE__); + } double NewMaximum = M3::_BAD_DOUBLE_; for(int iBin = 1; iBin < delta_chi2->GetNbinsX()+1; iBin++) { From b4e3440fd0aa83609e0f89d84d46ad9e0df43ff3 Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Tue, 19 May 2026 16:15:59 +0000 Subject: [PATCH 3/4] more tweaks after Henry --- Fitters/MCMCProcessor.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/Fitters/MCMCProcessor.cpp b/Fitters/MCMCProcessor.cpp index ea1461181..2eeafa43d 100644 --- a/Fitters/MCMCProcessor.cpp +++ b/Fitters/MCMCProcessor.cpp @@ -344,7 +344,7 @@ void MCMCProcessor::MakePostfit(const std::mapGetMaximum() == hpost[i]->Integral()*DrawRange) { MACH3LOG_WARN("Found fixed parameter: {} ({}), moving on", Title, i); - IamVaried[i] = false; + ParamVaried[i] = false; //KS:Set mean and error to prior for fixed parameters, it looks much better when fixed parameter has mean on prior rather than on 0 with 0 error. (*Means_HPD)(i) = Prior; (*Errors_HPD)(i) = PriorError; @@ -360,7 +360,7 @@ void MCMCProcessor::MakePostfit(const std::mapDraw(); @@ -746,7 +746,7 @@ void MCMCProcessor::MakeCredibleIntervals(const std::vector& CredibleInt for (int i = 0; i < nDraw; ++i) { - if(!IamVaried[i]) continue; + if(!ParamVaried[i]) continue; // Now make the TLine for the Asimov TString Title = ""; @@ -854,7 +854,7 @@ void MCMCProcessor::MakeViolin() { for (int x = 0; x < nDraw; ++x) { //KS: Consider another treatment for fixed params - //if (IamVaried[x] == false) continue; + //if (ParamVaried[x] == false) continue; for (int k = 0; k < nEntries; ++k) { //KS: Burn in cut @@ -985,7 +985,7 @@ void MCMCProcessor::MakeCovariance() { if (j == i) continue; // If this parameter isn't varied - if (IamVaried[j] == false) { + if (ParamVaried[j] == false) { (*Covariance)(i,j) = 0.0; (*Covariance)(j,i) = (*Covariance)(i,j); (*Correlation)(i,j) = 0.0; @@ -1228,7 +1228,7 @@ void MCMCProcessor::MakeCovariance_MP(const bool Mute) { if (j == i) continue; // If this parameter isn't varied - if (IamVaried[j] == false) { + if (ParamVaried[j] == false) { (*Covariance)(i,j) = 0.0; (*Covariance)(j,i) = (*Covariance)(i,j); (*Correlation)(i,j) = 0.0; @@ -1271,7 +1271,7 @@ void MCMCProcessor::MakeCovariance_MP(const bool Mute) { { // Skip the diagonal elements which we've already done above if (j == i) continue; - if (IamVaried[j] == false) continue; + if (ParamVaried[j] == false) continue; if(ParamType[i] == kXSecPar && ParamType[j] == kXSecPar) { @@ -1707,7 +1707,7 @@ void MCMCProcessor::DrawCorrelations1D() { for(int i = 0; i < nDraw; i++) { - if (IamVaried[i] == false) continue; + if (ParamVaried[i] == false) continue; Corr1DHist[i][0]->GetXaxis()->LabelsOption("v"); Corr1DHist[i][0]->SetMaximum(+1.); @@ -1783,7 +1783,7 @@ void MCMCProcessor::ProduceChi2(const std::string& GroupName) const { } else { GroupNameCurr = "Other"; // Use Other for all legacy params } - if (IamVaried[iPar] == false) continue; + if (ParamVaried[iPar] == false) continue; if (GroupName != "All" && GroupNameCurr != GroupName) continue; auto Chi2 = GetDeltaChi2(hpost[iPar]); @@ -1854,7 +1854,7 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector& CredibleRegio { // Skip the diagonal elements which we've already done above if (j == i) continue; - if (IamVaried[j] == false) continue; + if (ParamVaried[j] == false) continue; auto legend = std::make_unique(0.20, 0.7, 0.4, 0.92); legend->SetTextColor(kRed); @@ -2312,7 +2312,7 @@ void MCMCProcessor::ScanInput() { // Check order of parameter types ScanParameterOrder(); - IamVaried.resize(nDraw, true); + ParamVaried.resize(nDraw, true); // Print useful Info PrintInfo(); From 13cea0ec632734fe736ccbc4542043c6efa52832 Mon Sep 17 00:00:00 2001 From: Kamil Skwarczynski Date: Tue, 19 May 2026 16:25:18 +0000 Subject: [PATCH 4/4] missed header --- Fitters/MCMCProcessor.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Fitters/MCMCProcessor.h b/Fitters/MCMCProcessor.h index ef2bcef85..f97654471 100644 --- a/Fitters/MCMCProcessor.h +++ b/Fitters/MCMCProcessor.h @@ -463,7 +463,7 @@ class MCMCProcessor { std::vector ExcludedGroups; /// Is the ith parameter varied - std::vector IamVaried; + std::vector ParamVaried; /// Name of parameters which we are going to analyse std::vector> ParamNames; /// Parameters central values which we are going to analyse