MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
MCMCProcessor Class Reference

Class responsible for processing MCMC chains, performing diagnostics, generating plots, and managing Bayesian analysis. More...

#include <Fitters/MCMCProcessor.h>

Inheritance diagram for MCMCProcessor:
[legend]
Collaboration diagram for MCMCProcessor:
[legend]

Public Member Functions

 MCMCProcessor (const std::string &InputFile)
 Constructs an MCMCProcessor object with the specified input file and options.
 
virtual ~MCMCProcessor ()
 Destroys the MCMCProcessor object.
 
void Initialise ()
 Scan chain, what parameters we have and load information from covariance matrices.
 
void MakePostfit ()
 Make 1D projection for each parameter and prepare structure.
 
void MakeCovariance ()
 Calculate covariance by making 2D projection of each combination of parameters.
 
void CacheSteps ()
 KS:By caching each step we use multithreading.
 
void MakeCovariance_MP (const bool Mute=false)
 Calculate covariance by making 2D projection of each combination of parameters using multithreading.
 
void MakeSubOptimality (const int NIntervals=10)
 Make and Draw SubOptimality [24].
 
void ResetHistograms ()
 Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.
 
void DrawPostfit ()
 Draw the post-fit comparisons.
 
void MakeViolin ()
 Make and Draw Violin.
 
void MakeCredibleIntervals (const std::vector< double > &CredibleIntervals={0.99, 0.90, 0.68 }, const std::vector< Color_t > &CredibleIntervalsColours={kCyan+4, kCyan-2, kCyan-10}, const bool CredibleInSigmas=false)
 Make and Draw Credible intervals.
 
void DrawCovariance ()
 Draw the post-fit covariances.
 
void MakeCredibleRegions (const std::vector< double > &CredibleRegions={0.99, 0.90, 0.68}, const std::vector< Style_t > &CredibleRegionStyle={kDashed, kSolid, kDotted}, const std::vector< Color_t > &CredibleRegionColor={kGreen-3, kGreen-10, kGreen}, const bool CredibleInSigmas=false, const bool Draw2DPosterior=true, const bool DrawBestFit=true)
 Make and Draw Credible Regions.
 
void MakeTrianglePlot (const std::vector< std::string > &ParNames, const std::vector< double > &CredibleIntervals={0.99, 0.90, 0.68 }, const std::vector< Color_t > &CredibleIntervalsColours={kCyan+4, kCyan-2, kCyan-10}, const std::vector< double > &CredibleRegions={0.99, 0.90, 0.68}, const std::vector< Style_t > &CredibleRegionStyle={kDashed, kSolid, kDotted}, const std::vector< Color_t > &CredibleRegionColor={kGreen-3, kGreen-10, kGreen}, const bool CredibleInSigmas=false)
 Make fancy triangle plot for selected parameters.
 
void CheckCredibleIntervalsOrder (const std::vector< double > &CredibleIntervals, const std::vector< Color_t > &CredibleIntervalsColours)
 Checks the order and size consistency of the CredibleIntervals and CredibleIntervalsColours vectors.
 
void CheckCredibleRegionsOrder (const std::vector< double > &CredibleRegions, const std::vector< Style_t > &CredibleRegionStyle, const std::vector< Color_t > &CredibleRegionColor)
 Checks the order and size consistency of the CredibleRegions, CredibleRegionStyle, and CredibleRegionColor vectors.
 
void GetPolarPlot (const std::vector< std::string > &ParNames)
 Make funny polar plot.
 
void GetBayesFactor (const std::vector< std::string > &ParName, const std::vector< std::vector< double > > &Model1Bounds, const std::vector< std::vector< double > > &Model2Bounds, const std::vector< std::vector< std::string > > &ModelNames)
 Calculate Bayes factor for vector of params, and model boundaries.
 
void GetSavageDickey (const std::vector< std::string > &ParName, const std::vector< double > &EvaluationPoint, const std::vector< std::vector< double > > &Bounds)
 Calculate Bayes factor for point like hypothesis using SavageDickey.
 
void ReweightPrior (const std::vector< std::string > &Names, const std::vector< double > &NewCentral, const std::vector< double > &NewError)
 Reweight Prior by giving new central value and new error.
 
void ParameterEvolution (const std::vector< std::string > &Names, const std::vector< int > &NIntervals)
 Make .gif of parameter evolution.
 
void ThinMCMC (const int ThinningCut)
 Thin MCMC Chain, to save space and maintain low autocorrelations.
 
void DiagMCMC ()
 KS: Perform MCMC diagnostic including Autocorrelation, Trace etc.
 
int GetNParams ()
 Get total number of used parameters.
 
int GetNXSec ()
 
int GetNND ()
 
int GetNFD ()
 
int GetGroup (const std::string &name) const
 Number of params from a given group, for example flux.
 
TH1D * GetHpost (const int i)
 Get 1D posterior for a given parameter.
 
TH2D * GetHpost2D (const int i, const int j)
 Get 2D posterior for a given parameter combination.
 
TH2D * GetViolin ()
 Get Violin plot for all parameters with posterior values.
 
TH2D * GetViolinPrior ()
 Get Violin plot for all parameters with prior values.
 
std::vector< std::string > GetXSecCov () const
 
std::string GetNDCov () const
 
std::string GetFDCov () const
 
void GetPostfit (TVectorD *&Central, TVectorD *&Errors, TVectorD *&Central_Gauss, TVectorD *&Errors_Gauss, TVectorD *&Peaks)
 Get the post-fit results (arithmetic and Gaussian)
 
void GetCovariance (TMatrixDSym *&Cov, TMatrixDSym *&Corr)
 Get the post-fit covariances and correlations.
 
void GetPostfit_Ind (TVectorD *&Central, TVectorD *&Errors, TVectorD *&Peaks, ParameterEnum kParam)
 Or the individual post-fits.
 
const std::vector< TString > & GetBranchNames () const
 Get the vector of branch names from root file.
 
void GetNthParameter (const int param, double &Prior, double &PriorError, TString &Title) const
 Get properties of parameter by passing it number.
 
int GetParamIndexFromName (const std::string &Name)
 Get parameter number based on name.
 
Long64_t GetnEntries ()
 Get Number of entries that Chain has, for merged chains will not be the same Nsteps.
 
Long64_t GetnSteps ()
 Get Number of Steps that Chain has, for merged chains will not be the same nEntries.
 
void SetStepCut (const std::string &Cuts)
 Set the step cutting by string.
 
void SetStepCut (const int Cuts)
 Set the step cutting by int.
 
void SetPlotRelativeToPrior (const bool PlotOrNot)
 You can set relative to prior or relative to generated. It is advised to use relate to prior.
 
void SetPrintToPDF (const bool PlotOrNot)
 
void SetPlotErrorForFlatPrior (const bool PlotOrNot)
 Set whether you want to plot error for parameters which have flat prior.
 
void SetPlotBinValue (const bool PlotOrNot)
 
void SetFancyNames (const bool PlotOrNot)
 
void SetSmoothing (const bool PlotOrNot)
 Set whether want to use smoothing for histograms using ROOT algorithm.
 
void SetPost2DPlotThreshold (const double Threshold)
 Code will only plot 2D posteriors if Correlation are larger than defined threshold.
 
void SetUseFFTAutoCorrelation (const bool useFFT)
 Toggle using the FFT-based autocorrelation calculator.
 
void SetExcludedTypes (std::vector< std::string > Name)
 Setter related what parameters we want to exclude from analysis, for example if cross-section parameters look like xsec_, then passing "xsec_" will.
 
void SetExcludedNames (std::vector< std::string > Name)
 
void SetnBatches (const int Batches)
 Set value of Nbatches used for batched mean, this need to be done earlier as batches are made when reading tree.
 
void SetnLags (const int nLags)
 
void SetOutputSuffix (const std::string Suffix)
 Sett output suffix, this way jobs using the same file will have different names.
 
void SetPosterior1DCut (const std::string Cut)
 Allow to set addtional cuts based on ROOT TBrowser cut, for to only affect one mass ordering.
 

Protected Member Functions

std::unique_ptr< TH1D > MakePrefit ()
 Prepare prefit histogram for parameter overlay plot.
 
void MakeOutputFile ()
 prepare output root file and canvas to which we will save EVERYTHING
 
void DrawCorrelations1D ()
 Draw 1D correlations which might be more helpful than looking at huge 2D Corr matrix.
 
void ReadInputCov ()
 CW: Read the input Covariance matrix entries. Get stuff like parameter input errors, names, and so on.
 
void ReadInputCovLegacy ()
 
void FindInputFiles ()
 Read the output MCMC file and find what inputs were used.
 
void FindInputFilesLegacy ()
 
void ReadModelFile ()
 Read the xsec file and get the input central values and errors.
 
virtual void LoadAdditionalInfo ()
 allow loading additional info for example used for oscillation parameters
 
void ReadNDFile ()
 Read the ND cov file and get the input central values and errors.
 
void ReadFDFile ()
 Read the FD cov file and get the input central values and errors.
 
void PrintInfo () const
 Print info like how many params have been loaded etc.
 
void ScanInput ()
 Scan Input etc.
 
void ScanParameterOrder ()
 Scan order of params from a different groups.
 
void SetupOutput ()
 Prepare all objects used for output.
 
void PrepareDiagMCMC ()
 CW: Prepare branches etc. for DiagMCMC.
 
void ParamTraces ()
 CW: Draw trace plots of the parameters i.e. parameter vs step.
 
void AutoCorrelation ()
 KS: Calculate autocorrelations supports both OpenMP and CUDA :)
 
void AutoCorrelation_FFT ()
 MJR: Autocorrelation function using FFT algorithm for extra speed.
 
void CalculateESS (const int nLags, const std::vector< std::vector< double > > &LagL)
 KS: calc Effective Sample Size.
 
void BatchedAnalysis ()
 Get the batched means variance estimation and variable indicating if number of batches is sensible [4] [25].
 
void BatchedMeans ()
 CW: Batched means, literally read from an array and chuck into TH1D.
 
void GewekeDiagnostic ()
 Geweke Diagnostic based on the methods described by Fang (2014) and Karlsbakk (2011). [8] [18].
 
void AcceptanceProbabilities ()
 Acceptance Probability.
 
void PowerSpectrumAnalysis ()
 RC: Perform spectral analysis of MCMC [7].
 
std::vector< double > GetMargins (const std::unique_ptr< TCanvas > &Canv) const
 Get TCanvas margins, to be able to reset them if particular function need different margins.
 
void SetMargins (std::unique_ptr< TCanvas > &Canv, const std::vector< double > &margins)
 Set TCanvas margins to specified values.
 
void SetTLineStyle (TLine *Line, const Color_t Colour, const Width_t Width, const ELineStyle Style) const
 Configures a TLine object with the specified style parameters.
 
void SetLegendStyle (TLegend *Legend, const double size) const
 Configures the style of a TLegend object.
 

Protected Attributes

std::string MCMCFile
 Name of MCMC file.
 
std::string OutputSuffix
 Output file suffix useful when running over same file with different settings.
 
std::vector< std::vector< std::string > > CovPos
 Covariance matrix name position.
 
std::vector< YAML::Node > CovConfig
 Covariance matrix config.
 
TChain * Chain
 Main chain storing all steps etc.
 
std::string StepCut
 BurnIn Cuts.
 
std::string Posterior1DCut
 Cut used when making 1D Posterior distribution.
 
int UpperCut
 KS: Used only for SubOptimality.
 
int BurnInCut
 Value of burn in cut.
 
int nBranches
 Number of branches in a TTree.
 
int nEntries
 KS: For merged chains number of entries will be different from nSteps.
 
int nSteps
 KS: For merged chains number of entries will be different from nSteps.
 
int nSamples
 Number of sample PDF objects.
 
int nSysts
 Number of covariance objects.
 
int nDraw
 Number of all parameters used in the analysis.
 
std::vector< TString > BranchNames
 
std::vector< std::string > ExcludedTypes
 
std::vector< std::string > ExcludedNames
 
std::vector< bool > IamVaried
 Is the ith parameter varied.
 
std::vector< std::vector< TString > > ParamNames
 Name of parameters which we are going to analyse.
 
std::vector< std::vector< double > > ParamCentral
 Parameters central values which we are going to analyse.
 
std::vector< std::vector< double > > ParamNom
 
std::vector< std::vector< double > > ParamErrors
 Uncertainty on a single parameter.
 
std::vector< std::vector< bool > > ParamFlat
 Whether Param has flat prior or not.
 
std::vector< int > nParam
 Number of parameters per type.
 
std::vector< ParameterEnum > ParamType
 Make an enum for which class this parameter belongs to so we don't have to keep string comparing.
 
std::vector< int > ParamTypeStartPos
 
std::vector< std::string > ParameterGroup
 
std::vector< TString > SampleName_v
 Vector of each systematic.
 
std::vector< TString > SystName_v
 Vector of each sample PDF object.
 
std::string OutputName
 Name of output files.
 
TString CanvasName
 Name of canvas which help to save to the sample pdf.
 
bool PlotFlatPrior
 Whether we plot flat prior or not, we usually provide error even for flat prior params.
 
bool plotRelativeToPrior
 Whether we plot relative to prior or nominal, in most cases is prior.
 
bool MadePostfit
 Sanity check if Postfit is already done to not make several times.
 
bool printToPDF
 Will plot all plot to PDF not only to root file.
 
bool FancyPlotNames
 Whether we want fancy plot names or not.
 
bool plotBinValue
 If true it will print value on each bin of covariance matrix.
 
bool ApplySmoothing
 Apply smoothing for 2D histos using root algorithm.
 
double Post2DPlotThreshold
 KS: Set Threshold when to plot 2D posterior as by default we get a LOT of plots.
 
bool useFFTAutoCorrelation
 MJR: Use FFT-based autocorrelation algorithm (save time & resources)?
 
std::vector< int > NDSamplesBins
 
std::vector< std::string > NDSamplesNames
 
std::unique_ptr< TF1 > Gauss
 Gaussian fitter.
 
TFile * OutputFile
 The output file.
 
std::unique_ptr< TCanvas > Posterior
 Fancy canvas used for our beautiful plots.
 
TVectorD * Central_Value
 Vector with central value for each parameter.
 
TVectorD * Means
 Vector with mean values using Arithmetic Mean.
 
TVectorD * Errors
 Vector with errors values using RMS.
 
TVectorD * Means_Gauss
 Vector with mean values using Gaussian fit.
 
TVectorD * Errors_Gauss
 Vector with error values using Gaussian fit.
 
TVectorD * Means_HPD
 Vector with mean values using Highest Posterior Density.
 
TVectorD * Errors_HPD
 Vector with error values using Highest Posterior Density.
 
TVectorD * Errors_HPD_Positive
 Vector with positive error (right hand side) values using Highest Posterior Density.
 
TVectorD * Errors_HPD_Negative
 Vector with negative error (left hand side) values using Highest Posterior Density.
 
TMatrixDSym * Covariance
 Posterior Covariance Matrix.
 
TMatrixDSym * Correlation
 Posterior Correlation Matrix.
 
std::vector< TH1D * > hpost
 Holds 1D Posterior Distributions.
 
std::vector< std::vector< TH2D * > > hpost2D
 Holds 2D Posterior Distributions.
 
std::unique_ptr< TH2D > hviolin
 Holds violin plot for all dials.
 
std::unique_ptr< TH2D > hviolin_prior
 Holds prior violin plot for all dials,.
 
double ** ParStep
 Array holding values for all parameters.
 
int * StepNumber
 Step number for step, important if chains were merged.
 
int nBins
 Number of bins.
 
double DrawRange
 Drawrange for SetMaximum.
 
bool CacheMCMC
 MCMC Chain has been cached.
 
bool doDiagMCMC
 Doing MCMC Diagnostic.
 
int nBatches
 Number of batches for Batched Mean.
 
int AutoCorrLag
 LagL used in AutoCorrelation.
 
double * ParamSums
 Total parameter sum for each param.
 
double ** BatchedAverages
 Values of batched average for every param and batch.
 
double ** SampleValues
 Holds the sample values.
 
double ** SystValues
 Holds the systs values.
 
double * AccProbValues
 Holds all accProb.
 
double * AccProbBatchedAverages
 Holds all accProb in batches.
 

Detailed Description

Class responsible for processing MCMC chains, performing diagnostics, generating plots, and managing Bayesian analysis.

This class provides utilities to handle MCMC output generated by mcmc::runMCMC. It is particularly useful for extracting values from previous MCMC runs and initiating new MCMC runs with those values. Inspired by nd280_utils/DrawComp.cpp.

See also
For more details and examples, visit the Wiki.
Author
Clarence Wret
Kamil Skwarczynski

Definition at line 65 of file MCMCProcessor.h.

Constructor & Destructor Documentation

◆ MCMCProcessor()

_MaCh3_Safe_Include_Start_ _MaCh3_Safe_Include_End_ MCMCProcessor::MCMCProcessor ( const std::string &  InputFile)

Constructs an MCMCProcessor object with the specified input file and options.

Parameters
InputFileThe path to the input file containing MCMC data.

Definition at line 17 of file MCMCProcessor.cpp.

17 :
18 Chain(nullptr), StepCut(""), MadePostfit(false) {
19// ****************************
20 MCMCFile = InputFile;
21
22 SetMaCh3LoggerFormat();
23 MaCh3Utils::MaCh3Welcome();
24 MACH3LOG_INFO("Making post-fit processor for: {}", MCMCFile);
25
26 ParStep = nullptr;
27 StepNumber = nullptr;
28
29 Posterior = nullptr;
30 hviolin = nullptr;
31 hviolin_prior = nullptr;
32
33 OutputFile = nullptr;
34
35 ParamSums = nullptr;
36 BatchedAverages = nullptr;
37 SampleValues = nullptr;
38 SystValues = nullptr;
39 AccProbValues = nullptr;
40 AccProbBatchedAverages = nullptr;
41
42 //KS:Hardcoded should be a way to get it via config or something
43 plotRelativeToPrior = false;
44 printToPDF = false;
45 plotBinValue = false;
46 PlotFlatPrior = true;
47 CacheMCMC = false;
48 ApplySmoothing = true;
49 FancyPlotNames = true;
50 doDiagMCMC = false;
51
52 // KS: ROOT can compile FFT code but it will crash during run time. Turn off FFT dynamically
53#ifdef MaCh3_FFT
54 useFFTAutoCorrelation = true;
55#else
56 useFFTAutoCorrelation = false;
57#endif
58 OutputSuffix = "_Process";
59 Post2DPlotThreshold = 1.e-5;
60
61 nDraw = 0;
62 nEntries = 0;
63 UpperCut = _UNDEF_;
64 nSteps = 0;
65 nBatches = 0;
66 AutoCorrLag = 0;
67 nSysts = 0;
68 nSamples = 0;
69
70 nBins = 70;
71 DrawRange = 1.5;
72
73 Posterior1DCut = "";
74 //KS:Those keep basic information for ParameterEnum
75 ParamNames.resize(kNParameterEnum);
76 ParamCentral.resize(kNParameterEnum);
77 ParamNom.resize(kNParameterEnum);
78 ParamErrors.resize(kNParameterEnum);
79 ParamFlat.resize(kNParameterEnum);
80 ParamTypeStartPos.resize(kNParameterEnum);
81 nParam.resize(kNParameterEnum);
82 CovPos.resize(kNParameterEnum);
83 CovConfig.resize(kNParameterEnum);
84
85 for(int i = 0; i < kNParameterEnum; i++)
86 {
87 ParamTypeStartPos[i] = 0;
88 nParam[i] = 0;
89 }
90 //Only if GPU is enabled
91 #ifdef MaCh3_CUDA
92 ParStep_cpu = nullptr;
93 NumeratorSum_cpu = nullptr;
94 ParamSums_cpu = nullptr;
95 DenomSum_cpu = nullptr;
96
97 ParStep_gpu = nullptr;
98 NumeratorSum_gpu = nullptr;
99 ParamSums_gpu = nullptr;
100 DenomSum_gpu = nullptr;
101 #endif
102}
_UNDEF_
#define _UNDEF_
Definition: MCMCProcessor.h:4
kNParameterEnum
@ kNParameterEnum
Definition: MCMCProcessor.h:57
MACH3LOG_INFO
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
SetMaCh3LoggerFormat
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:30
MCMCProcessor::nBatches
int nBatches
Number of batches for Batched Mean.
Definition: MCMCProcessor.h:513
MCMCProcessor::StepNumber
int * StepNumber
Step number for step, important if chains were merged.
Definition: MCMCProcessor.h:499
MCMCProcessor::ParStep
double ** ParStep
Array holding values for all parameters.
Definition: MCMCProcessor.h:497
MCMCProcessor::Post2DPlotThreshold
double Post2DPlotThreshold
KS: Set Threshold when to plot 2D posterior as by default we get a LOT of plots.
Definition: MCMCProcessor.h:446
MCMCProcessor::SampleValues
double ** SampleValues
Holds the sample values.
Definition: MCMCProcessor.h:523
MCMCProcessor::ParamCentral
std::vector< std::vector< double > > ParamCentral
Parameters central values which we are going to analyse.
Definition: MCMCProcessor.h:403
MCMCProcessor::ParamErrors
std::vector< std::vector< double > > ParamErrors
Uncertainty on a single parameter.
Definition: MCMCProcessor.h:406
MCMCProcessor::hviolin_prior
std::unique_ptr< TH2D > hviolin_prior
Holds prior violin plot for all dials,.
Definition: MCMCProcessor.h:494
MCMCProcessor::nSamples
int nSamples
Number of sample PDF objects.
Definition: MCMCProcessor.h:387
MCMCProcessor::nParam
std::vector< int > nParam
Number of parameters per type.
Definition: MCMCProcessor.h:410
MCMCProcessor::DrawRange
double DrawRange
Drawrange for SetMaximum.
Definition: MCMCProcessor.h:504
MCMCProcessor::ParamFlat
std::vector< std::vector< bool > > ParamFlat
Whether Param has flat prior or not.
Definition: MCMCProcessor.h:408
MCMCProcessor::CovConfig
std::vector< YAML::Node > CovConfig
Covariance matrix config.
Definition: MCMCProcessor.h:368
MCMCProcessor::ParamNom
std::vector< std::vector< double > > ParamNom
Definition: MCMCProcessor.h:404
MCMCProcessor::SystValues
double ** SystValues
Holds the systs values.
Definition: MCMCProcessor.h:525
MCMCProcessor::AccProbBatchedAverages
double * AccProbBatchedAverages
Holds all accProb in batches.
Definition: MCMCProcessor.h:530
MCMCProcessor::useFFTAutoCorrelation
bool useFFTAutoCorrelation
MJR: Use FFT-based autocorrelation algorithm (save time & resources)?
Definition: MCMCProcessor.h:448
MCMCProcessor::StepCut
std::string StepCut
BurnIn Cuts.
Definition: MCMCProcessor.h:373
MCMCProcessor::AutoCorrLag
int AutoCorrLag
LagL used in AutoCorrelation.
Definition: MCMCProcessor.h:515
MCMCProcessor::Posterior
std::unique_ptr< TCanvas > Posterior
Fancy canvas used for our beautiful plots.
Definition: MCMCProcessor.h:460
MCMCProcessor::OutputFile
TFile * OutputFile
The output file.
Definition: MCMCProcessor.h:457
MCMCProcessor::ApplySmoothing
bool ApplySmoothing
Apply smoothing for 2D histos using root algorithm.
Definition: MCMCProcessor.h:444
MCMCProcessor::nBins
int nBins
Number of bins.
Definition: MCMCProcessor.h:502
MCMCProcessor::Chain
TChain * Chain
Main chain storing all steps etc.
Definition: MCMCProcessor.h:371
MCMCProcessor::MCMCFile
std::string MCMCFile
Name of MCMC file.
Definition: MCMCProcessor.h:362
MCMCProcessor::CovPos
std::vector< std::vector< std::string > > CovPos
Covariance matrix name position.
Definition: MCMCProcessor.h:366
MCMCProcessor::ParamSums
double * ParamSums
Total parameter sum for each param.
Definition: MCMCProcessor.h:518
MCMCProcessor::UpperCut
int UpperCut
KS: Used only for SubOptimality.
Definition: MCMCProcessor.h:377
MCMCProcessor::Posterior1DCut
std::string Posterior1DCut
Cut used when making 1D Posterior distribution.
Definition: MCMCProcessor.h:375
MCMCProcessor::AccProbValues
double * AccProbValues
Holds all accProb.
Definition: MCMCProcessor.h:528
MCMCProcessor::hviolin
std::unique_ptr< TH2D > hviolin
Holds violin plot for all dials.
Definition: MCMCProcessor.h:492
MCMCProcessor::nDraw
int nDraw
Number of all parameters used in the analysis.
Definition: MCMCProcessor.h:391
MCMCProcessor::OutputSuffix
std::string OutputSuffix
Output file suffix useful when running over same file with different settings.
Definition: MCMCProcessor.h:364
MCMCProcessor::plotBinValue
bool plotBinValue
If true it will print value on each bin of covariance matrix.
Definition: MCMCProcessor.h:442
MCMCProcessor::nSteps
int nSteps
KS: For merged chains number of entries will be different from nSteps.
Definition: MCMCProcessor.h:385
MCMCProcessor::BatchedAverages
double ** BatchedAverages
Values of batched average for every param and batch.
Definition: MCMCProcessor.h:520
MCMCProcessor::PlotFlatPrior
bool PlotFlatPrior
Whether we plot flat prior or not, we usually provide error even for flat prior params.
Definition: MCMCProcessor.h:430
MCMCProcessor::FancyPlotNames
bool FancyPlotNames
Whether we want fancy plot names or not.
Definition: MCMCProcessor.h:440
MCMCProcessor::printToPDF
bool printToPDF
Will plot all plot to PDF not only to root file.
Definition: MCMCProcessor.h:438
MCMCProcessor::ParamNames
std::vector< std::vector< TString > > ParamNames
Name of parameters which we are going to analyse.
Definition: MCMCProcessor.h:401
MCMCProcessor::doDiagMCMC
bool doDiagMCMC
Doing MCMC Diagnostic.
Definition: MCMCProcessor.h:509
MCMCProcessor::nSysts
int nSysts
Number of covariance objects.
Definition: MCMCProcessor.h:389
MCMCProcessor::CacheMCMC
bool CacheMCMC
MCMC Chain has been cached.
Definition: MCMCProcessor.h:507
MCMCProcessor::ParamTypeStartPos
std::vector< int > ParamTypeStartPos
Definition: MCMCProcessor.h:415
MCMCProcessor::MadePostfit
bool MadePostfit
Sanity check if Postfit is already done to not make several times.
Definition: MCMCProcessor.h:436
MCMCProcessor::nEntries
int nEntries
KS: For merged chains number of entries will be different from nSteps.
Definition: MCMCProcessor.h:383
MCMCProcessor::plotRelativeToPrior
bool plotRelativeToPrior
Whether we plot relative to prior or nominal, in most cases is prior.
Definition: MCMCProcessor.h:434
MaCh3Utils::MaCh3Welcome
void MaCh3Welcome()
KS: Prints welcome message with MaCh3 logo.
Definition: Monitor.cpp:11

◆ ~MCMCProcessor()

MCMCProcessor::~MCMCProcessor ( )
virtual

Destroys the MCMCProcessor object.

Definition at line 106 of file MCMCProcessor.cpp.

106 {
107// ****************************
108 // Close the pdf file
109 MACH3LOG_INFO("Closing pdf in MCMCProcessor: {}", CanvasName.Data());
110 CanvasName += "]";
111 if(printToPDF) Posterior->Print(CanvasName);
112
113 delete Covariance;
114 delete Correlation;
115 delete Central_Value;
116 delete Means;
117 delete Errors;
118 delete Means_Gauss;
119 delete Errors_Gauss;
120 delete Means_HPD;
121 delete Errors_HPD;
122 delete Errors_HPD_Positive;
123 delete Errors_HPD_Negative;
124
125 for (int i = 0; i < nDraw; ++i)
126 {
127 if(hpost[i] != nullptr) delete hpost[i];
128 }
129 if(CacheMCMC)
130 {
131 for (int i = 0; i < nDraw; ++i)
132 {
133 for (int j = 0; j < nDraw; ++j)
134 {
135 delete hpost2D[i][j];
136 }
137 delete[] ParStep[i];
138 }
139 delete[] ParStep;
140 }
141 if(StepNumber != nullptr) delete[] StepNumber;
142
143 if(OutputFile != nullptr) OutputFile->Close();
144 if(OutputFile != nullptr) delete OutputFile;
145 delete Chain;
146}
MCMCProcessor::Correlation
TMatrixDSym * Correlation
Posterior Correlation Matrix.
Definition: MCMCProcessor.h:485
MCMCProcessor::Means_HPD
TVectorD * Means_HPD
Vector with mean values using Highest Posterior Density.
Definition: MCMCProcessor.h:474
MCMCProcessor::Errors
TVectorD * Errors
Vector with errors values using RMS.
Definition: MCMCProcessor.h:468
MCMCProcessor::Means_Gauss
TVectorD * Means_Gauss
Vector with mean values using Gaussian fit.
Definition: MCMCProcessor.h:470
MCMCProcessor::hpost
std::vector< TH1D * > hpost
Holds 1D Posterior Distributions.
Definition: MCMCProcessor.h:488
MCMCProcessor::Errors_HPD_Negative
TVectorD * Errors_HPD_Negative
Vector with negative error (left hand side) values using Highest Posterior Density.
Definition: MCMCProcessor.h:480
MCMCProcessor::Errors_Gauss
TVectorD * Errors_Gauss
Vector with error values using Gaussian fit.
Definition: MCMCProcessor.h:472
MCMCProcessor::Central_Value
TVectorD * Central_Value
Vector with central value for each parameter.
Definition: MCMCProcessor.h:464
MCMCProcessor::CanvasName
TString CanvasName
Name of canvas which help to save to the sample pdf.
Definition: MCMCProcessor.h:427
MCMCProcessor::Errors_HPD
TVectorD * Errors_HPD
Vector with error values using Highest Posterior Density.
Definition: MCMCProcessor.h:476
MCMCProcessor::hpost2D
std::vector< std::vector< TH2D * > > hpost2D
Holds 2D Posterior Distributions.
Definition: MCMCProcessor.h:490
MCMCProcessor::Means
TVectorD * Means
Vector with mean values using Arithmetic Mean.
Definition: MCMCProcessor.h:466
MCMCProcessor::Errors_HPD_Positive
TVectorD * Errors_HPD_Positive
Vector with positive error (right hand side) values using Highest Posterior Density.
Definition: MCMCProcessor.h:478
MCMCProcessor::Covariance
TMatrixDSym * Covariance
Posterior Covariance Matrix.
Definition: MCMCProcessor.h:483

Member Function Documentation

◆ AcceptanceProbabilities()

void MCMCProcessor::AcceptanceProbabilities ( )
inlineprotected

Acceptance Probability.

Definition at line 4030 of file MCMCProcessor.cpp.

4030 {
4031// **************************
4032 if (AccProbBatchedAverages == nullptr) PrepareDiagMCMC();
4033
4034 MACH3LOG_INFO("Making AccProb plots...");
4035
4036 // Set the titles and limits for TH1Ds
4037 auto AcceptanceProbPlot = std::make_unique<TH1D>("AcceptanceProbability", "Acceptance Probability", nEntries, 0, nEntries);
4038 AcceptanceProbPlot->GetXaxis()->SetTitle("Step");
4039 AcceptanceProbPlot->GetYaxis()->SetTitle("Acceptance Probability");
4040
4041 auto BatchedAcceptanceProblot = std::make_unique<TH1D>("AcceptanceProbability_Batch", "AcceptanceProbability_Batch", nBatches, 0, nBatches);
4042 BatchedAcceptanceProblot->GetYaxis()->SetTitle("Acceptance Probability");
4043
4044 for (int i = 0; i < nBatches; ++i) {
4045 BatchedAcceptanceProblot->SetBinContent(i+1, AccProbBatchedAverages[i]);
4046 const int BatchRangeLow = double(i)*double(nEntries)/double(nBatches);
4047 const int BatchRangeHigh = double(i+1)*double(nEntries)/double(nBatches);
4048 std::stringstream ss;
4049 ss << BatchRangeLow << " - " << BatchRangeHigh;
4050 BatchedAcceptanceProblot->GetXaxis()->SetBinLabel(i+1, ss.str().c_str());
4051 }
4052
4053 #ifdef MULTITHREAD
4054 #pragma omp parallel for
4055 #endif
4056 for (int i = 0; i < nEntries; ++i) {
4057 // Set bin content for the i-th bin to the parameter values
4058 AcceptanceProbPlot->SetBinContent(i, AccProbValues[i]);
4059 }
4060
4061 TDirectory *probDir = OutputFile->mkdir("AccProb");
4062 probDir->cd();
4063
4064 AcceptanceProbPlot->Write();
4065 BatchedAcceptanceProblot->Write();
4066 delete[] AccProbValues;
4067 delete[] AccProbBatchedAverages;
4068
4069 probDir->Close();
4070 delete probDir;
4071
4072 OutputFile->cd();
4073}
MCMCProcessor::PrepareDiagMCMC
void PrepareDiagMCMC()
CW: Prepare branches etc. for DiagMCMC.
Definition: MCMCProcessor.cpp:2936

◆ AutoCorrelation()

void MCMCProcessor::AutoCorrelation ( )
inlineprotected

KS: Calculate autocorrelations supports both OpenMP and CUDA :)

Definition at line 3287 of file MCMCProcessor.cpp.

3287 {
3288// *********************************
3289 if (ParStep == nullptr) PrepareDiagMCMC();
3290
3291 TStopwatch clock;
3292 clock.Start();
3293 const int nLags = AutoCorrLag;
3294 MACH3LOG_INFO("Making auto-correlations for nLags = {}", nLags);
3295
3296 // The sum of (Y-Ymean)^2 over all steps for each parameter
3297 std::vector<std::vector<double>> DenomSum(nDraw);
3298 std::vector<std::vector<double>> NumeratorSum(nDraw);
3299 std::vector<std::vector<double>> LagL(nDraw);
3300 for (int j = 0; j < nDraw; ++j) {
3301 DenomSum[j].resize(nLags);
3302 NumeratorSum[j].resize(nLags);
3303 LagL[j].resize(nLags);
3304 }
3305 std::vector<TH1D*> LagKPlots(nDraw);
3306 // Loop over the parameters of interest
3307 for (int j = 0; j < nDraw; ++j)
3308 {
3309 // Loop over each lag
3310 for (int k = 0; k < nLags; ++k) {
3311 NumeratorSum[j][k] = 0.0;
3312 DenomSum[j][k] = 0.0;
3313 LagL[j][k] = 0.0;
3314 }
3315
3316 // Make TH1Ds for each parameter which hold the lag
3317 TString Title = "";
3318 double Prior = 1.0, PriorError = 1.0;
3319
3320 GetNthParameter(j, Prior, PriorError, Title);
3321 std::string HistName = Form("%s_%s_Lag", Title.Data(), BranchNames[j].Data());
3322 LagKPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nLags, 0.0, nLags);
3323 LagKPlots[j]->GetXaxis()->SetTitle("Lag");
3324 LagKPlots[j]->GetYaxis()->SetTitle("Auto-correlation function");
3325 }
3326//KS: If CUDA is not enabled do calculations on CPU
3327#ifndef MaCh3_CUDA
3328 // Loop over the lags
3329 //CW: Each lag is independent so might as well multi-thread them!
3330 #ifdef MULTITHREAD
3331 MACH3LOG_INFO("Using multi-threading...");
3332 #pragma omp parallel for collapse(2)
3333 #endif // Loop over the number of parameters
3334 for (int j = 0; j < nDraw; ++j) {
3335 for (int k = 0; k < nLags; ++k) {
3336 // Loop over the number of entries
3337 for (int i = 0; i < nEntries; ++i) {
3338 const double Diff = ParStep[j][i]-ParamSums[j];
3339
3340 // Only sum the numerator up to i = N-k
3341 if (i < nEntries-k) {
3342 const double LagTerm = ParStep[j][i+k]-ParamSums[j];
3343 const double Product = Diff*LagTerm;
3344 NumeratorSum[j][k] += Product;
3345 }
3346 // Square the difference to form the denominator
3347 const double Denom = Diff*Diff;
3348 DenomSum[j][k] += Denom;
3349 }
3350 }
3351 }
3352#else //NOW GPU specific code
3353 MACH3LOG_INFO("Using GPU");
3354 //KS: This allocates memory and copy data from CPU to GPU
3355 PrepareGPU_AutoCorr(nLags);
3356
3357 //KS: This runs the main kernel and copy results back to CPU
3358 RunGPU_AutoCorr(
3359 ParStep_gpu,
3360 ParamSums_gpu,
3361 NumeratorSum_gpu,
3362 DenomSum_gpu,
3363 NumeratorSum_cpu,
3364 DenomSum_cpu);
3365
3366 #ifdef MULTITHREAD
3367 #pragma omp parallel for collapse(2)
3368 #endif
3369 //KS: Now that that we received data from GPU convert it to CPU-like format
3370 for (int j = 0; j < nDraw; ++j)
3371 {
3372 for (int k = 0; k < nLags; ++k)
3373 {
3374 const int temp_index = j*nLags+k;
3375 NumeratorSum[j][k] = NumeratorSum_cpu[temp_index];
3376 DenomSum[j][k] = DenomSum_cpu[temp_index];
3377 }
3378 }
3379 //delete auxiliary variables
3380 delete[] NumeratorSum_cpu;
3381 delete[] DenomSum_cpu;
3382 delete[] ParStep_cpu;
3383 delete[] ParamSums_cpu;
3384
3385 //KS: Delete stuff at GPU as well
3386 CleanupGPU_AutoCorr(
3387 ParStep_gpu,
3388 NumeratorSum_gpu,
3389 ParamSums_gpu,
3390 DenomSum_gpu);
3391
3392//KS: End of GPU specific code
3393#endif
3394
3395 OutputFile->cd();
3396 TDirectory *AutoCorrDir = OutputFile->mkdir("Auto_corr");
3397 // Now fill the LagK auto-correlation plots
3398 for (int j = 0; j < nDraw; ++j) {
3399 for (int k = 0; k < nLags; ++k) {
3400 LagL[j][k] = NumeratorSum[j][k]/DenomSum[j][k];
3401 LagKPlots[j]->SetBinContent(k, NumeratorSum[j][k]/DenomSum[j][k]);
3402 }
3403 AutoCorrDir->cd();
3404 LagKPlots[j]->Write();
3405 delete LagKPlots[j];
3406 }
3407
3408 //KS: This is different diagnostic however it relies on calculated Lag, thus we call it before we delete LagKPlots
3409 CalculateESS(nLags, LagL);
3410
3411 delete[] ParamSums;
3412
3413 AutoCorrDir->Close();
3414 delete AutoCorrDir;
3415
3416 OutputFile->cd();
3417
3418 clock.Stop();
3419 MACH3LOG_INFO("Making auto-correlations took {:.2f}s", clock.RealTime());
3420}
MCMCProcessor::GetNthParameter
void GetNthParameter(const int param, double &Prior, double &PriorError, TString &Title) const
Get properties of parameter by passing it number.
Definition: MCMCProcessor.cpp:2425
MCMCProcessor::CalculateESS
void CalculateESS(const int nLags, const std::vector< std::vector< double > > &LagL)
KS: calc Effective Sample Size.
Definition: MCMCProcessor.cpp:3506
MCMCProcessor::BranchNames
std::vector< TString > BranchNames
Definition: MCMCProcessor.h:394
RunGPU_AutoCorr
__host__ void RunGPU_AutoCorr(float *ParStep_gpu, float *ParamSums_gpu, float *NumeratorSum_gpu, float *DenomSum_gpu, float *NumeratorSum_cpu, float *DenomSum_cpu)
KS: This call the main kernel responsible for calculating LagL and later copy results back to CPU.
Definition: gpuMCMCProcessorUtils.cu:160
CleanupGPU_AutoCorr
__host__ void CleanupGPU_AutoCorr(float *ParStep_gpu, float *NumeratorSum_gpu, float *ParamSums_gpu, float *DenomSum_gpu)
KS: free memory on gpu.
Definition: gpuMCMCProcessorUtils.cu:197

◆ AutoCorrelation_FFT()

void MCMCProcessor::AutoCorrelation_FFT ( )
inlineprotected

MJR: Autocorrelation function using FFT algorithm for extra speed.

Author
Michael Reh

Definition at line 3192 of file MCMCProcessor.cpp.

3192 {
3193// *********************************
3194 if (ParStep == nullptr) PrepareDiagMCMC();
3195
3196 TStopwatch clock;
3197 clock.Start();
3198 const int nLags = AutoCorrLag;
3199 MACH3LOG_INFO("Making auto-correlations for nLags = {}", nLags);
3200
3201 // Prep outputs
3202 OutputFile->cd();
3203 TDirectory* AutoCorrDir = OutputFile->mkdir("Auto_corr");
3204 std::vector<TH1D*> LagKPlots(nDraw);
3205 std::vector<std::vector<double>> LagL(nDraw);
3206
3207 // Arrays needed to perform FFT using ROOT
3208 std::vector<double> ACFFT(nEntries, 0.0); // Main autocorrelation array
3209 std::vector<double> ParVals(nEntries, 0.0); // Param values for full chain
3210 std::vector<double> ParValsFFTR(nEntries, 0.0); // FFT Real part
3211 std::vector<double> ParValsFFTI(nEntries, 0.0); // FFT Imaginary part
3212 std::vector<double> ParValsFFTSquare(nEntries, 0.0); // FFT Absolute square
3213 std::vector<double> ParValsComplex(nEntries, 0.0); // Input Imaginary values (0)
3214
3215 // Create forward and reverse FFT objects. I don't love using ROOT here,
3216 // but it works so I can't complain
3217 TVirtualFFT* fftf = TVirtualFFT::FFT(1, &nEntries, "C2CFORWARD");
3218 TVirtualFFT* fftb = TVirtualFFT::FFT(1, &nEntries, "C2CBACKWARD");
3219
3220 // Loop over all pars and calculate the full autocorrelation function using FFT
3221 for (int j = 0; j < nDraw; ++j) {
3222 // Initialize
3223 LagL[j].resize(nLags);
3224 for (int i = 0; i < nEntries; ++i) {
3225 ParVals[i] = ParStep[j][i]-ParamSums[j]; // Subtract the mean to make it numerically tractable
3226 ParValsComplex[i] = 0.; // Reset dummy array
3227 }
3228
3229 // Transform
3230 fftf->SetPointsComplex(ParVals.data(), ParValsComplex.data());
3231 fftf->Transform();
3232 fftf->GetPointsComplex(ParValsFFTR.data(), ParValsFFTI.data());
3233
3234 // Square the results to get the power spectrum
3235 for (int i = 0; i < nEntries; ++i) {
3236 ParValsFFTSquare[i] = ParValsFFTR[i]*ParValsFFTR[i] + ParValsFFTI[i]*ParValsFFTI[i];
3237 }
3238
3239 // Transforming back gives the autocovariance
3240 fftb->SetPointsComplex(ParValsFFTSquare.data(), ParValsComplex.data());
3241 fftb->Transform();
3242 fftb->GetPointsComplex(ACFFT.data(), ParValsComplex.data());
3243
3244 // Divide by norm to get autocorrelation
3245 double normAC = ACFFT[0];
3246 for (int i = 0; i < nEntries; ++i) {
3247 ACFFT[i] /= normAC;
3248 }
3249
3250 // Get plotting info
3251 TString Title = "";
3252 double Prior = 1.0, PriorError = 1.0;
3253 GetNthParameter(j, Prior, PriorError, Title);
3254 std::string HistName = Form("%s_%s_Lag", Title.Data(), BranchNames[j].Data());
3255
3256 // Initialize Lag plot
3257 LagKPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nLags, 0.0, nLags);
3258 LagKPlots[j]->GetXaxis()->SetTitle("Lag");
3259 LagKPlots[j]->GetYaxis()->SetTitle("Auto-correlation function");
3260
3261 // Fill plot
3262 for (int k = 0; k < nLags; ++k) {
3263 LagL[j][k] = ACFFT[k];
3264 LagKPlots[j]->SetBinContent(k, ACFFT[k]);
3265 }
3266
3267 // Write and clean up
3268 AutoCorrDir->cd();
3269 LagKPlots[j]->Write();
3270 delete LagKPlots[j];
3271 }
3272
3273 //KS: This is different diagnostic however it relies on calculated Lag, thus we call it before we delete LagKPlots
3274 CalculateESS(nLags, LagL);
3275
3276 AutoCorrDir->Close();
3277 delete AutoCorrDir;
3278
3279 OutputFile->cd();
3280
3281 clock.Stop();
3282 MACH3LOG_INFO("Making auto-correlations took {:.2f}s", clock.RealTime());
3283}

◆ BatchedAnalysis()

void MCMCProcessor::BatchedAnalysis ( )
inlineprotected

Get the batched means variance estimation and variable indicating if number of batches is sensible [4] [25].

Definition at line 3661 of file MCMCProcessor.cpp.

3661 {
3662// **************************
3663 if(BatchedAverages == nullptr)
3664 {
3665 MACH3LOG_ERROR("BatchedAverages haven't been initialises or have been deleted something is wrong");
3666 MACH3LOG_ERROR("I need it and refuse to go further");
3667 throw MaCh3Exception(__FILE__ , __LINE__ );
3668 }
3669
3670 // Calculate variance estimator using batched means following @cite chakraborty2019estimating see Eq. 1.2
3671 TVectorD* BatchedVariance = new TVectorD(nDraw);
3672 //KS: The hypothesis is rejected if C > z α for a given confidence level α. If the batch means do not pass the test, Correlated is reported for the half-width on the statistical reports following @cite rossetti2024batch alternatively for more old-school see Alexopoulos and Seila 1998 section 3.4.3
3673 TVectorD* C_Test_Statistics = new TVectorD(nDraw);
3674
3675 std::vector<double> OverallBatchMean(nDraw);
3676 std::vector<double> C_Rho_Nominator(nDraw);
3677 std::vector<double> C_Rho_Denominator(nDraw);
3678 std::vector<double> C_Nominator(nDraw);
3679 std::vector<double> C_Denominator(nDraw);
3680 const int BatchLength = nEntries/nBatches+1;
3681//KS: Start parallel region
3682#ifdef MULTITHREAD
3683#pragma omp parallel
3684{
3685#endif
3686 #ifdef MULTITHREAD
3687 #pragma omp for
3688 #endif
3689 //KS: First calculate mean of batched means for each param and Initialise everything to 0
3690 for (int j = 0; j < nDraw; ++j)
3691 {
3692 OverallBatchMean[j] = 0.0;
3693 C_Rho_Nominator[j] = 0.0;
3694 C_Rho_Denominator[j] = 0.0;
3695 C_Nominator[j] = 0.0;
3696 C_Denominator[j] = 0.0;
3697
3698 (*BatchedVariance)(j) = 0.0;
3699 (*C_Test_Statistics)(j) = 0.0;
3700 for (int i = 0; i < nBatches; ++i)
3701 {
3702 OverallBatchMean[j] += BatchedAverages[i][j];
3703 }
3704 OverallBatchMean[j] /= nBatches;
3705 }
3706
3707 #ifdef MULTITHREAD
3708 #pragma omp for nowait
3709 #endif
3710 //KS: next loop is completely independent thus nowait clause
3711 for (int j = 0; j < nDraw; ++j)
3712 {
3713 for (int i = 0; i < nBatches; ++i)
3714 {
3715 (*BatchedVariance)(j) += (OverallBatchMean[j] - BatchedAverages[i][j])*(OverallBatchMean[j] - BatchedAverages[i][j]);
3716 }
3717 (*BatchedVariance)(j) = (BatchLength/(nBatches-1))* (*BatchedVariance)(j);
3718 }
3719
3720 //KS: Now we focus on C test statistic, again use nowait as next is calculation is independent
3721 #ifdef MULTITHREAD
3722 #pragma omp for nowait
3723 #endif
3724 for (int j = 0; j < nDraw; ++j)
3725 {
3726 C_Nominator[j] = (OverallBatchMean[j] - BatchedAverages[0][j])*(OverallBatchMean[j] - BatchedAverages[0][j]) +
3727 (OverallBatchMean[j] - BatchedAverages[nBatches-1][j])*(OverallBatchMean[j] - BatchedAverages[nBatches-1][j]);
3728 for (int i = 0; i < nBatches; ++i)
3729 {
3730 C_Denominator[j] += (OverallBatchMean[j] - BatchedAverages[i][j])*(OverallBatchMean[j] - BatchedAverages[i][j]);
3731 }
3732 C_Denominator[j] = 2*C_Denominator[j];
3733 }
3734
3735 //KS: We still calculate C and for this we need rho wee need autocorrelations between batches
3736 #ifdef MULTITHREAD
3737 #pragma omp for
3738 #endif
3739 for (int j = 0; j < nDraw; ++j)
3740 {
3741 for (int i = 0; i < nBatches-1; ++i)
3742 {
3743 C_Rho_Nominator[j] += (OverallBatchMean[j] - BatchedAverages[i][j])*(OverallBatchMean[j] - BatchedAverages[i+1][j]);
3744 }
3745
3746 for (int i = 0; i < nBatches; ++i)
3747 {
3748 C_Rho_Denominator[j] += (OverallBatchMean[j] - BatchedAverages[i][j])*(OverallBatchMean[j] - BatchedAverages[i][j]);
3749 }
3750 }
3751
3752 //KS: Final calculations of C
3753 #ifdef MULTITHREAD
3754 #pragma omp for
3755 #endif
3756 for (int j = 0; j < nDraw; ++j)
3757 {
3758 (*C_Test_Statistics)(j) = std::sqrt((nBatches*nBatches - 1)/(nBatches-2)) * ( C_Rho_Nominator[j]/C_Rho_Denominator[j] + C_Nominator[j]/ C_Denominator[j]);
3759 }
3760#ifdef MULTITHREAD
3761} //End parallel region
3762#endif
3763
3764 //Save to file
3765 OutputFile->cd();
3766 BatchedVariance->Write("BatchedMeansVariance");
3767 C_Test_Statistics->Write("C_Test_Statistics");
3768
3769 //Delete all variables
3770 delete BatchedVariance;
3771 delete C_Test_Statistics;
3772}
MACH3LOG_ERROR
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
MaCh3Exception
Custom exception class for MaCh3 errors.
Definition: MaCh3Exception.h:13

◆ BatchedMeans()

void MCMCProcessor::BatchedMeans ( )
inlineprotected

CW: Batched means, literally read from an array and chuck into TH1D.

Definition at line 3603 of file MCMCProcessor.cpp.

3603 {
3604// **************************
3605 if (BatchedAverages == nullptr) PrepareDiagMCMC();
3606 MACH3LOG_INFO("Making BatchedMeans plots...");
3607
3608 std::vector<TH1D*> BatchedParamPlots(nDraw);
3609 for (int j = 0; j < nDraw; ++j) {
3610 TString Title = "";
3611 double Prior = 1.0, PriorError = 1.0;
3612
3613 GetNthParameter(j, Prior, PriorError, Title);
3614
3615 std::string HistName = Form("%s_%s_batch", Title.Data(), BranchNames[j].Data());
3616 BatchedParamPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nBatches, 0, nBatches);
3617 }
3618
3619 #ifdef MULTITHREAD
3620 #pragma omp parallel for
3621 #endif
3622 for (int j = 0; j < nDraw; ++j) {
3623 for (int i = 0; i < nBatches; ++i) {
3624 BatchedParamPlots[j]->SetBinContent(i+1, BatchedAverages[i][j]);
3625 const int BatchRangeLow = double(i)*double(nEntries)/double(nBatches);
3626 const int BatchRangeHigh = double(i+1)*double(nEntries)/double(nBatches);
3627 std::stringstream ss;
3628 ss << BatchRangeLow << " - " << BatchRangeHigh;
3629 BatchedParamPlots[j]->GetXaxis()->SetBinLabel(i+1, ss.str().c_str());
3630 }
3631 }
3632
3633 TDirectory *BatchDir = OutputFile->mkdir("Batched_means");
3634 BatchDir->cd();
3635 for (int j = 0; j < nDraw; ++j) {
3636 TF1 *Fitter = new TF1("Fitter","[0]", 0, nBatches);
3637 Fitter->SetLineColor(kRed);
3638 BatchedParamPlots[j]->Fit("Fitter","Rq");
3639 BatchedParamPlots[j]->Write();
3640 delete Fitter;
3641 delete BatchedParamPlots[j];
3642 }
3643
3644 //KS: Get the batched means variance estimation and variable indicating if number of batches is sensible
3645 // We do this before deleting BatchedAverages
3646 BatchedAnalysis();
3647
3648 for (int i = 0; i < nBatches; ++i) {
3649 delete BatchedAverages[i];
3650 }
3651 delete[] BatchedAverages;
3652
3653 BatchDir->Close();
3654 delete BatchDir;
3655
3656 OutputFile->cd();
3657}
MCMCProcessor::BatchedAnalysis
void BatchedAnalysis()
Get the batched means variance estimation and variable indicating if number of batches is sensible .
Definition: MCMCProcessor.cpp:3661

◆ CacheSteps()

void MCMCProcessor::CacheSteps ( )

KS:By caching each step we use multithreading.

Definition at line 1032 of file MCMCProcessor.cpp.

1032 {
1033// ***************
1034 if(CacheMCMC == true) return;
1035
1036 CacheMCMC = true;
1037
1038 if(ParStep != nullptr)
1039 {
1040 MACH3LOG_ERROR("It look like ParStep was already filled ");
1041 MACH3LOG_ERROR("Even though it is used for MakeCovariance_MP and for DiagMCMC ");
1042 MACH3LOG_ERROR("it has different structure in both for cache hits, sorry ");
1043 throw MaCh3Exception(__FILE__ , __LINE__ );
1044 }
1045
1046 MACH3LOG_INFO("Caching input tree...");
1047 MACH3LOG_INFO("Allocating {:.2f} MB", double(sizeof(double)*nDraw*nEntries)/1.E6);
1048 TStopwatch clock;
1049 clock.Start();
1050
1051 ParStep = new double*[nDraw];
1052 StepNumber = new int[nEntries];
1053
1054 hpost2D.resize(nDraw);
1055 for (int i = 0; i < nDraw; ++i)
1056 {
1057 ParStep[i] = new double[nEntries];
1058 hpost2D[i].resize(nDraw);
1059 for (int j = 0; j < nEntries; ++j)
1060 {
1061 ParStep[i][j] = -999.99;
1062 //KS: Set this only once
1063 if(i == 0) StepNumber[j] = -999.99;
1064 }
1065 }
1066
1067 // Set all the branches to off
1068 Chain->SetBranchStatus("*", false);
1069 int stepBranch = 0;
1070 double* ParValBranch = new double[nEntries]();
1071 // Turn on the branches which we want for parameters
1072 for (int i = 0; i < nDraw; ++i)
1073 {
1074 Chain->SetBranchStatus(BranchNames[i].Data(), true);
1075 Chain->SetBranchAddress(BranchNames[i].Data(), &ParValBranch[i]);
1076 }
1077 Chain->SetBranchStatus("step", true);
1078 Chain->SetBranchAddress("step", &stepBranch);
1079 const Long64_t countwidth = nEntries/10;
1080
1081 // Loop over the entries
1082 //KS: This is really a bottleneck right now, thus revisit with ROOT6 https://pep-root6.github.io/docs/analysis/parallell/root.html
1083 for (Long64_t j = 0; j < nEntries; ++j)
1084 {
1085 if (j % countwidth == 0) {
1086 MaCh3Utils::PrintProgressBar(j, nEntries);
1087 MaCh3Utils::EstimateDataTransferRate(Chain, j);
1088 } else {
1089 Chain->GetEntry(j);
1090 }
1091 StepNumber[j] = stepBranch;
1092 // Set the branch addresses for params
1093 for (int i = 0; i < nDraw; ++i)
1094 {
1095 ParStep[i][j] = ParValBranch[i];
1096 }
1097 }
1098 delete[] ParValBranch;
1099
1100 // Set all the branches to on
1101 Chain->SetBranchStatus("*", true);
1102
1103 // KS: Set temporary branch address to allow min/max, otherwise ROOT can segfaults
1104 double tempVal = 0.0;
1105 std::vector<double> Min_Chain(nDraw);
1106 std::vector<double> Max_Chain(nDraw);
1107 for (int i = 0; i < nDraw; ++i)
1108 {
1109 Chain->SetBranchAddress(BranchNames[i].Data(), &tempVal);
1110 Min_Chain[i] = Chain->GetMinimum(BranchNames[i]);
1111 Max_Chain[i] = Chain->GetMaximum(BranchNames[i]);
1112 }
1113
1114 // Cache max and min in chain for covariance matrix
1115 for (int i = 0; i < nDraw; ++i)
1116 {
1117 TString Title_i = "";
1118 double Prior_i, PriorError_i;
1119 GetNthParameter(i, Prior_i, PriorError_i, Title_i);
1120
1121 for (int j = 0; j <= i; ++j)
1122 {
1123 // TH2D to hold the Correlation
1124 hpost2D[i][j] = new TH2D(Form("hpost2D_%i_%i",i,j), Form("hpost2D_%i_%i",i,j), nBins, Min_Chain[i], Max_Chain[i], nBins, Min_Chain[j], Max_Chain[j]);
1125 TString Title_j = "";
1126 double Prior_j, PriorError_j;
1127 GetNthParameter(j, Prior_j, PriorError_j, Title_j);
1128
1129 hpost2D[i][j]->SetMinimum(0);
1130 hpost2D[i][j]->GetXaxis()->SetTitle(Title_i);
1131 hpost2D[i][j]->GetYaxis()->SetTitle(Title_j);
1132 hpost2D[i][j]->GetZaxis()->SetTitle("Steps");
1133 }
1134 }
1135 clock.Stop();
1136 MACH3LOG_INFO("Caching steps took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries );
1137}
MaCh3Utils::PrintProgressBar
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:212
MaCh3Utils::EstimateDataTransferRate
void EstimateDataTransferRate(TChain *chain, const Long64_t entry)
KS: Check what CPU you are using.
Definition: Monitor.cpp:195

◆ CalculateESS()

void MCMCProcessor::CalculateESS ( const int  nLags,
const std::vector< std::vector< double > > &  LagL 
)
inlineprotected

KS: calc Effective Sample Size.

Parameters
nLagsShould be the same nLags as used in AutoCorrelation()
LagLValue of LagL for each dial and each Lag

This function computes the Effective Sample Size (ESS) using the autocorrelations calculated by AutoCorrelation(). Ensure that the parameter nLags here matches the number of lags used in AutoCorrelation() to obtain accurate results. [27] [13] [9]

Definition at line 3506 of file MCMCProcessor.cpp.

3506 {
3507// **************************
3508 if(LagL.size() == 0)
3509 {
3510 MACH3LOG_ERROR("Size of LagL is 0");
3511 throw MaCh3Exception(__FILE__ , __LINE__ );
3512 }
3513 MACH3LOG_INFO("Making ESS plots...");
3514 TVectorD* EffectiveSampleSize = new TVectorD(nDraw);
3515 TVectorD* SamplingEfficiency = new TVectorD(nDraw);
3516 std::vector<double> TempDenominator(nDraw);
3517
3518 constexpr int Nhists = 5;
3519 constexpr double Thresholds[Nhists + 1] = {1, 0.02, 0.005, 0.001, 0.0001, 0.0};
3520 constexpr Color_t ESSColours[Nhists] = {kGreen, kGreen + 2, kYellow, kOrange, kRed};
3521
3522 //KS: This histogram is inspired by the following: @cite gabry2024visual
3523 std::vector<std::unique_ptr<TH1D>> EffectiveSampleSizeHist(Nhists);
3524 for(int i = 0; i < Nhists; ++i)
3525 {
3526 EffectiveSampleSizeHist[i] =
3527 std::make_unique<TH1D>(Form("EffectiveSampleSizeHist_%i", i), Form("EffectiveSampleSizeHist_%i", i), nDraw, 0, nDraw);
3528 EffectiveSampleSizeHist[i]->GetYaxis()->SetTitle("N_{eff}/N");
3529 EffectiveSampleSizeHist[i]->SetFillColor(ESSColours[i]);
3530 EffectiveSampleSizeHist[i]->SetLineColor(ESSColours[i]);
3531 EffectiveSampleSizeHist[i]->Sumw2();
3532 for (int j = 0; j < nDraw; ++j)
3533 {
3534 TString Title = "";
3535 double Prior = 1.0, PriorError = 1.0;
3536 GetNthParameter(j, Prior, PriorError, Title);
3537 EffectiveSampleSizeHist[i]->GetXaxis()->SetBinLabel(j+1, Title.Data());
3538 }
3539 }
3540
3541 #ifdef MULTITHREAD
3542 #pragma omp parallel for
3543 #endif
3544 //KS: Calculate ESS and MCMC efficiency for each parameter
3545 for (int j = 0; j < nDraw; ++j)
3546 {
3547 (*EffectiveSampleSize)(j) = _UNDEF_;
3548 (*SamplingEfficiency)(j) = _UNDEF_;
3549 TempDenominator[j] = 0.;
3550 //KS: Firs sum over all Calculated autocorrelations
3551 for (int k = 0; k < nLags; ++k)
3552 {
3553 TempDenominator[j] += LagL[j][k];
3554 }
3555 TempDenominator[j] = 1+2*TempDenominator[j];
3556 (*EffectiveSampleSize)(j) = double(nEntries)/TempDenominator[j];
3557 // 100 because we convert to percentage
3558 (*SamplingEfficiency)(j) = 100 * 1/TempDenominator[j];
3559
3560 for(int i = 0; i < Nhists; ++i)
3561 {
3562 EffectiveSampleSizeHist[i]->SetBinContent(j+1, 0);
3563 EffectiveSampleSizeHist[i]->SetBinError(j+1, 0);
3564
3565 const double TempEntry = std::fabs((*EffectiveSampleSize)(j)) / double(nEntries);
3566 if(Thresholds[i] >= TempEntry && TempEntry > Thresholds[i+1])
3567 {
3568 if( std::isnan((*EffectiveSampleSize)(j)) ) continue;
3569 EffectiveSampleSizeHist[i]->SetBinContent(j+1, TempEntry);
3570 }
3571 }
3572 }
3573
3574 //KS Write to the output tree
3575 //Save to file
3576 OutputFile->cd();
3577 EffectiveSampleSize->Write("EffectiveSampleSize");
3578 SamplingEfficiency->Write("SamplingEfficiency");
3579
3580 EffectiveSampleSizeHist[0]->SetTitle("Effective Sample Size");
3581 EffectiveSampleSizeHist[0]->Draw();
3582 for(int i = 1; i < Nhists; ++i)
3583 {
3584 EffectiveSampleSizeHist[i]->Draw("SAME");
3585 }
3586
3587 auto leg = std::make_unique<TLegend>(0.2, 0.7, 0.6, 0.95);
3588 SetLegendStyle(leg.get(), 0.03);
3589 for(int i = 0; i < Nhists; ++i)
3590 {
3591 leg->AddEntry(EffectiveSampleSizeHist[i].get(), Form("%.4f >= N_{eff}/N > %.4f", Thresholds[i], Thresholds[i+1]), "f");
3592 } leg->Draw("SAME");
3593
3594 Posterior->Write("EffectiveSampleSizeCanvas");
3595
3596 //Delete all variables
3597 delete EffectiveSampleSize;
3598 delete SamplingEfficiency;
3599}
EffectiveSampleSize
double * EffectiveSampleSize
Definition: RHat.cpp:67
MCMCProcessor::SetLegendStyle
void SetLegendStyle(TLegend *Legend, const double size) const
Configures the style of a TLegend object.
Definition: MCMCProcessor.cpp:4190

◆ CheckCredibleIntervalsOrder()

void MCMCProcessor::CheckCredibleIntervalsOrder ( const std::vector< double > &  CredibleIntervals,
const std::vector< Color_t > &  CredibleIntervalsColours 
)

Checks the order and size consistency of the CredibleIntervals and CredibleIntervalsColours vectors.

Parameters
CredibleIntervalsA vector of credible interval values.
CredibleIntervalsColoursA vector of colors associated with each credible interval.
Exceptions
MaCh3ExceptionIf the sizes are not equal or the intervals are not in decreasing order.

Definition at line 4076 of file MCMCProcessor.cpp.

4076 {
4077// **************************
4078 if (CredibleIntervals.size() != CredibleIntervalsColours.size()) {
4079 MACH3LOG_ERROR("size of CredibleIntervals is not equal to size of CredibleIntervalsColours");
4080 throw MaCh3Exception(__FILE__, __LINE__);
4081 }
4082 if (CredibleIntervals.size() > 1) {
4083 for (unsigned int i = 1; i < CredibleIntervals.size(); i++) {
4084 if (CredibleIntervals[i] > CredibleIntervals[i - 1]) {
4085 MACH3LOG_ERROR("Interval {} is smaller than {}", i, i - 1);
4086 MACH3LOG_ERROR("{:.2f} {:.2f}", CredibleIntervals[i], CredibleIntervals[i - 1]);
4087 MACH3LOG_ERROR("They should be grouped in decreasing order");
4088 throw MaCh3Exception(__FILE__, __LINE__);
4089 }
4090 }
4091 }
4092}

◆ CheckCredibleRegionsOrder()

void MCMCProcessor::CheckCredibleRegionsOrder ( const std::vector< double > &  CredibleRegions,
const std::vector< Style_t > &  CredibleRegionStyle,
const std::vector< Color_t > &  CredibleRegionColor 
)

Checks the order and size consistency of the CredibleRegions, CredibleRegionStyle, and CredibleRegionColor vectors.

Parameters
CredibleRegionsA vector of credible region values.
CredibleRegionStyleA vector of styles associated with each credible region.
CredibleRegionColorA vector of colors associated with each credible region.
Exceptions
MaCh3ExceptionIf the sizes are not equal or the regions are not in decreasing order.

Definition at line 4095 of file MCMCProcessor.cpp.

4097 {
4098// **************************
4099 if ((CredibleRegions.size() != CredibleRegionStyle.size()) || (CredibleRegionStyle.size() != CredibleRegionColor.size())) {
4100 MACH3LOG_ERROR("size of CredibleRegions is not equal to size of CredibleRegionStyle or CredibleRegionColor");
4101 throw MaCh3Exception(__FILE__, __LINE__);
4102 }
4103 for (unsigned int i = 1; i < CredibleRegions.size(); i++) {
4104 if (CredibleRegions[i] > CredibleRegions[i - 1]) {
4105 MACH3LOG_ERROR("Interval {} is smaller than {}", i, i - 1);
4106 MACH3LOG_ERROR("{:.2f} {:.2f}", CredibleRegions[i], CredibleRegions[i - 1]);
4107 MACH3LOG_ERROR("They should be grouped in decreasing order");
4108 throw MaCh3Exception(__FILE__, __LINE__);
4109 }
4110 }
4111}

◆ DiagMCMC()

void MCMCProcessor::DiagMCMC ( )

KS: Perform MCMC diagnostic including Autocorrelation, Trace etc.

Definition at line 2906 of file MCMCProcessor.cpp.

2906 {
2907// **************************
2908 // Prepare branches etc for DiagMCMC
2909 PrepareDiagMCMC();
2910
2911 // Draw the simple trace matrices
2912 ParamTraces();
2913
2914 // Get the batched means
2915 BatchedMeans();
2916
2917 // Draw the auto-correlations
2918 if (useFFTAutoCorrelation) {
2919 AutoCorrelation_FFT();
2920 } else {
2921 AutoCorrelation();
2922 }
2923
2924 // Calculate Power Spectrum for each param
2925 PowerSpectrumAnalysis();
2926
2927 // Get Geweke Z score helping select burn-in
2928 GewekeDiagnostic();
2929
2930 // Draw acceptance Probability
2931 AcceptanceProbabilities();
2932}
MCMCProcessor::GewekeDiagnostic
void GewekeDiagnostic()
Geweke Diagnostic based on the methods described by Fang (2014) and Karlsbakk (2011)....
Definition: MCMCProcessor.cpp:3898
MCMCProcessor::AcceptanceProbabilities
void AcceptanceProbabilities()
Acceptance Probability.
Definition: MCMCProcessor.cpp:4030
MCMCProcessor::AutoCorrelation
void AutoCorrelation()
KS: Calculate autocorrelations supports both OpenMP and CUDA :)
Definition: MCMCProcessor.cpp:3287
MCMCProcessor::ParamTraces
void ParamTraces()
CW: Draw trace plots of the parameters i.e. parameter vs step.
Definition: MCMCProcessor.cpp:3099
MCMCProcessor::AutoCorrelation_FFT
void AutoCorrelation_FFT()
MJR: Autocorrelation function using FFT algorithm for extra speed.
Definition: MCMCProcessor.cpp:3192
MCMCProcessor::PowerSpectrumAnalysis
void PowerSpectrumAnalysis()
RC: Perform spectral analysis of MCMC .
Definition: MCMCProcessor.cpp:3776
MCMCProcessor::BatchedMeans
void BatchedMeans()
CW: Batched means, literally read from an array and chuck into TH1D.
Definition: MCMCProcessor.cpp:3603

◆ DrawCorrelations1D()

void MCMCProcessor::DrawCorrelations1D ( )
inlineprotected

Draw 1D correlations which might be more helpful than looking at huge 2D Corr matrix.

Definition at line 1410 of file MCMCProcessor.cpp.

1410 {
1411// *********************
1412 //KS: Store it as we go back to them at the end
1413 const std::vector<double> Margins = GetMargins(Posterior);
1414 const int OptTitle = gStyle->GetOptTitle();
1415
1416 Posterior->SetTopMargin(0.1);
1417 Posterior->SetBottomMargin(0.2);
1418 gStyle->SetOptTitle(1);
1419
1420 constexpr int Nhists = 3;
1421 //KS: Highest value is just meant bo be sliglhy higher than 1 to catch >,
1422 constexpr double Thresholds[Nhists+1] = {0, 0.25, 0.5, 1.0001};
1423 constexpr Color_t CorrColours[Nhists] = {kRed-10, kRed-6, kRed};
1424
1425 //KS: This store necessary entries for stripped covariance which store only "meaningful correlations
1426 std::vector<std::vector<double>> CorrOfInterest;
1427 CorrOfInterest.resize(nDraw);
1428 std::vector<std::vector<std::string>> NameCorrOfInterest;
1429 NameCorrOfInterest.resize(nDraw);
1430
1431 std::vector<std::vector<std::unique_ptr<TH1D>>> Corr1DHist(nDraw);
1432 //KS: Initialising ROOT objects is never safe in MP loop
1433 for(int i = 0; i < nDraw; ++i)
1434 {
1435 TString Title = "";
1436 double Prior = 1.0, PriorError = 1.0;
1437 GetNthParameter(i, Prior, PriorError, Title);
1438
1439 Corr1DHist[i].resize(Nhists);
1440 for(int j = 0; j < Nhists; ++j)
1441 {
1442 Corr1DHist[i][j] = std::make_unique<TH1D>(Form("Corr1DHist_%i_%i", i, j), Form("Corr1DHist_%i_%i", i, j), nDraw, 0, nDraw);
1443 Corr1DHist[i][j]->SetTitle(Form("%s",Title.Data()));
1444 Corr1DHist[i][j]->GetYaxis()->SetTitle("Correlation");
1445 Corr1DHist[i][j]->SetFillColor(CorrColours[j]);
1446 Corr1DHist[i][j]->SetLineColor(kBlack);
1447
1448 for (int k = 0; k < nDraw; ++k)
1449 {
1450 TString Title_y = "";
1451 double Prior_y = 1.0;
1452 double PriorError_y = 1.0;
1453 GetNthParameter(k, Prior_y, PriorError_y, Title_y);
1454 Corr1DHist[i][j]->GetXaxis()->SetBinLabel(k+1, Title_y.Data());
1455 }
1456 }
1457 }
1458
1459 #ifdef MULTITHREAD
1460 #pragma omp parallel for
1461 #endif
1462 for(int i = 0; i < nDraw; ++i)
1463 {
1464 for(int j = 0; j < nDraw; ++j)
1465 {
1466 for(int k = 0; k < Nhists; ++k)
1467 {
1468 const double TempEntry = std::fabs((*Correlation)(i,j));
1469 if(Thresholds[k+1] > TempEntry && TempEntry >= Thresholds[k])
1470 {
1471 Corr1DHist[i][k]->SetBinContent(j+1, (*Correlation)(i,j));
1472 }
1473 }
1474 if(std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold && i != j)
1475 {
1476 CorrOfInterest[i].push_back((*Correlation)(i,j));
1477 NameCorrOfInterest[i].push_back(Corr1DHist[i][0]->GetXaxis()->GetBinLabel(j+1));
1478 }
1479 }
1480 }
1481
1482 TDirectory *CorrDir = OutputFile->mkdir("Corr1D");
1483 CorrDir->cd();
1484
1485 for(int i = 0; i < nDraw; i++)
1486 {
1487 if (IamVaried[i] == false) continue;
1488
1489 Corr1DHist[i][0]->SetMaximum(+1.);
1490 Corr1DHist[i][0]->SetMinimum(-1.);
1491 Corr1DHist[i][0]->Draw();
1492 for(int k = 1; k < Nhists; k++) {
1493 Corr1DHist[i][k]->Draw("SAME");
1494 }
1495
1496 auto leg = std::make_unique<TLegend>(0.3, 0.75, 0.6, 0.90);
1497 SetLegendStyle(leg.get(), 0.02);
1498 for(int k = 0; k < Nhists; k++) {
1499 leg->AddEntry(Corr1DHist[i][k].get(), Form("%.2f > |Corr| >= %.2f", Thresholds[k+1], Thresholds[k]), "f");
1500 }
1501 leg->Draw("SAME");
1502
1503 Posterior->Write(Corr1DHist[i][0]->GetTitle());
1504 if(printToPDF) Posterior->Print(CanvasName);
1505 }
1506
1507 //KS: Plot only meaningful correlations
1508 for(int i = 0; i < nDraw; i++)
1509 {
1510 const int size = int(CorrOfInterest[i].size());
1511
1512 if(size == 0) continue;
1513 auto Corr1DHist_Reduced = std::make_unique<TH1D>("Corr1DHist_Reduced", "Corr1DHist_Reduced", size, 0, size);
1514 Corr1DHist_Reduced->SetTitle(Corr1DHist[i][0]->GetTitle());
1515 Corr1DHist_Reduced->GetYaxis()->SetTitle("Correlation");
1516 Corr1DHist_Reduced->SetFillColor(kBlue);
1517 Corr1DHist_Reduced->SetLineColor(kBlue);
1518
1519 for (int j = 0; j < size; ++j)
1520 {
1521 Corr1DHist_Reduced->GetXaxis()->SetBinLabel(j+1, NameCorrOfInterest[i][j].c_str());
1522 Corr1DHist_Reduced->SetBinContent(j+1, CorrOfInterest[i][j]);
1523 }
1524 Corr1DHist_Reduced->GetXaxis()->LabelsOption("v");
1525
1526 Corr1DHist_Reduced->SetMaximum(+1.);
1527 Corr1DHist_Reduced->SetMinimum(-1.);
1528 Corr1DHist_Reduced->Draw();
1529
1530 Posterior->Write(Form("%s_Red", Corr1DHist_Reduced->GetTitle()));
1531 if(printToPDF) Posterior->Print(CanvasName);
1532 }
1533
1534 CorrDir->Close();
1535 delete CorrDir;
1536 OutputFile->cd();
1537
1538 SetMargins(Posterior, Margins);
1539 gStyle->SetOptTitle(OptTitle);
1540}
size
int size
Definition: GetPenaltyTerm.cpp:38
MCMCProcessor::GetMargins
std::vector< double > GetMargins(const std::unique_ptr< TCanvas > &Canv) const
Get TCanvas margins, to be able to reset them if particular function need different margins.
Definition: MCMCProcessor.cpp:4158
MCMCProcessor::SetMargins
void SetMargins(std::unique_ptr< TCanvas > &Canv, const std::vector< double > &margins)
Set TCanvas margins to specified values.
Definition: MCMCProcessor.cpp:4165
MCMCProcessor::IamVaried
std::vector< bool > IamVaried
Is the ith parameter varied.
Definition: MCMCProcessor.h:399

◆ DrawCovariance()

void MCMCProcessor::DrawCovariance ( )

Draw the post-fit covariances.

Definition at line 1317 of file MCMCProcessor.cpp.

1317 {
1318// *********************
1319 const double RightMargin = Posterior->GetRightMargin();
1320 Posterior->SetRightMargin(0.15);
1321
1322 // The Covariance matrix from the fit
1323 std::unique_ptr<TH2D> hCov = std::make_unique<TH2D>("hCov", "hCov", nDraw, 0, nDraw, nDraw, 0, nDraw);
1324 hCov->GetZaxis()->SetTitle("Covariance");
1325 // The Covariance matrix square root, with correct sign
1326 std::unique_ptr<TH2D> hCovSq = std::make_unique<TH2D>("hCovSq", "hCovSq", nDraw, 0, nDraw, nDraw, 0, nDraw);
1327 hCovSq->GetZaxis()->SetTitle("Covariance");
1328 // The Correlation
1329 std::unique_ptr<TH2D> hCorr = std::make_unique<TH2D>("hCorr", "hCorr", nDraw, 0, nDraw, nDraw, 0, nDraw);
1330 hCorr->GetZaxis()->SetTitle("Correlation");
1331 hCorr->SetMinimum(-1);
1332 hCorr->SetMaximum(1);
1333 hCov->GetXaxis()->SetLabelSize(0.015);
1334 hCov->GetYaxis()->SetLabelSize(0.015);
1335 hCovSq->GetXaxis()->SetLabelSize(0.015);
1336 hCovSq->GetYaxis()->SetLabelSize(0.015);
1337 hCorr->GetXaxis()->SetLabelSize(0.015);
1338 hCorr->GetYaxis()->SetLabelSize(0.015);
1339
1340 // Loop over the Covariance matrix entries
1341 for (int i = 0; i < nDraw; ++i)
1342 {
1343 TString titlex = "";
1344 double nom, err;
1345 GetNthParameter(i, nom, err, titlex);
1346
1347 hCov->GetXaxis()->SetBinLabel(i+1, titlex);
1348 hCovSq->GetXaxis()->SetBinLabel(i+1, titlex);
1349 hCorr->GetXaxis()->SetBinLabel(i+1, titlex);
1350
1351 for (int j = 0; j < nDraw; ++j)
1352 {
1353 // The value of the Covariance
1354 const double cov = (*Covariance)(i,j);
1355 const double corr = (*Correlation)(i,j);
1356
1357 hCov->SetBinContent(i+1, j+1, cov);
1358 hCovSq->SetBinContent(i+1, j+1, ((cov > 0) - (cov < 0))*std::sqrt(std::fabs(cov)));
1359 hCorr->SetBinContent(i+1, j+1, corr);
1360
1361 TString titley = "";
1362 double nom_j, err_j;
1363 GetNthParameter(j, nom_j, err_j, titley);
1364
1365 hCov->GetYaxis()->SetBinLabel(j+1, titley);
1366 hCovSq->GetYaxis()->SetBinLabel(j+1, titley);
1367 hCorr->GetYaxis()->SetBinLabel(j+1, titley);
1368 }
1369 }
1370
1371 // Take away the stat box
1372 gStyle->SetOptStat(0);
1373 if(plotBinValue)gStyle->SetPaintTextFormat("4.1f"); //Precision of value in matrix element
1374 // Make pretty Correlation colors (red to blue)
1375 const int NRGBs = 5;
1376 TColor::InitializeColors();
1377 Double_t stops[NRGBs] = { 0.00, 0.25, 0.50, 0.75, 1.00 };
1378 Double_t red[NRGBs] = { 0.00, 0.25, 1.00, 1.00, 0.50 };
1379 Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00 };
1380 Double_t blue[NRGBs] = { 0.50, 1.00, 1.00, 0.25, 0.00 };
1381 TColor::CreateGradientColorTable(5, stops, red, green, blue, 255);
1382 gStyle->SetNumberContours(255);
1383
1384 // cd into the correlation directory
1385 OutputFile->cd();
1386
1387 Posterior->cd();
1388 Posterior->Clear();
1389 if(plotBinValue) hCov->Draw("colz text");
1390 else hCov->Draw("colz");
1391 if(printToPDF) Posterior->Print(CanvasName);
1392
1393 Posterior->cd();
1394 Posterior->Clear();
1395 if(plotBinValue) hCorr->Draw("colz text");
1396 else hCorr->Draw("colz");
1397 if(printToPDF) Posterior->Print(CanvasName);
1398
1399 hCov->Write("Covariance_plot");
1400 hCovSq->Write("Covariance_sq_plot");
1401 hCorr->Write("Correlation_plot");
1402
1403 //Back to normal
1404 Posterior->SetRightMargin(RightMargin);
1405 DrawCorrelations1D();
1406}
MCMCProcessor::DrawCorrelations1D
void DrawCorrelations1D()
Draw 1D correlations which might be more helpful than looking at huge 2D Corr matrix.
Definition: MCMCProcessor.cpp:1410

◆ DrawPostfit()

void MCMCProcessor::DrawPostfit ( )

Draw the post-fit comparisons.

Todo:
this need revision

Definition at line 420 of file MCMCProcessor.cpp.

420 {
421// *******************
422 if (OutputFile == nullptr) MakeOutputFile();
423
424 // Make the prefit plot
425 std::unique_ptr<TH1D> prefit = MakePrefit();
426
427 prefit->GetXaxis()->SetTitle("");
428 // cd into the output file
429 OutputFile->cd();
430
431 std::string CutPosterior1D = "";
432 if(Posterior1DCut != "")
433 {
434 CutPosterior1D = StepCut +" && " + Posterior1DCut;
435 }
436 else CutPosterior1D = StepCut;
437
438 // Make a TH1D of the central values and the errors
439 TH1D *paramPlot = new TH1D("paramPlot", "paramPlot", nDraw, 0, nDraw);
440 paramPlot->SetName("mach3params");
441 paramPlot->SetTitle(CutPosterior1D.c_str());
442 paramPlot->SetFillStyle(3001);
443 paramPlot->SetFillColor(kBlue-1);
444 paramPlot->SetMarkerColor(paramPlot->GetFillColor());
445 paramPlot->SetMarkerStyle(20);
446 paramPlot->SetLineColor(paramPlot->GetFillColor());
447 paramPlot->SetMarkerSize(prefit->GetMarkerSize());
448 paramPlot->GetXaxis()->SetTitle("");
449
450 // Same but with Gaussian output
451 TH1D *paramPlot_Gauss = static_cast<TH1D*>(paramPlot->Clone());
452 paramPlot_Gauss->SetMarkerColor(kOrange-5);
453 paramPlot_Gauss->SetMarkerStyle(23);
454 paramPlot_Gauss->SetLineWidth(2);
455 paramPlot_Gauss->SetMarkerSize((prefit->GetMarkerSize())*0.75);
456 paramPlot_Gauss->SetFillColor(paramPlot_Gauss->GetMarkerColor());
457 paramPlot_Gauss->SetFillStyle(3244);
458 paramPlot_Gauss->SetLineColor(paramPlot_Gauss->GetMarkerColor());
459 paramPlot_Gauss->GetXaxis()->SetTitle("");
460
461 // Same but with Gaussian output
462 TH1D *paramPlot_HPD = static_cast<TH1D*>(paramPlot->Clone());
463 paramPlot_HPD->SetMarkerColor(kBlack);
464 paramPlot_HPD->SetMarkerStyle(25);
465 paramPlot_HPD->SetLineWidth(2);
466 paramPlot_HPD->SetMarkerSize((prefit->GetMarkerSize())*0.5);
467 paramPlot_HPD->SetFillColor(0);
468 paramPlot_HPD->SetFillStyle(0);
469 paramPlot_HPD->SetLineColor(paramPlot_HPD->GetMarkerColor());
470 paramPlot_HPD->GetXaxis()->SetTitle("");
471
472 // Set labels and data
473 for (int i = 0; i < nDraw; ++i)
474 {
475 //Those keep which parameter type we run currently and realtive number
476 int ParamEnu = ParamType[i];
477 int ParamNo = i - ParamTypeStartPos[ParameterEnum(ParamEnu)];
478
479 //KS: Slightly hacky way to get relative to prior or nominal as this is convention we use
480 //This only applies for xsec for other systematic types doesn't matter
481 double CentralValueTemp = 0;
482 double Central, Central_gauss, Central_HPD;
483 double Err, Err_Gauss, Err_HPD;
484
485 if(plotRelativeToPrior)
486 {
487 CentralValueTemp = ParamCentral[ParamEnu][ParamNo];
488 // Normalise the prior relative the nominal/prior, just the way we get our fit results in MaCh3
489 if ( CentralValueTemp != 0)
490 {
491 Central = (*Means)(i) / CentralValueTemp;
492 Err = (*Errors)(i) / CentralValueTemp;
493
494 Central_gauss = (*Means_Gauss)(i) / CentralValueTemp;
495 Err_Gauss = (*Errors_Gauss)(i) / CentralValueTemp;
496
497 Central_HPD = (*Means_HPD)(i) / CentralValueTemp;
498 Err_HPD = (*Errors_HPD)(i) / CentralValueTemp;
499 }
500 else {
501 Central = 1+(*Means)(i);
502 Err = (*Errors)(i);
503
504 Central_gauss = 1+(*Means_Gauss)(i);
505 Err_Gauss = (*Errors_Gauss)(i);
506
507 Central_HPD = 1+(*Means_HPD)(i) ;
508 Err_HPD = (*Errors_HPD)(i);
509 }
510 }
511 //KS: Just get value of each parameter without dividing by prior
512 else
513 {
514 Central = (*Means)(i);
515 Err = (*Errors)(i);
516
517 Central_gauss = (*Means_Gauss)(i);
518 Err_Gauss = (*Errors_Gauss)(i);
519
520 Central_HPD = (*Means_HPD)(i) ;
521 Err_HPD = (*Errors_HPD)(i);
522 }
523
524 paramPlot->SetBinContent(i+1, Central);
525 paramPlot->SetBinError(i+1, Err);
526
527 paramPlot_Gauss->SetBinContent(i+1, Central_gauss);
528 paramPlot_Gauss->SetBinError(i+1, Err_Gauss);
529
530 paramPlot_HPD->SetBinContent(i+1, Central_HPD);
531 paramPlot_HPD->SetBinError(i+1, Err_HPD);
532
533 paramPlot->GetXaxis()->SetBinLabel(i+1, prefit->GetXaxis()->GetBinLabel(i+1));
534 paramPlot_Gauss->GetXaxis()->SetBinLabel(i+1, prefit->GetXaxis()->GetBinLabel(i+1));
535 paramPlot_HPD->GetXaxis()->SetBinLabel(i+1, prefit->GetXaxis()->GetBinLabel(i+1));
536 }
537 prefit->GetXaxis()->LabelsOption("v");
538 paramPlot->GetXaxis()->LabelsOption("v");\
539 paramPlot_Gauss->GetXaxis()->LabelsOption("v");
540 paramPlot_HPD->GetXaxis()->LabelsOption("v");
541
542 // Make a TLegend
543 auto CompLeg = std::make_unique<TLegend>(0.33, 0.73, 0.76, 0.95);
544 CompLeg->AddEntry(prefit.get(), "Prefit", "fp");
545 CompLeg->AddEntry(paramPlot, "Postfit PDF", "fp");
546 CompLeg->AddEntry(paramPlot_Gauss, "Postfit Gauss", "fp");
547 CompLeg->AddEntry(paramPlot_HPD, "Postfit HPD", "lfep");
548 CompLeg->SetFillColor(0);
549 CompLeg->SetFillStyle(0);
550 CompLeg->SetLineWidth(0);
551 CompLeg->SetLineStyle(0);
552 CompLeg->SetBorderSize(0);
553
554 const std::vector<double> Margins = GetMargins(Posterior);
555 Posterior->SetBottomMargin(0.2);
556
557 OutputFile->cd();
558
559 //KS: Plot Xsec and Flux
561 if (nParam[kXSecPar] > 0)
562 {
563 const int Start = ParamTypeStartPos[kXSecPar];
564 const int nFlux = GetGroup("Flux");
565 // Plot the xsec parameters (0 to ~nXsec-nFlux) nXsec == xsec + flux, quite confusing I know
566 // Have already looked through the branches earlier
567 if(plotRelativeToPrior) prefit->GetYaxis()->SetTitle("Variation rel. prior");
568 else prefit->GetYaxis()->SetTitle("Parameter Value");
569 prefit->GetYaxis()->SetRangeUser(-2.5, 2.5);
570
571 prefit->GetXaxis()->SetRangeUser(Start, Start + nParam[kXSecPar]-nFlux);
572 paramPlot->GetXaxis()->SetRangeUser(Start, Start + nParam[kXSecPar]-nFlux);
573 paramPlot_Gauss->GetXaxis()->SetRangeUser(Start, Start+ nParam[kXSecPar]-nFlux);
574 paramPlot_HPD->GetXaxis()->SetRangeUser(Start, Start + nParam[kXSecPar]-nFlux);
575
576 // Write the individual ones
577 prefit->Write("param_xsec_prefit");
578 paramPlot->Write("param_xsec");
579 paramPlot_Gauss->Write("param_xsec_gaus");
580 paramPlot_HPD->Write("param_xsec_HPD");
581
582 // And the combined
583 prefit->Draw("e2");
584 paramPlot->Draw("e2, same");
585 paramPlot_Gauss->Draw("e2, same");
586 paramPlot_HPD->Draw("e1, same");
587 CompLeg->Draw("same");
588 Posterior->Write("param_xsec_canv");
589 if(printToPDF) Posterior->Print(CanvasName);
590 Posterior->Clear();
591
592 OutputFile->cd();
593 // Plot the flux parameters (nXSec to nxsec+nflux) if enabled
594 // Have already looked through the branches earlier
595 prefit->GetYaxis()->SetRangeUser(0.7, 1.3);
596 paramPlot->GetYaxis()->SetRangeUser(0.7, 1.3);
597 paramPlot_Gauss->GetYaxis()->SetRangeUser(0.7, 1.3);
598 paramPlot_HPD->GetYaxis()->SetRangeUser(0.7, 1.3);
599
600 prefit->GetXaxis()->SetRangeUser(Start + nParam[kXSecPar]-nFlux, Start + nParam[kXSecPar]);
601 paramPlot->GetXaxis()->SetRangeUser(Start + nParam[kXSecPar]-nFlux, Start + nParam[kXSecPar]);
602 paramPlot_Gauss->GetXaxis()->SetRangeUser(Start + nParam[kXSecPar]-nFlux, Start + nParam[kXSecPar]);
603 paramPlot_HPD->GetXaxis()->SetRangeUser(Start + nParam[kXSecPar]-nFlux, Start + nParam[kXSecPar]);
604
605 prefit->Write("param_flux_prefit");
606 paramPlot->Write("param_flux");
607 paramPlot_Gauss->Write("param_flux_gaus");
608 paramPlot_HPD->Write("param_flux_HPD");
609
610 prefit->Draw("e2");
611 paramPlot->Draw("e2, same");
612 paramPlot_Gauss->Draw("e1, same");
613 paramPlot_HPD->Draw("e1, same");
614 CompLeg->Draw("same");
615 Posterior->Write("param_flux_canv");
616 if(printToPDF) Posterior->Print(CanvasName);
617 Posterior->Clear();
618 }
619 if(nParam[kNDPar] > 0)
620 {
621 int Start = ParamTypeStartPos[kNDPar];
622 int NDbinCounter = Start;
623 //KS: Make prefit postfit for each ND sample, having all of them at the same plot is unreadable
624 for(unsigned int i = 0; i < NDSamplesNames.size(); i++ )
625 {
626 std::string NDname = NDSamplesNames[i];
627 NDbinCounter += NDSamplesBins[i];
628 OutputFile->cd();
629 prefit->GetYaxis()->SetTitle(("Variation for "+NDname).c_str());
630 prefit->GetYaxis()->SetRangeUser(0.6, 1.4);
631 prefit->GetXaxis()->SetRangeUser(Start, NDbinCounter);
632
633 paramPlot->GetYaxis()->SetTitle(("Variation for "+NDname).c_str());
634 paramPlot->GetYaxis()->SetRangeUser(0.6, 1.4);
635 paramPlot->GetXaxis()->SetRangeUser(Start, NDbinCounter);
636 paramPlot->SetTitle(CutPosterior1D.c_str());
637
638 paramPlot_Gauss->GetYaxis()->SetTitle(("Variation for "+NDname).c_str());
639 paramPlot_Gauss->GetYaxis()->SetRangeUser(0.6, 1.4);
640 paramPlot_Gauss->GetXaxis()->SetRangeUser(Start, NDbinCounter);
641 paramPlot_Gauss->SetTitle(CutPosterior1D.c_str());
642
643 paramPlot_HPD->GetYaxis()->SetTitle(("Variation for "+NDname).c_str());
644 paramPlot_HPD->GetYaxis()->SetRangeUser(0.6, 1.4);
645 paramPlot_HPD->GetXaxis()->SetRangeUser(Start, NDbinCounter);
646 paramPlot_HPD->SetTitle(CutPosterior1D.c_str());
647
648 prefit->Write(("param_"+NDname+"_prefit").c_str());
649 paramPlot->Write(("param_"+NDname).c_str());
650 paramPlot_Gauss->Write(("param_"+NDname+"_gaus").c_str());
651 paramPlot_HPD->Write(("param_"+NDname+"_HPD").c_str());
652
653 prefit->Draw("e2");
654 paramPlot->Draw("e2, same");
655 paramPlot_Gauss->Draw("e1, same");
656 paramPlot_HPD->Draw("e1, same");
657 CompLeg->Draw("same");
658 Posterior->Write(("param_"+NDname+"_canv").c_str());
659 if(printToPDF) Posterior->Print(CanvasName);
660 Posterior->Clear();
661 Start += NDSamplesBins[i];
662 }
663 }
664 delete paramPlot;
665 delete paramPlot_Gauss;
666 delete paramPlot_HPD;
667
668 //KS: Return Margin to default one
669 SetMargins(Posterior, Margins);
670}
ParameterEnum
ParameterEnum
KS: Enum for different covariance classes.
Definition: MCMCProcessor.h:52
kNDPar
@ kNDPar
Definition: MCMCProcessor.h:54
kXSecPar
@ kXSecPar
Definition: MCMCProcessor.h:53
MCMCProcessor::GetGroup
int GetGroup(const std::string &name) const
Number of params from a given group, for example flux.
Definition: MCMCProcessor.cpp:4114
MCMCProcessor::MakePrefit
std::unique_ptr< TH1D > MakePrefit()
Prepare prefit histogram for parameter overlay plot.
Definition: MCMCProcessor.cpp:2100
MCMCProcessor::MakeOutputFile
void MakeOutputFile()
prepare output root file and canvas to which we will save EVERYTHING
Definition: MCMCProcessor.cpp:201
MCMCProcessor::ParamType
std::vector< ParameterEnum > ParamType
Make an enum for which class this parameter belongs to so we don't have to keep string comparing.
Definition: MCMCProcessor.h:412
MCMCProcessor::NDSamplesNames
std::vector< std::string > NDSamplesNames
Definition: MCMCProcessor.h:451
MCMCProcessor::NDSamplesBins
std::vector< int > NDSamplesBins
Definition: MCMCProcessor.h:450

◆ FindInputFiles()

void MCMCProcessor::FindInputFiles ( )
inlineprotected

Read the output MCMC file and find what inputs were used.

Definition at line 2180 of file MCMCProcessor.cpp.

2180 {
2181// **************************
2182 // Now read the MCMC file
2183 TFile *TempFile = new TFile(MCMCFile.c_str(), "open");
2184 TDirectory* CovarianceFolder = TempFile->Get<TDirectory>("CovarianceFolder");
2185
2186 // Get the settings for the MCMC
2187 TMacro *Config = TempFile->Get<TMacro>("MaCh3_Config");
2188
2189 if (Config == nullptr) {
2190 MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", MCMCFile);
2191 TempFile->ls();
2192 throw MaCh3Exception(__FILE__ , __LINE__ );
2193 }
2194 //KS:Most inputs are in ${MACH3}/inputs/blarb.root
2195 if (std::getenv("MACH3") != nullptr) {
2196 MACH3LOG_INFO("Found MACH3 environment variable: {}", std::getenv("MACH3"));
2197 }
2198
2199 MACH3LOG_INFO("Loading YAML config from MCMC chain");
2200
2201 YAML::Node Settings = TMacroToYAML(*Config);
2202
2203 bool InputNotFound = false;
2204 //CW: Get the xsec Covariance matrix
2205 CovPos[kXSecPar] = GetFromManager<std::vector<std::string>>(Settings["General"]["Systematics"]["XsecCovFile"], {"none"});
2206 if(CovPos[kXSecPar].back() == "none")
2207 {
2208 MACH3LOG_WARN("Couldn't find XsecCov branch in output");
2209 InputNotFound = true;
2210 }
2211
2212 TMacro *XsecConfig = M3::GetConfigMacroFromChain(CovarianceFolder);
2213 if (XsecConfig == nullptr) {
2214 MACH3LOG_WARN("Didn't find Config_xsec_cov tree in MCMC file! {}", MCMCFile);
2215 } else {
2216 CovConfig[kXSecPar] = TMacroToYAML(*XsecConfig);
2217 }
2218 if(InputNotFound) MaCh3Utils::PrintConfig(Settings);
2219
2220 if (const char * mach3_env = std::getenv("MACH3"))
2221 {
2222 for(size_t i = 0; i < CovPos[kXSecPar].size(); i++)
2223 CovPos[kXSecPar][i].insert(0, std::string(mach3_env)+"/");
2224 }
2225
2226 // Delete the TTrees and the input file handle since we've now got the settings we need
2227 delete Config;
2228 delete XsecConfig;
2229
2230 // Delete the MCMCFile pointer we're reading
2231 CovarianceFolder->Close();
2232 delete CovarianceFolder;
2233 TempFile->Close();
2234 delete TempFile;
2235}
MACH3LOG_WARN
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
TMacroToYAML
YAML::Node TMacroToYAML(const TMacro &macro)
KS: Convert a ROOT TMacro object to a YAML node.
Definition: YamlHelper.h:146
M3::GetConfigMacroFromChain
TMacro * GetConfigMacroFromChain(TDirectory *CovarianceFolder)
KS: We store configuration macros inside the chain. In the past, multiple configs were stored,...
Definition: ParameterHandlerUtils.h:366
MaCh3Utils::PrintConfig
void PrintConfig(const YAML::Node &node)
KS: Print Yaml config using logger.
Definition: Monitor.cpp:294

◆ FindInputFilesLegacy()

void MCMCProcessor::FindInputFilesLegacy ( )
inlineprotected
Warning
This will no longer be supported in future

Definition at line 2239 of file MCMCProcessor.cpp.

2239 {
2240// **************************
2241 // Now read the MCMC file
2242 TFile *TempFile = new TFile(MCMCFile.c_str(), "open");
2243
2244 // Get the settings for the MCMC
2245 TMacro *Config = TempFile->Get<TMacro>("MaCh3_Config");
2246
2247 if (Config == nullptr) {
2248 MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", MCMCFile);
2249 TempFile->ls();
2250 throw MaCh3Exception(__FILE__ , __LINE__ );
2251 }
2252 YAML::Node Settings = TMacroToYAML(*Config);
2253
2254 //CW: And the ND Covariance matrix
2255 CovPos[kNDPar].push_back(GetFromManager<std::string>(Settings["General"]["Systematics"]["NDCovFile"], "none"));
2256 if(CovPos[kNDPar].back() == "none") {
2257 MACH3LOG_WARN("Couldn't find NDCov (legacy) branch in output");
2258 }
2259
2260 //CW: And the FD Covariance matrix
2261 CovPos[kFDDetPar].push_back(GetFromManager<std::string>(Settings["General"]["Systematics"]["FDCovFile"], "none"));
2262 if(CovPos[kFDDetPar].back() == "none") {
2263 MACH3LOG_WARN("Couldn't find FDCov (legacy) branch in output");
2264 }
2265
2266 if (const char * mach3_env = std::getenv("MACH3"))
2267 {
2268 for(size_t i = 0; i < CovPos[kNDPar].size(); i++)
2269 CovPos[kNDPar][i].insert(0, std::string(mach3_env)+"/");
2270
2271 for(size_t i = 0; i < CovPos[kFDDetPar].size(); i++)
2272 CovPos[kFDDetPar][i].insert(0, std::string(mach3_env)+"/");
2273 }
2274 TempFile->Close();
2275 delete TempFile;
2276}
kFDDetPar
@ kFDDetPar
Definition: MCMCProcessor.h:55

◆ GetBayesFactor()

void MCMCProcessor::GetBayesFactor ( const std::vector< std::string > &  ParName,
const std::vector< std::vector< double > > &  Model1Bounds,
const std::vector< std::vector< double > > &  Model2Bounds,
const std::vector< std::vector< std::string > > &  ModelNames 
)

Calculate Bayes factor for vector of params, and model boundaries.

Parameters
ParNameVector with parameter names for which we calculate Bayes factor
Model1BoundsLower and upper bound for hypothesis 1. Within this bound we calculate integral used later for Bayes Factor
Model2BoundsLower and upper bound for hypothesis 2. Within this bound we calculate integral used later for Bayes Factor
ModelNamesNames for hypothesis 1 and 2

Definition at line 2545 of file MCMCProcessor.cpp.

2548 {
2549// **************************
2550 if(hpost[0] == nullptr) MakePostfit();
2551
2552 MACH3LOG_INFO("Calculating Bayes Factor");
2553 if((ParNames.size() != Model1Bounds.size()) || (Model2Bounds.size() != Model1Bounds.size()) || (Model2Bounds.size() != ModelNames.size()))
2554 {
2555 MACH3LOG_ERROR("Size doesn't match");
2556 throw MaCh3Exception(__FILE__ , __LINE__ );
2557 }
2558 for(unsigned int k = 0; k < ParNames.size(); ++k)
2559 {
2560 //KS: First we need to find parameter number based on name
2561 int ParamNo = GetParamIndexFromName(ParNames[k]);
2562 if(ParamNo == _UNDEF_)
2563 {
2564 MACH3LOG_WARN("Couldn't find param {}. Will not calculate Bayes Factor", ParNames[k]);
2565 continue;
2566 }
2567
2568 const double M1_min = Model1Bounds[k][0];
2569 const double M2_min = Model2Bounds[k][0];
2570 const double M1_max = Model1Bounds[k][1];
2571 const double M2_max = Model2Bounds[k][1];
2572
2573 long double IntegralMode1 = hpost[ParamNo]->Integral(hpost[ParamNo]->FindFixBin(M1_min), hpost[ParamNo]->FindFixBin(M1_max));
2574 long double IntegralMode2 = hpost[ParamNo]->Integral(hpost[ParamNo]->FindFixBin(M2_min), hpost[ParamNo]->FindFixBin(M2_max));
2575
2576 double BayesFactor = 0.;
2577 std::string Name = "";
2578 //KS: Calc Bayes Factor
2579 //If M1 is more likely
2580 if(IntegralMode1 >= IntegralMode2)
2581 {
2582 BayesFactor = IntegralMode1/IntegralMode2;
2583 Name = "\\mathfrak{B}(" + ModelNames[k][0]+ "/" + ModelNames[k][1] + ") = " + std::to_string(BayesFactor);
2584 }
2585 else //If M2 is more likely
2586 {
2587 BayesFactor = IntegralMode2/IntegralMode1;
2588 Name = "\\mathfrak{B}(" + ModelNames[k][1]+ "/" + ModelNames[k][0] + ") = " + std::to_string(BayesFactor);
2589 }
2590 std::string JeffreysScale = GetJeffreysScale(BayesFactor);
2591 std::string DunneKabothScale = GetDunneKaboth(BayesFactor);
2592
2593 MACH3LOG_INFO("{} for {}", Name, ParNames[k]);
2594 MACH3LOG_INFO("Following Jeffreys Scale = {}", JeffreysScale);
2595 MACH3LOG_INFO("Following Dunne-Kaboth Scale = {}", DunneKabothScale);
2596 MACH3LOG_INFO("");
2597 }
2598}
GetJeffreysScale
std::string GetJeffreysScale(const double BayesFactor)
KS: Following H. Jeffreys .
Definition: StatisticalUtils.cpp:5
GetDunneKaboth
std::string GetDunneKaboth(const double BayesFactor)
KS: Based on Table 1 in https://www.t2k.org/docs/technotes/435.
Definition: StatisticalUtils.cpp:20
MCMCProcessor::GetParamIndexFromName
int GetParamIndexFromName(const std::string &Name)
Get parameter number based on name.
Definition: MCMCProcessor.cpp:2438
MCMCProcessor::MakePostfit
void MakePostfit()
Make 1D projection for each parameter and prepare structure.
Definition: MCMCProcessor.cpp:235

◆ GetBranchNames()

const std::vector< TString > & MCMCProcessor::GetBranchNames ( ) const
inline

Get the vector of branch names from root file.

Definition at line 227 of file MCMCProcessor.h.

227{ return BranchNames;};

◆ GetCovariance()

void MCMCProcessor::GetCovariance ( TMatrixDSym *&  Cov,
TMatrixDSym *&  Corr 
)

Get the post-fit covariances and correlations.

Parameters
CovCovariance matrix
CorrCorrelation matrix

Definition at line 192 of file MCMCProcessor.cpp.

192 {
193// ***************
194 if (CacheMCMC) MakeCovariance_MP();
195 else MakeCovariance();
196 Cov = static_cast<TMatrixDSym*>(Covariance->Clone());
197 Corr = static_cast<TMatrixDSym*>(Correlation->Clone());
198}
MCMCProcessor::MakeCovariance_MP
void MakeCovariance_MP(const bool Mute=false)
Calculate covariance by making 2D projection of each combination of parameters using multithreading.
Definition: MCMCProcessor.cpp:1141
MCMCProcessor::MakeCovariance
void MakeCovariance()
Calculate covariance by making 2D projection of each combination of parameters.
Definition: MCMCProcessor.cpp:925

◆ GetFDCov()

std::string MCMCProcessor::GetFDCov ( ) const
inline

Definition at line 215 of file MCMCProcessor.h.

215{ return CovPos[kFDDetPar].back(); };

◆ GetGroup()

int MCMCProcessor::GetGroup ( const std::string &  name) const

Number of params from a given group, for example flux.

Definition at line 4114 of file MCMCProcessor.cpp.

4114 {
4115// **************************
4116 // Lambda to compare strings case-insensitively
4117 auto caseInsensitiveCompare = [](const std::string& a, const std::string& b) {
4118 return std::equal(a.begin(), a.end(), b.begin(), b.end(),
4119 [](char c1, char c2) { return std::tolower(c1) == std::tolower(c2); });
4120 };
4121 int numerator = 0;
4122 for (const auto& groupName : ParameterGroup) {
4123 if (caseInsensitiveCompare(groupName, name)) {
4124 numerator++;
4125 }
4126 }
4127 return numerator;
4128}
MCMCProcessor::ParameterGroup
std::vector< std::string > ParameterGroup
Definition: MCMCProcessor.h:417

◆ GetHpost()

TH1D * MCMCProcessor::GetHpost ( const int  i)
inline

Get 1D posterior for a given parameter.

Parameters
iparameter index

Definition at line 202 of file MCMCProcessor.h.

202{ return hpost[i]; };

◆ GetHpost2D()

TH2D * MCMCProcessor::GetHpost2D ( const int  i,
const int  j 
)
inline

Get 2D posterior for a given parameter combination.

Parameters
iparameter index X
jparameter index Y

Definition at line 206 of file MCMCProcessor.h.

206{ return hpost2D[i][j]; };

◆ GetMargins()

std::vector< double > MCMCProcessor::GetMargins ( const std::unique_ptr< TCanvas > &  Canv) const
protected

Get TCanvas margins, to be able to reset them if particular function need different margins.

Definition at line 4158 of file MCMCProcessor.cpp.

4158 {
4159// **************************
4160 return std::vector<double>{Canv->GetTopMargin(), Canv->GetBottomMargin(),
4161 Canv->GetLeftMargin(), Canv->GetRightMargin()};
4162}

◆ GetNDCov()

std::string MCMCProcessor::GetNDCov ( ) const
inline

Definition at line 214 of file MCMCProcessor.h.

214{ return CovPos[kNDPar].back(); };

◆ GetnEntries()

Long64_t MCMCProcessor::GetnEntries ( )
inline

Get Number of entries that Chain has, for merged chains will not be the same Nsteps.

Definition at line 233 of file MCMCProcessor.h.

233{return nEntries;};

◆ GetNFD()

int MCMCProcessor::GetNFD ( )
inline

Definition at line 196 of file MCMCProcessor.h.

196{ return nParam[kFDDetPar]; };

◆ GetNND()

int MCMCProcessor::GetNND ( )
inline

Definition at line 195 of file MCMCProcessor.h.

195{ return nParam[kNDPar]; };

◆ GetNParams()

int MCMCProcessor::GetNParams ( )
inline

Get total number of used parameters.

Definition at line 193 of file MCMCProcessor.h.

193{ return nDraw; };

◆ GetnSteps()

Long64_t MCMCProcessor::GetnSteps ( )
inline

Get Number of Steps that Chain has, for merged chains will not be the same nEntries.

Definition at line 235 of file MCMCProcessor.h.

235{return nSteps;};

◆ GetNthParameter()

void MCMCProcessor::GetNthParameter ( const int  param,
double &  Prior,
double &  PriorError,
TString &  Title 
) const

Get properties of parameter by passing it number.

Definition at line 2425 of file MCMCProcessor.cpp.

2425 {
2426// **************************
2427 ParameterEnum ParType = ParamType[param];
2428 int ParamNo = _UNDEF_;
2429 ParamNo = param - ParamTypeStartPos[ParType];
2430
2431 Prior = ParamCentral[ParType][ParamNo];
2432 PriorError = ParamErrors[ParType][ParamNo];
2433 Title = ParamNames[ParType][ParamNo];
2434}

◆ GetNXSec()

int MCMCProcessor::GetNXSec ( )
inline

Definition at line 194 of file MCMCProcessor.h.

194{ return nParam[kXSecPar]; };

◆ GetParamIndexFromName()

int MCMCProcessor::GetParamIndexFromName ( const std::string &  Name)

Get parameter number based on name.

Definition at line 2438 of file MCMCProcessor.cpp.

2438 {
2439// **************************
2440 int ParamNo = _UNDEF_;
2441 for (int i = 0; i < nDraw; ++i)
2442 {
2443 TString Title = "";
2444 double Prior = 1.0, PriorError = 1.0;
2445 GetNthParameter(i, Prior, PriorError, Title);
2446
2447 if(Name == Title)
2448 {
2449 ParamNo = i;
2450 break;
2451 }
2452 }
2453 return ParamNo;
2454}

◆ GetPolarPlot()

void MCMCProcessor::GetPolarPlot ( const std::vector< std::string > &  ParNames)

Make funny polar plot.

Parameters
ParNamesVector with parameter names for which Polar Plot will be made

Definition at line 2476 of file MCMCProcessor.cpp.

2476 {
2477// **************************
2478 if(hpost[0] == nullptr) MakePostfit();
2479
2480 std::vector<double> Margins = GetMargins(Posterior);
2481
2482 Posterior->SetTopMargin(0.1);
2483 Posterior->SetBottomMargin(0.1);
2484 Posterior->SetLeftMargin(0.1);
2485 Posterior->SetRightMargin(0.1);
2486 Posterior->Update();
2487
2488 MACH3LOG_INFO("Calculating Polar Plot");
2489 TDirectory *PolarDir = OutputFile->mkdir("PolarDir");
2490 PolarDir->cd();
2491
2492 for(unsigned int k = 0; k < ParNames.size(); ++k)
2493 {
2494 //KS: First we need to find parameter number based on name
2495 int ParamNo = GetParamIndexFromName(ParNames[k]);
2496 if(ParamNo == _UNDEF_)
2497 {
2498 MACH3LOG_WARN("Couldn't find param {}. Will not calculate Polar Plot", ParNames[k]);
2499 continue;
2500 }
2501
2502 TString Title = "";
2503 double Prior = 1.0, PriorError = 1.0;
2504 GetNthParameter(ParamNo, Prior, PriorError, Title);
2505
2506 std::vector<double> x_val(nBins);
2507 std::vector<double> y_val(nBins);
2508
2509 double xmin = 0;
2510 double xmax = 2*TMath::Pi();
2511
2512 double Integral = hpost[ParamNo]->Integral();
2513 for (Int_t ipt = 0; ipt < nBins; ipt++)
2514 {
2515 x_val[ipt] = ipt*(xmax-xmin)/nBins+xmin;
2516 y_val[ipt] = hpost[ParamNo]->GetBinContent(ipt+1)/Integral;
2517 }
2518
2519 auto PolarGraph = std::make_unique<TGraphPolar>(nBins, x_val.data(), y_val.data());
2520 PolarGraph->SetLineWidth(2);
2521 PolarGraph->SetFillStyle(3001);
2522 PolarGraph->SetLineColor(kRed);
2523 PolarGraph->SetFillColor(kRed);
2524 PolarGraph->Draw("AFL");
2525
2526 auto Text = std::make_unique<TText>(0.6, 0.1, Title);
2527 Text->SetTextSize(0.04);
2528 Text->SetNDC(true);
2529 Text->Draw("");
2530
2531 Posterior->Print(CanvasName);
2532 Posterior->Write(Title);
2533 } //End loop over parameters
2534
2535 PolarDir->Close();
2536 delete PolarDir;
2537
2538 OutputFile->cd();
2539
2540 SetMargins(Posterior, Margins);
2541}

◆ GetPostfit()

void MCMCProcessor::GetPostfit ( TVectorD *&  Central,
TVectorD *&  Errors,
TVectorD *&  Central_Gauss,
TVectorD *&  Errors_Gauss,
TVectorD *&  Peaks 
)

Get the post-fit results (arithmetic and Gaussian)

Definition at line 159 of file MCMCProcessor.cpp.

159 {
160// ***************
161 // Make the post fit
162 MakePostfit();
163
164 // We now have the private members
165 Central_PDF = Means;
166 Errors_PDF = Errors;
167 Central_G = Means_Gauss;
168 Errors_G = Errors_Gauss;
169 Peak_Values = Means_HPD;
170}

◆ GetPostfit_Ind()

void MCMCProcessor::GetPostfit_Ind ( TVectorD *&  Central,
TVectorD *&  Errors,
TVectorD *&  Peaks,
ParameterEnum  kParam 
)

Or the individual post-fits.

Definition at line 174 of file MCMCProcessor.cpp.

174 {
175// ***************
176 // Make the post fit
177 MakePostfit();
178
179 // Loop over the loaded param types
180 const int ParamTypeSize = int(ParamType.size());
181 int ParamNumber = 0;
182 for (int i = 0; i < ParamTypeSize; ++i) {
183 if (ParamType[i] != kParam) continue;
184 (*PDF_Central)(ParamNumber) = (*Means)(i);
185 (*PDF_Errors)(ParamNumber) = (*Errors)(i);
186 (*Peak_Values)(ParamNumber) = (*Means_HPD)(i);
187 ++ParamNumber;
188 }
189}

◆ GetSavageDickey()

void MCMCProcessor::GetSavageDickey ( const std::vector< std::string > &  ParName,
const std::vector< double > &  EvaluationPoint,
const std::vector< std::vector< double > > &  Bounds 
)

Calculate Bayes factor for point like hypothesis using SavageDickey.

Definition at line 2602 of file MCMCProcessor.cpp.

2604 {
2605// **************************
2606 if((ParNames.size() != EvaluationPoint.size()) || (Bounds.size() != EvaluationPoint.size()))
2607 {
2608 MACH3LOG_ERROR("Size doesn't match");
2609 throw MaCh3Exception(__FILE__ , __LINE__ );
2610 }
2611
2612 if(hpost[0] == nullptr) MakePostfit();
2613
2614 MACH3LOG_INFO("Calculating Savage Dickey");
2615 TDirectory *SavageDickeyDir = OutputFile->mkdir("SavageDickey");
2616 SavageDickeyDir->cd();
2617
2618 for(unsigned int k = 0; k < ParNames.size(); ++k)
2619 {
2620 //KS: First we need to find parameter number based on name
2621 int ParamNo = GetParamIndexFromName(ParNames[k]);
2622 if(ParamNo == _UNDEF_)
2623 {
2624 MACH3LOG_WARN("Couldn't find param {}. Will not calculate SavageDickey", ParNames[k]);
2625 continue;
2626 }
2627
2628 TString Title = "";
2629 double Prior = 1.0, PriorError = 1.0;
2630 bool FlatPrior = false;
2631 GetNthParameter(ParamNo, Prior, PriorError, Title);
2632
2633 ParameterEnum ParType = ParamType[ParamNo];
2634 int ParamTemp = ParamNo - ParamTypeStartPos[ParType];
2635 FlatPrior = ParamFlat[ParType][ParamTemp];
2636
2637 TH1D* PosteriorHist = static_cast<TH1D *>(hpost[ParamNo]->Clone(Title));
2638 RemoveFitter(PosteriorHist, "Gauss");
2639
2640 TH1D* PriorHist = nullptr;
2641 //KS: If flat prior we need to have well defined bounds otherwise Prior distribution will not make sense
2642 if(FlatPrior)
2643 {
2644 int NBins = PosteriorHist->GetNbinsX();
2645 if(Bounds[k][0] > Bounds[k][1])
2646 {
2647 MACH3LOG_ERROR("Lower bound is higher than upper bound");
2648 throw MaCh3Exception(__FILE__ , __LINE__ );
2649 }
2650 PriorHist = new TH1D("PriorHist", Title, NBins, Bounds[k][0], Bounds[k][1]);
2651
2652 double FlatProb = ( Bounds[k][1] - Bounds[k][0]) / NBins;
2653 for (int g = 0; g < NBins + 1; ++g)
2654 {
2655 PriorHist->SetBinContent(g+1, FlatProb);
2656 }
2657 }
2658 else //KS: Otherwise throw from Gaussian
2659 {
2660 PriorHist = static_cast<TH1D*>(PosteriorHist->Clone("Prior"));
2661 PriorHist->Reset("");
2662 PriorHist->Fill(0.0, 0.0);
2663
2664 auto rand = std::make_unique<TRandom3>(0);
2665 //KS: Throw nice gaussian, just need big number to have smooth distribution
2666 for(int g = 0; g < 1000000; ++g)
2667 {
2668 PriorHist->Fill(rand->Gaus(Prior, PriorError));
2669 }
2670 }
2671 // Area normalise the distributions
2672 PriorHist->Scale(1./PriorHist->Integral(), "width");
2673 PosteriorHist->Scale(1./PosteriorHist->Integral(), "width");
2674
2675 PriorHist->SetLineColor(kRed);
2676 PriorHist->SetMarkerColor(kRed);
2677 PriorHist->SetFillColorAlpha(kRed, 0.35);
2678 PriorHist->SetFillStyle(1001);
2679 PriorHist->GetXaxis()->SetTitle(Title);
2680 PriorHist->GetYaxis()->SetTitle("Posterior Probability");
2681 PriorHist->SetMaximum(PosteriorHist->GetMaximum()*1.5);
2682 PriorHist->GetYaxis()->SetLabelOffset(999);
2683 PriorHist->GetYaxis()->SetLabelSize(0);
2684 PriorHist->SetLineWidth(2);
2685 PriorHist->SetLineStyle(kSolid);
2686
2687 PosteriorHist->SetLineColor(kBlue);
2688 PosteriorHist->SetMarkerColor(kBlue);
2689 PosteriorHist->SetFillColorAlpha(kBlue, 0.35);
2690 PosteriorHist->SetFillStyle(1001);
2691
2692 PriorHist->Draw("hist");
2693 PosteriorHist->Draw("hist same");
2694
2695 double ProbPrior = PriorHist->GetBinContent(PriorHist->FindBin(EvaluationPoint[k]));
2696 //KS: In case we go so far away that prior is 0, set this to small value to avoid dividing by 0
2697 if(ProbPrior < 0) ProbPrior = 0.00001;
2698 double ProbPosterior = PosteriorHist->GetBinContent(PosteriorHist->FindBin(EvaluationPoint[k]));
2699 double SavageDickey = ProbPosterior/ProbPrior;
2700
2701 std::string DunneKabothScale = GetDunneKaboth(SavageDickey);
2702 //Get Best point
2703 std::unique_ptr<TGraph> PostPoint(new TGraph(1));
2704 PostPoint->SetPoint(0, EvaluationPoint[k], ProbPosterior);
2705 PostPoint->SetMarkerStyle(20);
2706 PostPoint->SetMarkerSize(1);
2707 PostPoint->Draw("P same");
2708
2709 std::unique_ptr<TGraph> PriorPoint(new TGraph(1));
2710 PriorPoint->SetPoint(0, EvaluationPoint[k], ProbPrior);
2711 PriorPoint->SetMarkerStyle(20);
2712 PriorPoint->SetMarkerSize(1);
2713 PriorPoint->Draw("P same");
2714
2715 auto legend = std::make_unique<TLegend>(0.12, 0.6, 0.6, 0.97);
2716 SetLegendStyle(legend.get(), 0.04);
2717 legend->AddEntry(PriorHist, "Prior", "l");
2718 legend->AddEntry(PosteriorHist, "Posterior", "l");
2719 legend->AddEntry(PostPoint.get(), Form("SavageDickey = %.2f, (%s)", SavageDickey, DunneKabothScale.c_str()),"");
2720 legend->Draw("same");
2721
2722 Posterior->Print(CanvasName);
2723 Posterior->Write(Title);
2724
2725 delete PosteriorHist;
2726 delete PriorHist;
2727 } //End loop over parameters
2728
2729 SavageDickeyDir->Close();
2730 delete SavageDickeyDir;
2731
2732 OutputFile->cd();
2733}
RemoveFitter
void RemoveFitter(TH1D *hist, const std::string &name)
KS: Remove fitted TF1 from hist to make comparison easier.
Definition: HistogramUtils.cpp:394

◆ GetViolin()

TH2D * MCMCProcessor::GetViolin ( )
inline

Get Violin plot for all parameters with posterior values.

Definition at line 208 of file MCMCProcessor.h.

208{ return hviolin.get(); };

◆ GetViolinPrior()

TH2D * MCMCProcessor::GetViolinPrior ( )
inline

Get Violin plot for all parameters with prior values.

Definition at line 210 of file MCMCProcessor.h.

210{ return hviolin_prior.get(); };

◆ GetXSecCov()

std::vector< std::string > MCMCProcessor::GetXSecCov ( ) const
inline

Definition at line 213 of file MCMCProcessor.h.

213{ return CovPos[kXSecPar]; };

◆ GewekeDiagnostic()

void MCMCProcessor::GewekeDiagnostic ( )
inlineprotected

Geweke Diagnostic based on the methods described by Fang (2014) and Karlsbakk (2011). [8] [18].

Definition at line 3898 of file MCMCProcessor.cpp.

3898 {
3899// **************************
3900 MACH3LOG_INFO("Making Geweke Diagnostic");
3901 //KS: Up refers to upper limit we check, it stays constant, in literature it is mostly 50% thus using 0.5 for threshold
3902 std::vector<double> MeanUp(nDraw, 0.0);
3903 std::vector<double> SpectralVarianceUp(nDraw, 0.0);
3904 std::vector<int> DenomCounterUp(nDraw, 0);
3905 const double Threshold = 0.5 * nSteps;
3906
3907 //KS: Select values between which you want to scan, for example 0 means 0% burn in and 1 100% burn in.
3908 constexpr double LowerThreshold = 0;
3909 constexpr double UpperThreshold = 1.0;
3910 // Tells how many intervals between thresholds we want to check
3911 constexpr int NChecks = 100;
3912 constexpr double Division = (UpperThreshold - LowerThreshold)/NChecks;
3913
3914 std::vector<std::unique_ptr<TH1D>> GewekePlots(nDraw);
3915 for (int j = 0; j < nDraw; ++j)
3916 {
3917 TString Title = "";
3918 double Prior = 1.0, PriorError = 1.0;
3919 GetNthParameter(j, Prior, PriorError, Title);
3920 std::string HistName = Form("%s_%s_Geweke", Title.Data(), BranchNames[j].Data());
3921 GewekePlots[j] = std::make_unique<TH1D>(HistName.c_str(), HistName.c_str(), NChecks, 0.0, 100 * UpperThreshold);
3922 GewekePlots[j]->GetXaxis()->SetTitle("Burn-In (%)");
3923 GewekePlots[j]->GetYaxis()->SetTitle("Geweke T score");
3924 }
3925
3926//KS: Start parallel region
3927#ifdef MULTITHREAD
3928#pragma omp parallel
3929{
3930#endif
3931 //KS: First we calculate mean and spectral variance for the upper limit, this doesn't change and in literature is most often 50%
3932 #ifdef MULTITHREAD
3933 #pragma omp for
3934 #endif
3935 for (int j = 0; j < nDraw; ++j)
3936 {
3937 for(int i = 0; i < nEntries; ++i)
3938 {
3939 if(StepNumber[i] > Threshold)
3940 {
3941 MeanUp[j] += ParStep[j][i];
3942 DenomCounterUp[j]++;
3943 }
3944 }
3945 MeanUp[j] = MeanUp[j]/DenomCounterUp[j];
3946 }
3947
3948 //KS: now Spectral variance which in this case is sample variance
3949 #ifdef MULTITHREAD
3950 #pragma omp for collapse(2)
3951 #endif
3952 for (int j = 0; j < nDraw; ++j)
3953 {
3954 for(int i = 0; i < nEntries; ++i)
3955 {
3956 if(StepNumber[i] > Threshold)
3957 {
3958 SpectralVarianceUp[j] += (ParStep[j][i] - MeanUp[j])*(ParStep[j][i] - MeanUp[j]);
3959 }
3960 }
3961 }
3962
3963 //Loop over how many intervals we calculate
3964 #ifdef MULTITHREAD
3965 #pragma omp for
3966 #endif
3967 for (int k = 1; k < NChecks+1; ++k)
3968 {
3969 //KS each thread has it's own
3970 std::vector<double> MeanDown(nDraw, 0.0);
3971 std::vector<double> SpectralVarianceDown(nDraw, 0.0);
3972 std::vector<int> DenomCounterDown(nDraw, 0);
3973
3974 const int ThresholsCheck = Division*k*nSteps;
3975 //KS: First mean
3976 for (int j = 0; j < nDraw; ++j)
3977 {
3978 for(int i = 0; i < nEntries; ++i)
3979 {
3980 if(StepNumber[i] < ThresholsCheck)
3981 {
3982 MeanDown[j] += ParStep[j][i];
3983 DenomCounterDown[j]++;
3984 }
3985 }
3986 MeanDown[j] = MeanDown[j]/DenomCounterDown[j];
3987 }
3988 //Now spectral variance
3989 for (int j = 0; j < nDraw; ++j)
3990 {
3991 for(int i = 0; i < nEntries; ++i)
3992 {
3993 if(StepNumber[i] < ThresholsCheck)
3994 {
3995 SpectralVarianceDown[j] += (ParStep[j][i] - MeanDown[j])*(ParStep[j][i] - MeanDown[j]);
3996 }
3997 }
3998 }
3999 //Lastly calc T score and fill histogram entry
4000 for (int j = 0; j < nDraw; ++j)
4001 {
4002 double T_score = std::fabs((MeanDown[j] - MeanUp[j])/std::sqrt(SpectralVarianceDown[j]/DenomCounterDown[j] + SpectralVarianceUp[j]/DenomCounterUp[j]));
4003 GewekePlots[j]->SetBinContent(k, T_score);
4004 }
4005 } //end loop over intervals
4006#ifdef MULTITHREAD
4007} //End parallel region
4008#endif
4009
4010 //Finally save it to TFile
4011 OutputFile->cd();
4012 TDirectory *GewekeDir = OutputFile->mkdir("Geweke");
4013 for (int j = 0; j < nDraw; ++j)
4014 {
4015 GewekeDir->cd();
4016 GewekePlots[j]->Write();
4017 }
4018 for (int i = 0; i < nDraw; ++i) {
4019 delete[] ParStep[i];
4020 }
4021 delete[] ParStep;
4022
4023 GewekeDir->Close();
4024 delete GewekeDir;
4025 OutputFile->cd();
4026}

◆ Initialise()

void MCMCProcessor::Initialise ( )

Scan chain, what parameters we have and load information from covariance matrices.

Definition at line 149 of file MCMCProcessor.cpp.

149 {
150// ***************
151 // Scan the ROOT file for useful branches
152 ScanInput();
153
154 // Setup the output
155 SetupOutput();
156}
MCMCProcessor::SetupOutput
void SetupOutput()
Prepare all objects used for output.
Definition: MCMCProcessor.cpp:2025
MCMCProcessor::ScanInput
void ScanInput()
Scan Input etc.
Definition: MCMCProcessor.cpp:1920

◆ LoadAdditionalInfo()

virtual void MCMCProcessor::LoadAdditionalInfo ( )
inlineprotectedvirtual

allow loading additional info for example used for oscillation parameters

Reimplemented in OscProcessor.

Definition at line 293 of file MCMCProcessor.h.

293{};

◆ MakeCovariance()

void MCMCProcessor::MakeCovariance ( )

Calculate covariance by making 2D projection of each combination of parameters.

Definition at line 925 of file MCMCProcessor.cpp.

925 {
926// *********************
927 if (OutputFile == nullptr) MakeOutputFile();
928
929 bool HaveMadeDiagonal = false;
930 MACH3LOG_INFO("Making post-fit covariances...");
931 // Check that the diagonal entries have been filled
932 // i.e. MakePostfit() has been called
933 for (int i = 0; i < nDraw; ++i) {
934 if ((*Covariance)(i,i) == _UNDEF_) {
935 HaveMadeDiagonal = false;
936 MACH3LOG_INFO("Have not run diagonal elements in covariance, will do so now by calling MakePostfit()");
937 break;
938 } else {
939 HaveMadeDiagonal = true;
940 }
941 }
942
943 if (HaveMadeDiagonal == false) {
944 MakePostfit();
945 }
946 gStyle->SetPalette(55);
947 // Now we are sure we have the diagonal elements, let's make the off-diagonals
948 for (int i = 0; i < nDraw; ++i)
949 {
950 if (i % (nDraw/5) == 0)
951 MaCh3Utils::PrintProgressBar(i, nDraw);
952
953 TString Title_i = "";
954 double Prior_i, PriorError;
955
956 GetNthParameter(i, Prior_i, PriorError, Title_i);
957
958 const double min_i = Chain->GetMinimum(BranchNames[i]);
959 const double max_i = Chain->GetMaximum(BranchNames[i]);
960
961 // Loop over the other parameters to get the correlations
962 for (int j = 0; j <= i; ++j) {
963 // Skip the diagonal elements which we've already done above
964 if (j == i) continue;
965
966 // If this parameter isn't varied
967 if (IamVaried[j] == false) {
968 (*Covariance)(i,j) = 0.0;
969 (*Covariance)(j,i) = (*Covariance)(i,j);
970 (*Correlation)(i,j) = 0.0;
971 (*Correlation)(j,i) = (*Correlation)(i,j);
972 continue;
973 }
974
975 TString Title_j = "";
976 double Prior_j, PriorError_j;
977 GetNthParameter(j, Prior_j, PriorError_j, Title_j);
978
979 OutputFile->cd();
980
981 // The draw which we want to perform
982 TString DrawMe = BranchNames[j]+":"+BranchNames[i];
983
984 const double max_j = Chain->GetMaximum(BranchNames[j]);
985 const double min_j = Chain->GetMinimum(BranchNames[j]);
986
987 // TH2F to hold the Correlation
988 std::unique_ptr<TH2D> hpost_2D = std::make_unique<TH2D>(DrawMe, DrawMe, nBins, min_i, max_i, nBins, min_j, max_j);
989 hpost_2D->SetMinimum(0);
990 hpost_2D->GetXaxis()->SetTitle(Title_i);
991 hpost_2D->GetYaxis()->SetTitle(Title_j);
992 hpost_2D->GetZaxis()->SetTitle("Steps");
993
994 // The draw command we want, i.e. draw param j vs param i
995 Chain->Project(DrawMe, DrawMe, StepCut.c_str());
996
997 if(ApplySmoothing) hpost_2D->Smooth();
998 // Get the Covariance for these two parameters
999 (*Covariance)(i,j) = hpost_2D->GetCovariance();
1000 (*Covariance)(j,i) = (*Covariance)(i,j);
1001
1002 (*Correlation)(i,j) = hpost_2D->GetCorrelationFactor();
1003 (*Correlation)(j,i) = (*Correlation)(i,j);
1004
1005 if(printToPDF)
1006 {
1007 //KS: Skip Flux Params
1008 if(ParamType[i] == kXSecPar && ParamType[j] == kXSecPar)
1009 {
1010 if(std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold)
1011 {
1012 Posterior->cd();
1013 hpost_2D->Draw("colz");
1014 Posterior->SetName(hpost_2D->GetName());
1015 Posterior->SetTitle(hpost_2D->GetTitle());
1016 Posterior->Print(CanvasName);
1017 }
1018 }
1019 }
1020 // Write it to root file
1021 //OutputFile->cd();
1022 //if( std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold ) hpost_2D->Write();
1023 } // End j loop
1024 } // End i loop
1025 OutputFile->cd();
1026 Covariance->Write("Covariance");
1027 Correlation->Write("Correlation");
1028}

◆ MakeCovariance_MP()

void MCMCProcessor::MakeCovariance_MP ( const bool  Mute = false)

Calculate covariance by making 2D projection of each combination of parameters using multithreading.

Parameters
MuteAllow silencing many messages, especially important if we calculate matrix many times

Definition at line 1141 of file MCMCProcessor.cpp.

1141 {
1142// *********************
1143 if (OutputFile == nullptr) MakeOutputFile();
1144
1145 if(!CacheMCMC) CacheSteps();
1146
1147 bool HaveMadeDiagonal = false;
1148 // Check that the diagonal entries have been filled
1149 // i.e. MakePostfit() has been called
1150 for (int i = 0; i < nDraw; ++i) {
1151 if ((*Covariance)(i,i) == _UNDEF_) {
1152 HaveMadeDiagonal = false;
1153 MACH3LOG_WARN("Have not run diagonal elements in covariance, will do so now by calling MakePostfit()");
1154 break;
1155 } else {
1156 HaveMadeDiagonal = true;
1157 }
1158 }
1159
1160 if (HaveMadeDiagonal == false) MakePostfit();
1161 if(!Mute) MACH3LOG_INFO("Calculating covariance matrix");
1162 TStopwatch clock;
1163 if(!Mute) clock.Start();
1164
1165 gStyle->SetPalette(55);
1166 // Now we are sure we have the diagonal elements, let's make the off-diagonals
1167 #ifdef MULTITHREAD
1168 #pragma omp parallel for
1169 #endif
1170 for (int i = 0; i < nDraw; ++i)
1171 {
1172 for (int j = 0; j <= i; ++j)
1173 {
1174 // Skip the diagonal elements which we've already done above
1175 if (j == i) continue;
1176
1177 // If this parameter isn't varied
1178 if (IamVaried[j] == false) {
1179 (*Covariance)(i,j) = 0.0;
1180 (*Covariance)(j,i) = (*Covariance)(i,j);
1181 (*Correlation)(i,j) = 0.0;
1182 (*Correlation)(j,i) = (*Correlation)(i,j);
1183 continue;
1184 }
1185 hpost2D[i][j]->SetMinimum(0);
1186
1187 for (int k = 0; k < nEntries; ++k)
1188 {
1189 //KS: Burn in cut
1190 if(StepNumber[k] < BurnInCut) continue;
1191
1192 //KS: Fill histogram with cached steps
1193 hpost2D[i][j]->Fill(ParStep[i][k], ParStep[j][k]);
1194 }
1195 if(ApplySmoothing) hpost2D[i][j]->Smooth();
1196
1197 // Get the Covariance for these two parameters
1198 (*Covariance)(i,j) = hpost2D[i][j]->GetCovariance();
1199 (*Covariance)(j,i) = (*Covariance)(i,j);
1200
1201 //KS: Since we already have covariance consider calculating correlation using it, right now we effectively calculate covariance twice
1202 //https://root.cern.ch/doc/master/TH2_8cxx_source.html#l01099
1203 (*Correlation)(i,j) = hpost2D[i][j]->GetCorrelationFactor();
1204 (*Correlation)(j,i) = (*Correlation)(i,j);
1205 }// End j loop
1206 }// End i loop
1207
1208 if(!Mute) {
1209 clock.Stop();
1210 MACH3LOG_INFO("Making Covariance took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries);
1211 }
1212 OutputFile->cd();
1213 if(printToPDF)
1214 {
1215 Posterior->cd();
1216 for (int i = 0; i < nDraw; ++i)
1217 {
1218 for (int j = 0; j <= i; ++j)
1219 {
1220 // Skip the diagonal elements which we've already done above
1221 if (j == i) continue;
1222 if (IamVaried[j] == false) continue;
1223
1224 if(ParamType[i] == kXSecPar && ParamType[j] == kXSecPar)
1225 {
1226 if(std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold)
1227 {
1228 hpost2D[i][j]->Draw("colz");
1229 Posterior->SetName(hpost2D[i][j]->GetName());
1230 Posterior->SetTitle(hpost2D[i][j]->GetTitle());
1231 Posterior->Print(CanvasName);
1232 }
1233 }
1234 //if( std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold) hpost2D[i][j]->Write();
1235 }// End j loop
1236 }// End i loop
1237 } //end if pdf
1238 if(!Mute) {
1239 Covariance->Write("Covariance");
1240 Correlation->Write("Correlation");
1241 }
1242}
MCMCProcessor::BurnInCut
int BurnInCut
Value of burn in cut.
Definition: MCMCProcessor.h:379
MCMCProcessor::GetCovariance
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
Definition: MCMCProcessor.cpp:192
MCMCProcessor::CacheSteps
void CacheSteps()
KS:By caching each step we use multithreading.
Definition: MCMCProcessor.cpp:1032

◆ MakeCredibleIntervals()

void MCMCProcessor::MakeCredibleIntervals ( const std::vector< double > &  CredibleIntervals = {0.99, 0.90, 0.68 },
const std::vector< Color_t > &  CredibleIntervalsColours = {kCyan+4, kCyan-2, kCyan-10},
const bool  CredibleInSigmas = false 
)

Make and Draw Credible intervals.

Parameters
CredibleIntervalsVector with values of credible intervals, must be in descending order
CredibleIntervalsColoursColor_t telling what colour to use for each Interval line
CredibleInSigmasBool telling whether intervals are in percentage or in sigmas, then special conversions is used

Scale the histograms so it shows the posterior probability

Definition at line 674 of file MCMCProcessor.cpp.

676 {
677// *********************
678 if(hpost[0] == nullptr) MakePostfit();
679
680 MACH3LOG_INFO("Making Credible Intervals ");
681 const double LeftMargin = Posterior->GetLeftMargin();
682 Posterior->SetLeftMargin(0.15);
683
684 // KS: Sanity check of size and ordering is correct
685 CheckCredibleIntervalsOrder(CredibleIntervals, CredibleIntervalsColours);
686 const int nCredible = int(CredibleIntervals.size());
687 std::vector<std::unique_ptr<TH1D>> hpost_copy(nDraw);
688 std::vector<std::vector<std::unique_ptr<TH1D>>> hpost_cl(nDraw);
689
690 //KS: Copy all histograms to be thread safe
691 for (int i = 0; i < nDraw; ++i)
692 {
693 hpost_copy[i] = M3::Clone<TH1D>(hpost[i], Form("hpost_copy_%i", i));
694 hpost_cl[i].resize(nCredible);
695 for (int j = 0; j < nCredible; ++j)
696 {
697 hpost_cl[i][j] = M3::Clone<TH1D>(hpost[i], Form("hpost_copy_%i_CL_%f", i, CredibleIntervals[j]));
698
699 //KS: Reset to get rid to TF1 otherwise we run into segfault :(
700 hpost_cl[i][j]->Reset("");
701 hpost_cl[i][j]->Fill(0.0, 0.0);
702 }
703 }
704
705 #ifdef MULTITHREAD
706 #pragma omp parallel for
707 #endif
708 for (int i = 0; i < nDraw; ++i)
709 {
711 hpost_copy[i]->Scale(1. / hpost_copy[i]->Integral());
712 for (int j = 0; j < nCredible; ++j)
713 {
714 // Scale the histograms before gettindg credible intervals
715 hpost_cl[i][j]->Scale(1. / hpost_cl[i][j]->Integral());
716 GetCredibleIntervalSig(hpost_copy[i], hpost_cl[i][j], CredibleInSigmas, CredibleIntervals[j]);
717
718 hpost_cl[i][j]->SetFillColor(CredibleIntervalsColours[j]);
719 hpost_cl[i][j]->SetLineWidth(1);
720 }
721 hpost_copy[i]->GetYaxis()->SetTitleOffset(1.8);
722 hpost_copy[i]->SetLineWidth(1);
723 hpost_copy[i]->SetMaximum(hpost_copy[i]->GetMaximum()*1.2);
724 hpost_copy[i]->SetLineWidth(2);
725 hpost_copy[i]->SetLineColor(kBlack);
726 hpost_copy[i]->GetYaxis()->SetTitle("Posterior Probability");
727 }
728
729 OutputFile->cd();
730 TDirectory *CredibleDir = OutputFile->mkdir("Credible");
731
732 for (int i = 0; i < nDraw; ++i)
733 {
734 if(!IamVaried[i]) continue;
735
736 // Now make the TLine for the Asimov
737 TString Title = "";
738 double Prior = 1.0, PriorError = 1.0;
739 GetNthParameter(i, Prior, PriorError, Title);
740
741 auto Asimov = std::make_unique<TLine>(Prior, hpost_copy[i]->GetMinimum(), Prior, hpost_copy[i]->GetMaximum());
742 SetTLineStyle(Asimov.get(), kRed-3, 2, kDashed);
743
744 auto legend = std::make_unique<TLegend>(0.20, 0.7, 0.4, 0.92);
745 SetLegendStyle(legend.get(), 0.03);
746 hpost_copy[i]->Draw("HIST");
747
748 for (int j = 0; j < nCredible; ++j)
749 hpost_cl[i][j]->Draw("HIST SAME");
750 for (int j = nCredible-1; j >= 0; --j)
751 {
752 if(CredibleInSigmas)
753 legend->AddEntry(hpost_cl[i][j].get(), Form("%.0f#sigma Credible Interval", CredibleIntervals[j]), "f");
754 else
755 legend->AddEntry(hpost_cl[i][j].get(), Form("%.0f%% Credible Interval", CredibleIntervals[j]*100), "f");
756 }
757 legend->AddEntry(Asimov.get(), Form("#splitline{Prior}{x = %.2f , #sigma = %.2f}", Prior, PriorError), "l");
758 legend->Draw("SAME");
759 Asimov->Draw("SAME");
760
761 // Write to file
762 Posterior->SetName(hpost[i]->GetName());
763 Posterior->SetTitle(hpost[i]->GetTitle());
764
765 if(printToPDF) Posterior->Print(CanvasName);
766 // cd into directory in root file
767 CredibleDir->cd();
768 Posterior->Write();
769 }
770 CredibleDir->Close();
771 delete CredibleDir;
772
773 OutputFile->cd();
774
775 //Set back to normal
776 Posterior->SetLeftMargin(LeftMargin);
777}
GetCredibleIntervalSig
void GetCredibleIntervalSig(const std::unique_ptr< TH1D > &hist, std::unique_ptr< TH1D > &hpost_copy, const bool CredibleInSigmas, const double coverage)
KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning,...
Definition: StatisticalUtils.cpp:367
MCMCProcessor::SetTLineStyle
void SetTLineStyle(TLine *Line, const Color_t Colour, const Width_t Width, const ELineStyle Style) const
Configures a TLine object with the specified style parameters.
Definition: MCMCProcessor.cpp:4182
MCMCProcessor::CheckCredibleIntervalsOrder
void CheckCredibleIntervalsOrder(const std::vector< double > &CredibleIntervals, const std::vector< Color_t > &CredibleIntervalsColours)
Checks the order and size consistency of the CredibleIntervals and CredibleIntervalsColours vectors.
Definition: MCMCProcessor.cpp:4076

◆ MakeCredibleRegions()

void MCMCProcessor::MakeCredibleRegions ( const std::vector< double > &  CredibleRegions = {0.99, 0.90, 0.68},
const std::vector< Style_t > &  CredibleRegionStyle = {kDashed, kSolid, kDotted},
const std::vector< Color_t > &  CredibleRegionColor = {kGreen-3, kGreen-10, kGreen},
const bool  CredibleInSigmas = false,
const bool  Draw2DPosterior = true,
const bool  DrawBestFit = true 
)

Make and Draw Credible Regions.

Parameters
CredibleRegionsVector with values of credible intervals, must be in descending order
CredibleRegionStyleStyle_t telling what line style to use for each Interval line
CredibleRegionColorColor_t telling what colour to use for each Interval line
CredibleInSigmasBool telling whether intervals are in percentage or in sigmas, then special conversions is used

Definition at line 1544 of file MCMCProcessor.cpp.

1549 {
1550// *********************
1551 if(hpost2D.size() == 0) MakeCovariance_MP();
1552 MACH3LOG_INFO("Making Credible Regions");
1553
1554 CheckCredibleRegionsOrder(CredibleRegions, CredibleRegionStyle, CredibleRegionColor);
1555 const int nCredible = int(CredibleRegions.size());
1556
1557 std::vector<std::vector<std::unique_ptr<TH2D>>> hpost_2D_copy(nDraw);
1558 std::vector<std::vector<std::vector<std::unique_ptr<TH2D>>>> hpost_2D_cl(nDraw);
1559 //KS: Copy all histograms to be thread safe
1560 for (int i = 0; i < nDraw; ++i)
1561 {
1562 hpost_2D_copy[i].resize(nDraw);
1563 hpost_2D_cl[i].resize(nDraw);
1564 for (int j = 0; j <= i; ++j)
1565 {
1566 hpost_2D_copy[i][j] = M3::Clone<TH2D>(hpost2D[i][j], Form("hpost_copy_%i_%i", i, j));
1567 hpost_2D_cl[i][j].resize(nCredible);
1568 for (int k = 0; k < nCredible; ++k)
1569 {
1570 hpost_2D_cl[i][j][k] = M3::Clone<TH2D>(hpost2D[i][j], Form("hpost_copy_%i_%i_CL_%f", i, j, CredibleRegions[k]));
1571 }
1572 }
1573 }
1574
1575 #ifdef MULTITHREAD
1576 #pragma omp parallel for
1577 #endif
1578 //Calculate credible histogram
1579 for (int i = 0; i < nDraw; ++i)
1580 {
1581 for (int j = 0; j <= i; ++j)
1582 {
1583 for (int k = 0; k < nCredible; ++k)
1584 {
1585 GetCredibleRegionSig(hpost_2D_cl[i][j][k], CredibleInSigmas, CredibleRegions[k]);
1586 hpost_2D_cl[i][j][k]->SetLineColor(CredibleRegionColor[k]);
1587 hpost_2D_cl[i][j][k]->SetLineWidth(2);
1588 hpost_2D_cl[i][j][k]->SetLineStyle(CredibleRegionStyle[k]);
1589 }
1590 }
1591 }
1592
1593 gStyle->SetPalette(51);
1594 for (int i = 0; i < nDraw; ++i)
1595 {
1596 for (int j = 0; j <= i; ++j)
1597 {
1598 // Skip the diagonal elements which we've already done above
1599 if (j == i) continue;
1600 if (IamVaried[j] == false) continue;
1601
1602 auto legend = std::make_unique<TLegend>(0.20, 0.7, 0.4, 0.92);
1603 legend->SetTextColor(kRed);
1604 SetLegendStyle(legend.get(), 0.03);
1605
1606 //Get Best point
1607 auto bestfitM = std::make_unique<TGraph>(1);
1608 const int MaxBin = hpost_2D_copy[i][j]->GetMaximumBin();
1609 int Mbx, Mby, Mbz;
1610 hpost_2D_copy[i][j]->GetBinXYZ(MaxBin, Mbx, Mby, Mbz);
1611 const double Mx = hpost_2D_copy[i][j]->GetXaxis()->GetBinCenter(Mbx);
1612 const double My = hpost_2D_copy[i][j]->GetYaxis()->GetBinCenter(Mby);
1613
1614 bestfitM->SetPoint(0, Mx, My);
1615 bestfitM->SetMarkerStyle(22);
1616 bestfitM->SetMarkerSize(1);
1617 bestfitM->SetMarkerColor(kMagenta);
1618
1619 //Plot default 2D posterior
1620
1621 if(Draw2DPosterior){
1622 hpost_2D_copy[i][j]->Draw("COLZ");
1623 }
1624 else{
1625 hpost_2D_copy[i][j]->Draw("AXIS");
1626 }
1627
1628 //Now credible regions
1629 for (int k = 0; k < nCredible; ++k)
1630 hpost_2D_cl[i][j][k]->Draw("CONT3 SAME");
1631 for (int k = nCredible-1; k >= 0; --k)
1632 {
1633 if(CredibleInSigmas)
1634 legend->AddEntry(hpost_2D_cl[i][j][k].get(), Form("%.0f#sigma Credible Interval", CredibleRegions[k]), "l");
1635 else
1636 legend->AddEntry(hpost_2D_cl[i][j][k].get(), Form("%.0f%% Credible Region", CredibleRegions[k]*100), "l");
1637 }
1638 legend->Draw("SAME");
1639
1640 if(DrawBestFit){
1641 legend->AddEntry(bestfitM.get(),"Best Fit","p");
1642 bestfitM->Draw("SAME.P");
1643 }
1644
1645 // Write to file
1646 Posterior->SetName(hpost2D[i][j]->GetName());
1647 Posterior->SetTitle(hpost2D[i][j]->GetTitle());
1648
1649 //KS: Print only regions with correlation greater than specified value, by default 0.2. This is done to avoid dumping thousands of plots
1650 if(printToPDF && std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold) Posterior->Print(CanvasName);
1651 // Write it to root file
1652 //OutputFile->cd();
1653 //if( std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold ) Posterior->Write();
1654 }
1655 }
1656
1657 OutputFile->cd();
1658}
GetCredibleRegionSig
void GetCredibleRegionSig(std::unique_ptr< TH2D > &hist2D, const bool CredibleInSigmas, const double coverage)
KS: Set 2D contour within some coverage.
Definition: StatisticalUtils.cpp:431
MCMCProcessor::CheckCredibleRegionsOrder
void CheckCredibleRegionsOrder(const std::vector< double > &CredibleRegions, const std::vector< Style_t > &CredibleRegionStyle, const std::vector< Color_t > &CredibleRegionColor)
Checks the order and size consistency of the CredibleRegions, CredibleRegionStyle,...
Definition: MCMCProcessor.cpp:4095

◆ MakeOutputFile()

void MCMCProcessor::MakeOutputFile ( )
inlineprotected

prepare output root file and canvas to which we will save EVERYTHING

Definition at line 201 of file MCMCProcessor.cpp.

201 {
202// ***************
203 //KS: ROOT hates me... but we can create several instances of MCMC Processor, each with own TCanvas ROOT is mad and will delete if there is more than one canvas with the same name, so we add random number to avoid issue
204 auto rand = std::make_unique<TRandom3>(0);
205 const int uniform = int(rand->Uniform(0, 10000));
206 // Open a TCanvas to write the posterior onto
207 Posterior = std::make_unique<TCanvas>(("Posterior" + std::to_string(uniform)).c_str(), ("Posterior" + std::to_string(uniform)).c_str(), 0, 0, 1024, 1024);
208 //KS: No idea why but ROOT changed treatment of violin in R6. If you have non uniform binning this will results in very hard to see violin plots.
209 TCandle::SetScaledViolin(false);
210
211 Posterior->SetGrid();
212 gStyle->SetOptStat(0);
213 gStyle->SetOptTitle(0);
214 Posterior->SetTickx();
215 Posterior->SetTicky();
216
217 Posterior->SetBottomMargin(0.1);
218 Posterior->SetTopMargin(0.05);
219 Posterior->SetRightMargin(0.03);
220 Posterior->SetLeftMargin(0.15);
221
222 //To avoid TCanvas::Print> messages
223 gErrorIgnoreLevel = kWarning;
224
225 // Output file to write to
226 OutputName = MCMCFile + OutputSuffix +".root";
227
228 // Output file
229 OutputFile = new TFile(OutputName.c_str(), "recreate");
230 OutputFile->cd();
231}
MCMCProcessor::OutputName
std::string OutputName
Name of output files.
Definition: MCMCProcessor.h:425

◆ MakePostfit()

void MCMCProcessor::MakePostfit ( )

Make 1D projection for each parameter and prepare structure.

Definition at line 235 of file MCMCProcessor.cpp.

235 {
236// ****************************
237 // Check if we've already made post-fit
238 if (MadePostfit == true) return;
239 MadePostfit = true;
240
241 // Check if the output file is ready
242 if (OutputFile == nullptr) MakeOutputFile();
243
244 MACH3LOG_INFO("MCMCProcessor is making post-fit plots...");
245
246 int originalErrorLevel = gErrorIgnoreLevel;
247 gErrorIgnoreLevel = kFatal;
248
249 // Directory for posteriors
250 TDirectory *PostDir = OutputFile->mkdir("Post");
251 TDirectory *PostHistDir = OutputFile->mkdir("Post_1d_hists");
252
253 // nDraw is number of draws we want to do
254 for (int i = 0; i < nDraw; ++i)
255 {
256 if (i % (nDraw/5) == 0) {
257 MaCh3Utils::PrintProgressBar(i, nDraw);
258 }
259 OutputFile->cd();
260 TString Title = "";
261 double Prior = 1.0, PriorError = 1.0;
262 GetNthParameter(i, Prior, PriorError, Title);
263
264 // This holds the posterior density
265 const double maxi = Chain->GetMaximum(BranchNames[i]);
266 const double mini = Chain->GetMinimum(BranchNames[i]);
267 // This holds the posterior density
268 hpost[i] = new TH1D(BranchNames[i], BranchNames[i], nBins, mini, maxi);
269 hpost[i]->SetMinimum(0);
270 hpost[i]->GetYaxis()->SetTitle("Steps");
271 hpost[i]->GetYaxis()->SetNoExponent(false);
272
273 //KS: Apply additional Cuts, like mass ordering
274 std::string CutPosterior1D = "";
275 if(Posterior1DCut != "")
276 {
277 CutPosterior1D = StepCut +" && " + Posterior1DCut;
278 }
279 else CutPosterior1D = StepCut;
280
281 // Project BranchNames[i] onto hpost, applying stepcut
282 Chain->Project(BranchNames[i], BranchNames[i], CutPosterior1D.c_str());
283
284 if(ApplySmoothing) hpost[i]->Smooth();
285
286 (*Central_Value)(i) = Prior;
287
288 double Mean, Err, Err_p, Err_m;
289 GetArithmetic(hpost[i], Mean, Err);
290 (*Means)(i) = Mean;
291 (*Errors)(i) = Err;
292
293 GetGaussian(hpost[i], Gauss.get(), Mean, Err);
294 (*Means_Gauss)(i) = Mean;
295 (*Errors_Gauss)(i) = Err;
296
297 GetHPD(hpost[i], Mean, Err, Err_p, Err_m);
298 (*Means_HPD)(i) = Mean;
299 (*Errors_HPD)(i) = Err;
300 (*Errors_HPD_Positive)(i) = Err_p;
301 (*Errors_HPD_Negative)(i) = Err_m;
302
303 // Write the results from the projection into the TVectors and TMatrices
304 (*Covariance)(i,i) = (*Errors)(i)*(*Errors)(i);
305 (*Correlation)(i,i) = 1.0;
306
307 //KS: This need to be before SetMaximum(), this way plot is nicer as line end at the maximum
308 auto hpd = std::make_unique<TLine>((*Means_HPD)(i), hpost[i]->GetMinimum(), (*Means_HPD)(i), hpost[i]->GetMaximum());
309 SetTLineStyle(hpd.get(), kBlack, 2, kSolid);
310
311 hpost[i]->SetLineWidth(2);
312 hpost[i]->SetLineColor(kBlue-1);
313 hpost[i]->SetMaximum(hpost[i]->GetMaximum()*DrawRange);
314 hpost[i]->SetTitle(Title);
315 hpost[i]->GetXaxis()->SetTitle(hpost[i]->GetTitle());
316
317 // Now make the TLine for the Asimov
318 auto Asimov = std::make_unique<TLine>(Prior, hpost[i]->GetMinimum(), Prior, hpost[i]->GetMaximum());
319 SetTLineStyle(Asimov.get(), kRed-3, 2, kDashed);
320
321 auto leg = std::make_unique<TLegend>(0.12, 0.6, 0.6, 0.97);
322 SetLegendStyle(leg.get(), 0.04);
323 leg->AddEntry(hpost[i], Form("#splitline{PDF}{#mu = %.2f, #sigma = %.2f}", hpost[i]->GetMean(), hpost[i]->GetRMS()), "l");
324 leg->AddEntry(Gauss.get(), Form("#splitline{Gauss}{#mu = %.2f, #sigma = %.2f}", Gauss->GetParameter(1), Gauss->GetParameter(2)), "l");
325 leg->AddEntry(hpd.get(), Form("#splitline{HPD}{#mu = %.2f, #sigma = %.2f (+%.2f-%.2f)}", (*Means_HPD)(i), (*Errors_HPD)(i), (*Errors_HPD_Positive)(i), (*Errors_HPD_Negative)(i)), "l");
326 leg->AddEntry(Asimov.get(), Form("#splitline{Prior}{x = %.2f , #sigma = %.2f}", Prior, PriorError), "l");
327
328 //CW: Don't plot if this is a fixed histogram (i.e. the peak is the whole integral)
329 if (hpost[i]->GetMaximum() == hpost[i]->Integral()*DrawRange)
330 {
331 MACH3LOG_WARN("Found fixed parameter: {} ({}), moving on", Title, i);
332 IamVaried[i] = false;
333 //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.
334 (*Means_HPD)(i) = Prior;
335 (*Errors_HPD)(i) = PriorError;
336 continue;
337 }
338
339 // Store that this parameter is indeed being varied
340 IamVaried[i] = true;
341
342 // Write to file
343 Posterior->SetName(Title);
344 Posterior->SetTitle(Title);
345
346 // Draw onto the TCanvas
347 hpost[i]->Draw();
348 hpd->Draw("same");
349 Asimov->Draw("same");
350 leg->Draw("same");
351
352 if(printToPDF) Posterior->Print(CanvasName);
353
354 // cd into params directory in root file
355 PostDir->cd();
356 Posterior->Write();
357
358 hpost[i]->SetName(Title);
359 hpost[i]->SetTitle(Title);
360 PostHistDir->cd();
361 hpost[i]->Write();
362 } // end the for loop over nDraw
363
364 OutputFile->cd();
365 TTree *SettingsBranch = new TTree("Settings", "Settings");
366 int CrossSectionParameters = nParam[kXSecPar];
367 SettingsBranch->Branch("CrossSectionParameters", &CrossSectionParameters);
368 int CrossSectionParametersStartingPos = ParamTypeStartPos[kXSecPar];
369 SettingsBranch->Branch("CrossSectionParametersStartingPos", &CrossSectionParametersStartingPos);
370 int FluxParameters = GetGroup("Flux");
371 SettingsBranch->Branch("FluxParameters", &FluxParameters);
372
373 int NDParameters = nParam[kNDPar];
374 SettingsBranch->Branch("NDParameters", &NDParameters);
375 int NDParametersStartingPos = ParamTypeStartPos[kNDPar];
376 SettingsBranch->Branch("NDParametersStartingPos", &NDParametersStartingPos);
377
378 int FDParameters = nParam[kFDDetPar];
379 SettingsBranch->Branch("FDParameters", &FDParameters);
380 int FDParametersStartingPos = ParamTypeStartPos[kFDDetPar];
381 SettingsBranch->Branch("FDParametersStartingPos", &FDParametersStartingPos);
382
383 SettingsBranch->Branch("NDSamplesBins", &NDSamplesBins);
384 SettingsBranch->Branch("NDSamplesNames", &NDSamplesNames);
385
386 SettingsBranch->Fill();
387 SettingsBranch->Write();
388 delete SettingsBranch;
389
390 TDirectory *Names = OutputFile->mkdir("Names");
391 Names->cd();
392 for (std::vector<TString>::iterator it = BranchNames.begin(); it != BranchNames.end(); ++it) {
393 TObjString((*it)).Write();
394 }
395 Names->Close();
396 delete Names;
397
398 OutputFile->cd();
399 Central_Value->Write("Central_Value");
400 Means->Write("PDF_Means");
401 Errors->Write("PDF_Error");
402 Means_Gauss->Write("Gauss_Means");
403 Errors_Gauss->Write("Gauss_Errors");
404 Means_HPD->Write("Means_HPD");
405 Errors_HPD->Write("Errors_HPD");
406 Errors_HPD_Positive->Write("Errors_HPD_Positive");
407 Errors_HPD_Negative->Write("Errors_HPD_Negative");
408
409 PostDir->Close();
410 delete PostDir;
411 PostHistDir->Close();
412 delete PostHistDir;
413
414 // restore original warning setting
415 gErrorIgnoreLevel = originalErrorLevel;
416} // Have now written the postfit projections
FDParameters
int FDParameters
Definition: GetPostfitParamPlots.cpp:59
CrossSectionParameters
int CrossSectionParameters
Definition: GetPostfitParamPlots.cpp:62
FDParametersStartingPos
int FDParametersStartingPos
Definition: GetPostfitParamPlots.cpp:60
NDParametersStartingPos
int NDParametersStartingPos
Definition: GetPostfitParamPlots.cpp:57
NDParameters
int NDParameters
Definition: GetPostfitParamPlots.cpp:56
FluxParameters
int FluxParameters
Definition: GetPostfitParamPlots.cpp:65
Mean
double ** Mean
Definition: RHat.cpp:58
GetGaussian
void GetGaussian(TH1D *&hist, TF1 *gauss, double &Mean, double &Error)
CW: Fit Gaussian to posterior.
Definition: StatisticalUtils.cpp:208
GetHPD
void GetHPD(TH1D *const hist, double &Mean, double &Error, double &Error_p, double &Error_m, const double coverage)
Get Highest Posterior Density (HPD)
Definition: StatisticalUtils.cpp:228
GetArithmetic
void GetArithmetic(TH1D *const hist, double &Mean, double &Error)
CW: Get Arithmetic mean from posterior.
Definition: StatisticalUtils.cpp:201
MCMCProcessor::Gauss
std::unique_ptr< TF1 > Gauss
Gaussian fitter.
Definition: MCMCProcessor.h:454

◆ MakePrefit()

std::unique_ptr< TH1D > MCMCProcessor::MakePrefit ( )
inlineprotected

Prepare prefit histogram for parameter overlay plot.

Definition at line 2100 of file MCMCProcessor.cpp.

2100 {
2101// *****************************
2102 if (OutputFile == nullptr) MakeOutputFile();
2103
2104 auto PreFitPlot = std::make_unique<TH1D>("Prefit", "Prefit", nDraw, 0, nDraw);
2105 PreFitPlot->SetDirectory(nullptr);
2106 for (int i = 0; i < PreFitPlot->GetNbinsX() + 1; ++i) {
2107 PreFitPlot->SetBinContent(i+1, 0);
2108 PreFitPlot->SetBinError(i+1, 0);
2109 }
2110
2111 //KS: Slightly hacky way to get relative to prior or nominal as this is convention we use,
2112 //Only applies for xsec, for other systematic it make no difference
2113 double CentralValueTemp, Central, Error;
2114
2115 // Set labels and data
2116 for (int i = 0; i < nDraw; ++i)
2117 {
2118 //Those keep which parameter type we run currently and relative number
2119 int ParamEnum = ParamType[i];
2120 int ParamNo = i - ParamTypeStartPos[ParameterEnum(ParamEnum)];
2121 CentralValueTemp = ParamCentral[ParamEnum][ParamNo];
2122 if(plotRelativeToPrior)
2123 {
2124 // Normalise the prior relative the nominal/prior, just the way we get our fit results in MaCh3
2125 if ( CentralValueTemp != 0) {
2126 Central = ParamCentral[ParamEnum][ParamNo] / CentralValueTemp;
2127 Error = ParamErrors[ParamEnum][ParamNo]/CentralValueTemp;
2128 } else {
2129 Central = CentralValueTemp + 1.0;
2130 Error = ParamErrors[ParamEnum][ParamNo];
2131 }
2132 }
2133 else
2134 {
2135 Central = CentralValueTemp;
2136 Error = ParamErrors[ParamEnum][ParamNo];
2137 }
2138 //KS: If plotting error for param with flat prior is turned off and given param really has flat prior set error to 0
2139 if(!PlotFlatPrior && ParamFlat[ParamEnum][ParamNo]) {
2140 Error = 0.;
2141 }
2142 PreFitPlot->SetBinContent(i+1, Central);
2143 PreFitPlot->SetBinError(i+1, Error);
2144 PreFitPlot->GetXaxis()->SetBinLabel(i+1, ParamNames[ParamEnum][ParamNo]);
2145 }
2146 PreFitPlot->SetDirectory(nullptr);
2147
2148 PreFitPlot->SetFillStyle(1001);
2149 PreFitPlot->SetFillColor(kRed-3);
2150 PreFitPlot->SetMarkerStyle(21);
2151 PreFitPlot->SetMarkerSize(2.4);
2152 PreFitPlot->SetMarkerColor(kWhite);
2153 PreFitPlot->SetLineColor(PreFitPlot->GetFillColor());
2154 PreFitPlot->GetXaxis()->LabelsOption("v");
2155
2156 return PreFitPlot;
2157}

◆ MakeSubOptimality()

void MCMCProcessor::MakeSubOptimality ( const int  NIntervals = 10)

Make and Draw SubOptimality [24].

Author
Henry Wallace

Definition at line 1247 of file MCMCProcessor.cpp.

1247 {
1248// *********************
1249
1250 //Save burn in cut, at the end of the loop we will return to default values
1251 const int DefaultUpperCut = UpperCut;
1252 const int DefaultBurnInCut = BurnInCut;
1253 bool defaultPrintToPDF = printToPDF;
1254 BurnInCut = 0;
1255 UpperCut = 0;
1256 printToPDF = false;
1257
1258 //Set via config in future
1259 int MaxStep = nSteps;
1260 int MinStep = 0;
1261 const int IntervalsSize = nSteps/NIntervals;
1262
1263 MACH3LOG_INFO("Making Suboptimality");
1264 TStopwatch clock;
1265 clock.Start();
1266
1267 std::unique_ptr<TH1D> SubOptimality = std::make_unique<TH1D>("Suboptimality", "Suboptimality", NIntervals, MinStep, MaxStep);
1268 SubOptimality->GetXaxis()->SetTitle("Step");
1269 SubOptimality->GetYaxis()->SetTitle("Suboptimality");
1270 SubOptimality->SetLineWidth(2);
1271 SubOptimality->SetLineColor(kBlue);
1272
1273 for(int i = 0; i < NIntervals; ++i)
1274 {
1275 //Reset our cov matrix
1276 ResetHistograms();
1277
1278 //Set threshold for calculating new matrix
1279 UpperCut = i*IntervalsSize;
1280 //Calculate cov matrix
1281 MakeCovariance_MP(true);
1282
1283 //Calculate eigen values
1284 TMatrixDSymEigen eigen(*Covariance);
1285 TVectorD eigen_values;
1286 eigen_values.ResizeTo(eigen.GetEigenValues());
1287 eigen_values = eigen.GetEigenValues();
1288
1289 //KS: Converting from ROOT to vector as to make using other libraires (Eigen) easier in future
1290 std::vector<double> EigenValues(eigen_values.GetNrows());
1291 for(unsigned int j = 0; j < EigenValues.size(); j++)
1292 {
1293 EigenValues[j] = eigen_values(j);
1294 }
1295 const double SubOptimalityValue = GetSubOptimality(EigenValues, nDraw);
1296 SubOptimality->SetBinContent(i+1, SubOptimalityValue);
1297 }
1298 clock.Stop();
1299 MACH3LOG_INFO("Making Suboptimality took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries);
1300
1301 UpperCut = DefaultUpperCut;
1302 BurnInCut = DefaultBurnInCut;
1303 printToPDF = defaultPrintToPDF;
1304
1305 SubOptimality->Draw("l");
1306 Posterior->SetName(SubOptimality->GetName());
1307 Posterior->SetTitle(SubOptimality->GetTitle());
1308
1309 if(printToPDF) Posterior->Print(CanvasName);
1310 // Write it to root file
1311 OutputFile->cd();
1312 Posterior->Write();
1313}
GetSubOptimality
double GetSubOptimality(const std::vector< double > &EigenValues, const int TotalTarameters)
Based on .
Definition: StatisticalUtils.cpp:184
MCMCProcessor::ResetHistograms
void ResetHistograms()
Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.
Definition: MCMCProcessor.cpp:2458

◆ MakeTrianglePlot()

void MCMCProcessor::MakeTrianglePlot ( const std::vector< std::string > &  ParNames,
const std::vector< double > &  CredibleIntervals = {0.99, 0.90, 0.68 },
const std::vector< Color_t > &  CredibleIntervalsColours = {kCyan+4, kCyan-2, kCyan-10},
const std::vector< double > &  CredibleRegions = {0.99, 0.90, 0.68},
const std::vector< Style_t > &  CredibleRegionStyle = {kDashed, kSolid, kDotted},
const std::vector< Color_t > &  CredibleRegionColor = {kGreen-3, kGreen-10, kGreen},
const bool  CredibleInSigmas = false 
)

Make fancy triangle plot for selected parameters.

Parameters
CredibleIntervalsVector with values of credible intervals, must be in descending order
CredibleIntervalsColoursColor_t telling what colour to use for each Interval line
CredibleInSigmasBool telling whether intervals are in percentage or in sigmas, then special conversions is used
CredibleRegionsVector with values of credible intervals, must be in descending order
CredibleRegionStyleStyle_t telling what line style to use for each Interval line
CredibleRegionColorColor_t telling what colour to use for each Interval line

Scale the histograms so it shows the posterior probability

Definition at line 1662 of file MCMCProcessor.cpp.

1671 {
1672// *********************
1673 if(hpost2D.size() == 0) MakeCovariance_MP();
1674 MACH3LOG_INFO("Making Triangle Plot");
1675
1676 const int nParamPlot = int(ParNames.size());
1677 std::vector<int> ParamNumber;
1678 for(int j = 0; j < nParamPlot; ++j)
1679 {
1680 int ParamNo = GetParamIndexFromName(ParNames[j]);
1681 if(ParamNo == _UNDEF_)
1682 {
1683 MACH3LOG_WARN("Couldn't find param {}. Will not plot Triangle plot", ParNames[j]);
1684 return;
1685 }
1686 ParamNumber.push_back(ParamNo);
1687 }
1688
1689 //KS: Store it as we go back to them at the end
1690 const std::vector<double> Margins = GetMargins(Posterior);
1691 Posterior->SetTopMargin(0.001);
1692 Posterior->SetBottomMargin(0.001);
1693 Posterior->SetLeftMargin(0.001);
1694 Posterior->SetRightMargin(0.001);
1695
1696 // KS: We later format hist several times so make one unfired lambda
1697 auto FormatHistogram = [](auto& hist) {
1698 hist->GetXaxis()->SetTitle("");
1699 hist->GetYaxis()->SetTitle("");
1700 hist->SetTitle("");
1701
1702 hist->GetXaxis()->SetLabelSize(0.1);
1703 hist->GetYaxis()->SetLabelSize(0.1);
1704
1705 hist->GetXaxis()->SetNdivisions(4);
1706 hist->GetYaxis()->SetNdivisions(4);
1707 };
1708
1709 Posterior->cd();
1710 Posterior->Clear();
1711 Posterior->Update();
1712
1713 //KS: We sort to have parameters from highest to lowest, this is related to how we make 2D projections in MakeCovariance_MP
1714 std::sort(ParamNumber.begin(), ParamNumber.end(), std::greater<int>());
1715
1716 //KS: Calculate how many pads/plots we need
1717 int Npad = 0;
1718 for(int j = 1; j < nParamPlot+1; j++) Npad += j;
1719 Posterior->cd();
1720 // KS: Sanity check of size and ordering is correct
1721 CheckCredibleIntervalsOrder(CredibleIntervals, CredibleIntervalsColours);
1722 CheckCredibleRegionsOrder(CredibleRegions, CredibleRegionStyle, CredibleRegionColor);
1723
1724 const int nCredibleIntervals = int(CredibleIntervals.size());
1725 const int nCredibleRegions = int(CredibleRegions.size());
1726
1727 //KS: Initialise Tpad histograms etc we will need
1728 std::vector<TPad*> TrianglePad(Npad);
1729 //KS: 1D copy of posterior, we need it as we modify them
1730 std::vector<std::unique_ptr<TH1D>> hpost_copy(nParamPlot);
1731 std::vector<std::vector<std::unique_ptr<TH1D>>> hpost_cl(nParamPlot);
1732 std::vector<std::unique_ptr<TText>> TriangleText(nParamPlot * 2);
1733 std::vector<std::unique_ptr<TH2D>> hpost_2D_copy(Npad-nParamPlot);
1734 std::vector<std::vector<std::unique_ptr<TH2D>>> hpost_2D_cl(Npad-nParamPlot);
1735 gStyle->SetPalette(51);
1736
1737 //KS: Super convoluted way of calculating ranges for our pads, trust me it works...
1738 std::vector<double> X_Min(nParamPlot);
1739 std::vector<double> X_Max(nParamPlot);
1740 X_Min[0] = 0.10;
1741 double xScale = (0.95 - (X_Min[0]+0.05))/nParamPlot;
1742 //KS: 0.05 is because we need additional offset for labels
1743 X_Max[0] = X_Min[0]+xScale+0.05;
1744 for(int i = 1; i < nParamPlot; i++)
1745 {
1746 X_Min[i] = X_Max[i-1];
1747 X_Max[i] = X_Min[i]+xScale;
1748 }
1749 std::vector<double> Y_Min(nParamPlot);
1750 std::vector<double> Y_Max(nParamPlot);
1751 Y_Max[0] = 0.95;
1752 //KS: 0.10 is becasue we need additional offset for labels
1753 double yScale = std::fabs(0.10 - (Y_Max[0]))/nParamPlot;
1754 Y_Min[0] = Y_Max[0]-yScale;
1755 for(int i = 1; i < nParamPlot; i++)
1756 {
1757 Y_Max[i] = Y_Min[i-1];
1758 Y_Min[i] = Y_Max[i]-yScale;
1759 }
1760
1761 //KS: We store as numbering of isn't straightforward
1762 int counterPad = 0, counterText = 0, counterPost = 0, counter2DPost = 0;
1763 //KS: We start from top of the plot, might be confusing but works very well
1764 for(int y = 0; y < nParamPlot; y++)
1765 {
1766 //KS: start from left and go right, depending on y
1767 for(int x = 0; x <= y; x++)
1768 {
1769 //KS: Need to go to canvas every time to have our pads in the same canvas, not pads in the pads
1770 Posterior->cd();
1771 TrianglePad[counterPad] = new TPad(Form("TPad_%i", counterPad), Form("TPad_%i", counterPad),
1772 X_Min[x], Y_Min[y], X_Max[x], Y_Max[y]);
1773
1774 TrianglePad[counterPad]->SetTopMargin(0);
1775 TrianglePad[counterPad]->SetRightMargin(0);
1776
1777 TrianglePad[counterPad]->SetGrid();
1778 TrianglePad[counterPad]->SetFrameBorderMode(0);
1779 TrianglePad[counterPad]->SetBorderMode(0);
1780 TrianglePad[counterPad]->SetBorderSize(0);
1781
1782 //KS: Corresponds to bottom part of the plot, need margins for labels
1783 TrianglePad[counterPad]->SetBottomMargin(y == (nParamPlot - 1) ? 0.1 : 0);
1784 //KS: Corresponds to left part, need margins for labels
1785 TrianglePad[counterPad]->SetLeftMargin(x == 0 ? 0.15 : 0);
1786
1787 TrianglePad[counterPad]->Draw();
1788 TrianglePad[counterPad]->cd();
1789
1790 //KS:if diagonal plot main posterior
1791 if(x == y)
1792 {
1793 hpost_copy[counterPost] = M3::Clone<TH1D>(hpost[ParamNumber[x]], Form("hpost_copy_%i", ParamNumber[x]));
1794 hpost_cl[counterPost].resize(nCredibleIntervals);
1796 hpost_copy[counterPost]->Scale(1. / hpost_copy[counterPost]->Integral());
1797 for (int j = 0; j < nCredibleIntervals; ++j)
1798 {
1799 hpost_cl[counterPost][j] = M3::Clone<TH1D>(hpost[ParamNumber[x]], Form("hpost_copy_%i_CL_%f", ParamNumber[x], CredibleIntervals[j]));
1800 //KS: Reset to get rid to TF1 otherwise we run into segfault :(
1801 hpost_cl[counterPost][j]->Reset("");
1802 hpost_cl[counterPost][j]->Fill(0.0, 0.0);
1803
1804 // Scale the histograms before gettindg credible intervals
1805 hpost_cl[counterPost][j]->Scale(1. / hpost_cl[counterPost][j]->Integral());
1806 GetCredibleIntervalSig(hpost_copy[counterPost], hpost_cl[counterPost][j], CredibleInSigmas, CredibleIntervals[j]);
1807
1808 hpost_cl[counterPost][j]->SetFillColor(CredibleIntervalsColours[j]);
1809 hpost_cl[counterPost][j]->SetLineWidth(1);
1810 }
1811
1812 hpost_copy[counterPost]->SetMaximum(hpost_copy[counterPost]->GetMaximum()*1.2);
1813 hpost_copy[counterPost]->SetLineWidth(2);
1814 hpost_copy[counterPost]->SetLineColor(kBlack);
1815
1816 //KS: Don't want any titles
1817 FormatHistogram(hpost_copy[counterPost]);
1818
1819 hpost_copy[counterPost]->Draw("HIST");
1820 for (int j = 0; j < nCredibleIntervals; ++j){
1821 hpost_cl[counterPost][j]->Draw("HIST SAME");
1822 }
1823 counterPost++;
1824 }
1825 //KS: Here we plot 2D credible regions
1826 else
1827 {
1828 hpost_2D_copy[counter2DPost] = M3::Clone<TH2D>(hpost2D[ParamNumber[x]][ParamNumber[y]],
1829 Form("hpost_copy_%i_%i", ParamNumber[x], ParamNumber[y]));
1830 hpost_2D_cl[counter2DPost].resize(nCredibleRegions);
1831 //KS: Now copy for every credible region
1832 for (int k = 0; k < nCredibleRegions; ++k)
1833 {
1834 hpost_2D_cl[counter2DPost][k] = M3::Clone<TH2D>(hpost2D[ParamNumber[x]][ParamNumber[y]],
1835 Form("hpost_copy_%i_%i_CL_%f", ParamNumber[x], ParamNumber[y], CredibleRegions[k]));
1836 GetCredibleRegionSig(hpost_2D_cl[counter2DPost][k], CredibleInSigmas, CredibleRegions[k]);
1837
1838 hpost_2D_cl[counter2DPost][k]->SetLineColor(CredibleRegionColor[k]);
1839 hpost_2D_cl[counter2DPost][k]->SetLineWidth(2);
1840 hpost_2D_cl[counter2DPost][k]->SetLineStyle(CredibleRegionStyle[k]);
1841 }
1842 //KS: Don't want any titles
1843 FormatHistogram(hpost_2D_copy[counter2DPost]);
1844
1845 hpost_2D_copy[counter2DPost]->Draw("COL");
1846 //Now credible regions
1847 for (int k = 0; k < nCredibleRegions; ++k){
1848 hpost_2D_cl[counter2DPost][k]->Draw("CONT3 SAME");
1849 }
1850 counter2DPost++;
1851 }
1852 //KS: Corresponds to bottom part of the plot
1853 if(y == (nParamPlot-1))
1854 {
1855 Posterior->cd();
1856 TriangleText[counterText] = std::make_unique<TText>(X_Min[x]+ (X_Max[x]-X_Min[x])/4, 0.04, hpost[ParamNumber[x]]->GetTitle());
1857 //KS: Unfortunately for many plots or long names this can go out of bounds :(
1858 TriangleText[counterText]->SetTextSize(0.015);
1859 TriangleText[counterText]->SetNDC(true);
1860 TriangleText[counterText]->Draw();
1861 counterText++;
1862 }
1863 //KS: Corresponds to left part
1864 if(x == 0)
1865 {
1866 Posterior->cd();
1867 TriangleText[counterText] = std::make_unique<TText>(0.04, Y_Min[y] + (Y_Max[y]-Y_Min[y])/4, hpost[ParamNumber[y]]->GetTitle());
1868 //KS: Rotate as this is y axis
1869 TriangleText[counterText]->SetTextAngle(90);
1870 //KS: Unfortunately for many plots or long names this can go out of bounds :(
1871 TriangleText[counterText]->SetTextSize(0.015);
1872 TriangleText[counterText]->SetNDC(true);
1873 TriangleText[counterText]->Draw();
1874 counterText++;
1875 }
1876 Posterior->Update();
1877 counterPad++;
1878 }
1879 }
1880
1881 Posterior->cd();
1882 auto legend = std::make_unique<TLegend>(0.60, 0.7, 0.9, 0.9);
1883 SetLegendStyle(legend.get(), 0.03);
1884 //KS: Legend is shared so just take first histograms
1885 for (int j = nCredibleIntervals-1; j >= 0; --j)
1886 {
1887 if(CredibleInSigmas)
1888 legend->AddEntry(hpost_cl[0][j].get(), Form("%.0f#sigma Credible Interval", CredibleIntervals[j]), "f");
1889 else
1890 legend->AddEntry(hpost_cl[0][j].get(), Form("%.0f%% Credible Interval", CredibleRegions[j]*100), "f");
1891 }
1892 for (int k = nCredibleRegions-1; k >= 0; --k)
1893 {
1894 if(CredibleInSigmas)
1895 legend->AddEntry(hpost_2D_cl[0][k].get(), Form("%.0f#sigma Credible Region", CredibleRegions[k]), "l");
1896 else
1897 legend->AddEntry(hpost_2D_cl[0][k].get(), Form("%.0f%% Credible Region", CredibleRegions[k]*100), "l");
1898 }
1899 legend->Draw("SAME");
1900 Posterior->Update();
1901
1902 // Write to file
1903 Posterior->SetName("TrianglePlot");
1904 Posterior->SetTitle("TrianglePlot");
1905
1906 if(printToPDF) Posterior->Print(CanvasName);
1907 // Write it to root file
1908 OutputFile->cd();
1909 Posterior->Write();
1910
1911 //KS: Remove allocated structures
1912 for(int i = 0; i < Npad; i++) delete TrianglePad[i];
1913
1914 //KS: Restore margin
1915 SetMargins(Posterior, Margins);
1916}

◆ MakeViolin()

void MCMCProcessor::MakeViolin ( )

Make and Draw Violin.

Definition at line 781 of file MCMCProcessor.cpp.

781 {
782// *********************
783 //KS: Make sure we have steps
784 if(!CacheMCMC) CacheSteps();
785
786 MACH3LOG_INFO("Producing Violin Plot");
787
788 //KS: Find min and max to make histogram in range
789 double maxi_y = Chain->GetMaximum(BranchNames[0]);
790 double mini_y = Chain->GetMinimum(BranchNames[0]);
791 for (int i = 1; i < nDraw; ++i)
792 {
793 const double max_val = Chain->GetMaximum(BranchNames[i]);
794 const double min_val = Chain->GetMinimum(BranchNames[i]);
795
796 maxi_y = std::max(maxi_y, max_val);
797 mini_y = std::min(mini_y, min_val);
798 }
799
800 const int vBins = (maxi_y-mini_y)*25;
801 hviolin = std::make_unique<TH2D>("hviolin", "hviolin", nDraw, 0, nDraw, vBins, mini_y, maxi_y);
802 hviolin->SetDirectory(nullptr);
803 //KS: Prior has larger errors so we increase range and number of bins
804 constexpr int PriorFactor = 4;
805 hviolin_prior = std::make_unique<TH2D>("hviolin_prior", "hviolin_prior", nDraw, 0, nDraw, PriorFactor*vBins, PriorFactor*mini_y, PriorFactor*maxi_y);
806 hviolin_prior->SetDirectory(nullptr);
807
808 auto rand = std::make_unique<TRandom3>(0);
809 std::vector<double> PriorVec(nDraw);
810 std::vector<double> PriorErrorVec(nDraw);
811 std::vector<bool> PriorFlatVec(nDraw);
812
813 for (int x = 0; x < nDraw; ++x)
814 {
815 TString Title;
816 double Prior, PriorError;
817
818 GetNthParameter(x, Prior, PriorError, Title);
819 //Set fancy labels
820 hviolin->GetXaxis()->SetBinLabel(x+1, Title);
821 hviolin_prior->GetXaxis()->SetBinLabel(x+1, Title);
822 PriorVec[x] = Prior;
823 PriorErrorVec[x] = PriorError;
824
825 ParameterEnum ParType = ParamType[x];
826 int ParamTemp = x - ParamTypeStartPos[ParType];
827 PriorFlatVec[x] = ParamFlat[ParType][ParamTemp];
828 }
829
830 TStopwatch clock;
831 clock.Start();
832
833 // nDraw is number of draws we want to do
834 #ifdef MULTITHREAD
835 #pragma omp parallel for
836 #endif
837 for (int x = 0; x < nDraw; ++x)
838 {
839 //KS: Consider another treatment for fixed params
840 //if (IamVaried[x] == false) continue;
841 for (int k = 0; k < nEntries; ++k)
842 {
843 //KS: Burn in cut
844 if(StepNumber[k] < BurnInCut) continue;
845
846 //Only used for Suboptimatlity
847 if(StepNumber[k] > UpperCut) continue;
848
849 //KS: We know exactly which x bin we will end up, find y bin. This allow to avoid coslty Fill() and enable multithreading becasue I am master of faster
850 const double y = hviolin->GetYaxis()->FindBin(ParStep[x][k]);
851 hviolin->SetBinContent(x+1, y, hviolin->GetBinContent(x+1, y)+1);
852 }
853
854 //KS: If we set option to not plot flat prior and param has flat prior then we skip this step
855 if(!(!PlotFlatPrior && PriorFlatVec[x]))
856 {
857 for (int k = 0; k < nEntries; ++k)
858 {
859 const double Entry = rand->Gaus(PriorVec[x], PriorErrorVec[x]);
860 const double y = hviolin_prior->GetYaxis()->FindBin(Entry);
861 hviolin_prior->SetBinContent(x+1, y, hviolin_prior->GetBinContent(x+1, y)+1);
862 }
863 }
864 } // end the for loop over nDraw
865 clock.Stop();
866 MACH3LOG_INFO("Making Violin plot took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries);
867
868 //KS: Tells how many parameters in one canvas we want
869 constexpr int IntervalsSize = 10;
870 const int NIntervals = nDraw/IntervalsSize;
871
872 hviolin->GetYaxis()->SetTitle("Parameter Value");
873 hviolin->GetXaxis()->SetTitle();
874 hviolin->GetXaxis()->LabelsOption("v");
875
876 hviolin_prior->GetYaxis()->SetTitle("Parameter Value");
877 hviolin_prior->GetXaxis()->SetTitle();
878 hviolin_prior->GetXaxis()->LabelsOption("v");
879
880 hviolin_prior->SetLineColor(kRed);
881 hviolin_prior->SetMarkerColor(kRed);
882 hviolin_prior->SetFillColorAlpha(kRed, 0.35);
883 hviolin_prior->SetMarkerStyle(20);
884 hviolin_prior->SetMarkerSize(0.5);
885
886 // These control violin width, if you use larger then 1 they will most likely overlay, so be cautious
887 hviolin_prior->SetBarWidth(1.0);
888 hviolin_prior->SetBarOffset(0);
889
890 hviolin->SetLineColor(kBlue);
891 hviolin->SetMarkerColor(kBlue);
892 hviolin->SetFillColorAlpha(kBlue, 0.35);
893 hviolin->SetMarkerStyle(20);
894 hviolin->SetMarkerSize(1.0);
895
896 const double BottomMargin = Posterior->GetBottomMargin();
897 Posterior->SetBottomMargin(0.2);
898
899 OutputFile->cd();
900 hviolin->Write("param_violin");
901 hviolin_prior->Write("param_violin_prior");
902 //KS: This is mostly for example plots, we have full file in the ROOT file so can do much better plot later
903 hviolin->GetYaxis()->SetRangeUser(-1, +2);
904 hviolin_prior->GetYaxis()->SetRangeUser(-1, +2);
905 for (int i = 0; i < NIntervals+1; ++i)
906 {
907 hviolin->GetXaxis()->SetRangeUser(i*IntervalsSize, i*IntervalsSize+IntervalsSize);
908 hviolin_prior->GetXaxis()->SetRangeUser(i*IntervalsSize, i*IntervalsSize+IntervalsSize);
909 if(i == NIntervals+1)
910 {
911 hviolin->GetXaxis()->SetRangeUser(i*IntervalsSize, nDraw);
912 hviolin_prior->GetXaxis()->SetRangeUser(i*IntervalsSize, nDraw);
913 }
914 //KS: ROOT6 has some additional options, consider updating it. more https://root.cern/doc/master/classTHistPainter.html#HP140b
915 hviolin_prior->Draw("violinX(03100300)");
916 hviolin->Draw("violinX(03100300) SAME");
917 if(printToPDF) Posterior->Print(CanvasName);
918 }
919 //KS: Return Margin to default one
920 Posterior->SetBottomMargin(BottomMargin);
921}

◆ ParameterEvolution()

void MCMCProcessor::ParameterEvolution ( const std::vector< std::string > &  Names,
const std::vector< int > &  NIntervals 
)

Make .gif of parameter evolution.

Parameters
NamesParameter names for which we do .gif
NIntervalsNumber of intervals for a gif

Definition at line 2836 of file MCMCProcessor.cpp.

2837 {
2838// **************************
2839 MACH3LOG_INFO("Parameter Evolution gif");
2840
2841 //KS: First we need to find parameter number based on name
2842 for(unsigned int k = 0; k < Names.size(); ++k)
2843 {
2844 //KS: First we need to find parameter number based on name
2845 int ParamNo = GetParamIndexFromName(Names[k]);
2846 if(ParamNo == _UNDEF_)
2847 {
2848 MACH3LOG_WARN("Couldn't find param {}. Can't reweight Prior", Names[k]);
2849 continue;
2850 }
2851
2852 const int IntervalsSize = nSteps/NIntervals[k];
2853 // ROOT won't overwrite gifs so we need to delete the file if it's there already
2854 int ret = system(fmt::format("rm {}.gif",Names[k]).c_str());
2855 if (ret != 0){
2856 MACH3LOG_WARN("Error: system call to delete {} failed with code {}", Names[k], ret);
2857 }
2858
2859 // This holds the posterior density
2860 const double maxi = Chain->GetMaximum(BranchNames[ParamNo]);
2861 const double mini = Chain->GetMinimum(BranchNames[ParamNo]);
2862
2863 int Counter = 0;
2864 for(int i = NIntervals[k]-1; i >= 0; --i)
2865 {
2866 // This holds the posterior density
2867 TH1D* EvePlot = new TH1D(BranchNames[ParamNo], BranchNames[ParamNo], nBins, mini, maxi);
2868 EvePlot->SetMinimum(0);
2869 EvePlot->GetYaxis()->SetTitle("PDF");
2870 EvePlot->GetYaxis()->SetNoExponent(false);
2871
2872 //KS: Apply additional Cuts, like mass ordering
2873 std::string CutPosterior1D = "step > " + std::to_string(i*IntervalsSize+IntervalsSize);
2874
2875 std::string TextTitle = "Steps = 0 - "+std::to_string(Counter*IntervalsSize+IntervalsSize);
2876 // Project BranchNames[ParamNo] onto hpost, applying stepcut
2877 Chain->Project(BranchNames[ParamNo], BranchNames[ParamNo], CutPosterior1D.c_str());
2878
2879 EvePlot->SetLineWidth(2);
2880 EvePlot->SetLineColor(kBlue-1);
2881 EvePlot->SetTitle(Names[k].c_str());
2882 EvePlot->GetXaxis()->SetTitle(EvePlot->GetTitle());
2883 EvePlot->GetYaxis()->SetLabelOffset(1000);
2884 if(ApplySmoothing) EvePlot->Smooth();
2885
2886 EvePlot->Scale(1. / EvePlot->Integral());
2887 EvePlot->Draw("HIST");
2888
2889 TText text(0.3, 0.8, TextTitle.c_str());
2890 text.SetTextFont (43);
2891 text.SetTextSize (40);
2892 text.SetNDC(true);
2893 text.Draw("SAME");
2894
2895 if(i == 0) Posterior->Print((Names[k] + ".gif++20").c_str()); // produces infinite loop animated GIF
2896 else Posterior->Print((Names[k] + ".gif+20").c_str()); // add picture to .gif
2897
2898 delete EvePlot;
2899 Counter++;
2900 }
2901 }
2902}

◆ ParamTraces()

void MCMCProcessor::ParamTraces ( )
inlineprotected

CW: Draw trace plots of the parameters i.e. parameter vs step.

Definition at line 3099 of file MCMCProcessor.cpp.

3099 {
3100// *****************
3101 if (ParStep == nullptr) PrepareDiagMCMC();
3102 MACH3LOG_INFO("Making trace plots...");
3103 // Make the TH1Ds
3104 std::vector<TH1D*> TraceParamPlots(nDraw);
3105 std::vector<TH1D*> TraceSamplePlots(nSamples);
3106 std::vector<TH1D*> TraceSystsPlots(nSysts);
3107
3108 // Set the titles and limits for TH2Ds
3109 for (int j = 0; j < nDraw; ++j) {
3110 TString Title = "";
3111 double Prior = 1.0, PriorError = 1.0;
3112
3113 GetNthParameter(j, Prior, PriorError, Title);
3114 std::string HistName = Form("%s_%s_Trace", Title.Data(), BranchNames[j].Data());
3115 TraceParamPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nEntries, 0, nEntries);
3116 TraceParamPlots[j]->GetXaxis()->SetTitle("Step");
3117 TraceParamPlots[j]->GetYaxis()->SetTitle("Parameter Variation");
3118 }
3119
3120 for (int j = 0; j < nSamples; ++j) {
3121 std::string HistName = SampleName_v[j].Data();
3122 TraceSamplePlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nEntries, 0, nEntries);
3123 TraceSamplePlots[j]->GetXaxis()->SetTitle("Step");
3124 TraceSamplePlots[j]->GetYaxis()->SetTitle("Sample -logL");
3125 }
3126
3127 for (int j = 0; j < nSysts; ++j) {
3128 std::string HistName = SystName_v[j].Data();
3129 TraceSystsPlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), nEntries, 0, nEntries);
3130 TraceSystsPlots[j]->GetXaxis()->SetTitle("Step");
3131 TraceSystsPlots[j]->GetYaxis()->SetTitle("Systematic -logL");
3132 }
3133
3134 // Have now made the empty TH1Ds, now for writing content to them!
3135 // Loop over the number of parameters to draw their traces
3136 // Each histogram
3137#ifdef MULTITHREAD
3138 MACH3LOG_INFO("Using multi-threading...");
3139 #pragma omp parallel for
3140#endif
3141 for (int i = 0; i < nEntries; ++i) {
3142 // Set bin content for the ith bin to the parameter values
3143 for (int j = 0; j < nDraw; ++j) {
3144 TraceParamPlots[j]->SetBinContent(i, ParStep[j][i]);
3145 }
3146 for (int j = 0; j < nSamples; ++j) {
3147 TraceSamplePlots[j]->SetBinContent(i, SampleValues[i][j]);
3148 }
3149 for (int j = 0; j < nSysts; ++j) {
3150 TraceSystsPlots[j]->SetBinContent(i, SystValues[i][j]);
3151 }
3152 }
3153
3154 // Write the output and delete the TH2Ds
3155 TDirectory *TraceDir = OutputFile->mkdir("Trace");
3156 TraceDir->cd();
3157 for (int j = 0; j < nDraw; ++j) {
3158 // Fit a linear function to the traces
3159 TF1 *Fitter = new TF1("Fitter","[0]", nEntries/2, nEntries);
3160 Fitter->SetLineColor(kRed);
3161 TraceParamPlots[j]->Fit("Fitter","Rq");
3162 TraceParamPlots[j]->Write();
3163 delete Fitter;
3164 delete TraceParamPlots[j];
3165 }
3166
3167 TDirectory *LLDir = OutputFile->mkdir("LogL");
3168 LLDir->cd();
3169 for (int j = 0; j < nSamples; ++j) {
3170 TraceSamplePlots[j]->Write();
3171 delete TraceSamplePlots[j];
3172 delete[] SampleValues[j];
3173 }
3174 delete[] SampleValues;
3175
3176 for (int j = 0; j < nSysts; ++j) {
3177 TraceSystsPlots[j]->Write();
3178 delete TraceSystsPlots[j];
3179 delete SystValues[j];
3180 }
3181 delete[] SystValues;
3182
3183 TraceDir->Close();
3184 delete TraceDir;
3185
3186 OutputFile->cd();
3187}
MCMCProcessor::SampleName_v
std::vector< TString > SampleName_v
Vector of each systematic.
Definition: MCMCProcessor.h:420
MCMCProcessor::SystName_v
std::vector< TString > SystName_v
Vector of each sample PDF object.
Definition: MCMCProcessor.h:422

◆ PowerSpectrumAnalysis()

void MCMCProcessor::PowerSpectrumAnalysis ( )
inlineprotected

RC: Perform spectral analysis of MCMC [7].

Author
Richard Calland
Todo:
KS: Code is awfully slow... I know how to make it faster (GPU scream in a distant) but for now just make it for two params, bit hacky sry...

Definition at line 3776 of file MCMCProcessor.cpp.

3776 {
3777// **************************
3778 TStopwatch clock;
3779 clock.Start();
3780
3781 //KS: Store it as we go back to them at the end
3782 const double TopMargin = Posterior->GetTopMargin();
3783 const int OptTitle = gStyle->GetOptTitle();
3784
3785 Posterior->SetTopMargin(0.1);
3786 gStyle->SetOptTitle(1);
3787
3788 MACH3LOG_INFO("Making Power Spectrum plots...");
3789
3790 // This is only to reduce number of computations...
3791 const int N_Coeffs = std::min(10000, nEntries);
3792 const int start = -(N_Coeffs/2-1);
3793 const int end = N_Coeffs/2-1;
3794 const int v_size = end - start;
3795
3796 int nPrams = nDraw;
3798 nPrams = 2;
3799
3800 std::vector<std::vector<float>> k_j(nPrams, std::vector<float>(v_size, 0.0));
3801 std::vector<std::vector<float>> P_j(nPrams, std::vector<float>(v_size, 0.0));
3802
3803 int _N = nEntries;
3804 if (_N % 2 != 0) _N -= 1; // N must be even
3805
3806 //This is being used a lot so calculate it once to increase performance
3807 const double two_pi_over_N = 2 * TMath::Pi() / static_cast<double>(_N);
3808
3809 // KS: This could be moved to GPU I guess
3810 #ifdef MULTITHREAD
3811 #pragma omp parallel for collapse(2)
3812 #endif
3813 // RC: equation 11: for each value of j coef, from range -N/2 -> N/2
3814 for (int j = 0; j < nPrams; ++j)
3815 {
3816 for (int jj = start; jj < end; ++jj)
3817 {
3818 std::complex<double> a_j = 0.0;
3819 for (int n = 0; n < _N; ++n)
3820 {
3821 //if(StepNumber[n] < BurnInCut) continue;
3822 std::complex<double> exp_temp(0, two_pi_over_N * jj * n);
3823 a_j += ParStep[j][n] * std::exp(exp_temp);
3824 }
3825 a_j /= std::sqrt(float(_N));
3826 const int _c = jj - start;
3827
3828 k_j[j][_c] = two_pi_over_N * jj;
3829 // Equation 13
3830 P_j[j][_c] = std::norm(a_j);
3831 }
3832 }
3833
3834 TDirectory *PowerDir = OutputFile->mkdir("PowerSpectrum");
3835 PowerDir->cd();
3836
3837 TVectorD* PowerSpectrumStepSize = new TVectorD(nPrams);
3838 for (int j = 0; j < nPrams; ++j)
3839 {
3840 auto plot = std::make_unique<TGraph>(v_size, k_j[j].data(), P_j[j].data());
3841
3842 TString Title = "";
3843 double Prior = 1.0, PriorError = 1.0;
3844 GetNthParameter(j, Prior, PriorError, Title);
3845
3846 std::string name = Form("Power Spectrum of %s;k;P(k)", Title.Data());
3847
3848 plot->SetTitle(name.c_str());
3849 name = Form("%s_power_spectrum", Title.Data());
3850 plot->SetName(name.c_str());
3851 plot->SetMarkerStyle(7);
3852
3853 // Equation 18
3854 TF1 *func = new TF1("power_template", "[0]*( ([1] / x)^[2] / (([1] / x)^[2] +1) )", 0.0, 1.0);
3855 // P0 gives the amplitude of the white noise spectrum in the k → 0 limit
3856 func->SetParameter(0, 10.0);
3857 // k* indicates the position of the turnover to a different power law behaviour
3858 func->SetParameter(1, 0.1);
3859 // alpha free parameter
3860 func->SetParameter(2, 2.0);
3861
3862 // Set parameter limits for stability
3863 func->SetParLimits(0, 0.0, 100.0); // Amplitude should be non-negative
3864 func->SetParLimits(1, 0.001, 1.0); // k* should be within a reasonable range
3865 func->SetParLimits(2, 0.0, 5.0); // alpha should be positive
3866
3867 plot->Fit("power_template","Rq");
3868
3869 Posterior->SetLogx();
3870 Posterior->SetLogy();
3871 Posterior->SetGrid();
3872 plot->Write(plot->GetName());
3873 plot->Draw("AL");
3874 func->Draw("SAME");
3875 if(printToPDF) Posterior->Print(CanvasName);
3876
3877 //KS: I have no clue what is the reason behind this. Found this in Rick Calland code...
3878 (*PowerSpectrumStepSize)(j) = std::sqrt(func->GetParameter(0)/float(v_size*0.5));
3879 delete func;
3880 }
3881
3882 PowerSpectrumStepSize->Write("PowerSpectrumStepSize");
3883 delete PowerSpectrumStepSize;
3884 PowerDir->Close();
3885 delete PowerDir;
3886
3887 clock.Stop();
3888 MACH3LOG_INFO("Making Power Spectrum took {:.2f}s", clock.RealTime());
3889
3890 Posterior->SetTopMargin(TopMargin);
3891 gStyle->SetOptTitle(OptTitle);
3892}

◆ PrepareDiagMCMC()

void MCMCProcessor::PrepareDiagMCMC ( )
inlineprotected

CW: Prepare branches etc. for DiagMCMC.

Definition at line 2936 of file MCMCProcessor.cpp.

2936 {
2937// **************************
2938 doDiagMCMC = true;
2939
2940 if(ParStep != nullptr) {
2941 MACH3LOG_ERROR("It look like ParStep was already filled ");
2942 MACH3LOG_ERROR("Even though it is used for MakeCovariance_MP and for DiagMCMC");
2943 MACH3LOG_ERROR("it has different structure in both for cache hits, sorry ");
2944 throw MaCh3Exception(__FILE__ , __LINE__ );
2945 }
2946 if(nBatches == 0) {
2947 MACH3LOG_ERROR("nBatches is equal to 0");
2948 MACH3LOG_ERROR("please use SetnBatches to set other value fore example 20");
2949 throw MaCh3Exception(__FILE__ , __LINE__ );
2950 }
2951
2952 // Initialise ParStep
2953 ParStep = new double*[nDraw]();
2954 for (int j = 0; j < nDraw; ++j) {
2955 ParStep[j] = new double[nEntries]();
2956 for (int i = 0; i < nEntries; ++i) {
2957 ParStep[j][i] = -999.99;
2958 }
2959 }
2960
2961 SampleValues = new double*[nEntries]();
2962 SystValues = new double*[nEntries]();
2963 AccProbValues = new double[nEntries]();
2964 StepNumber = new int[nEntries]();
2965 for (int i = 0; i < nEntries; ++i) {
2966 SampleValues[i] = new double[nSamples]();
2967 SystValues[i] = new double[nSysts]();
2968
2969 for (int j = 0; j < nSamples; ++j) {
2970 SampleValues[i][j] = -999.99;
2971 }
2972 for (int j = 0; j < nSysts; ++j) {
2973 SystValues[i][j] = -999.99;
2974 }
2975 AccProbValues[i] = -999.99;
2976 StepNumber[i] = -999.99;
2977 }
2978
2979 // Initialise the sums
2980 ParamSums = new double[nDraw]();
2981 for (int i = 0; i < nDraw; ++i) {
2982 ParamSums[i] = 0.0;
2983 }
2984 MACH3LOG_INFO("Reading input tree...");
2985 TStopwatch clock;
2986 clock.Start();
2987
2988 // Set all the branches to off
2989 Chain->SetBranchStatus("*", false);
2990
2991 // 10 entries output
2992 const int countwidth = nEntries/10;
2993
2994 // Can also do the batched means here to minimize excessive loops
2995 // The length of each batch
2996 const int BatchLength = nEntries/nBatches+1;
2997 BatchedAverages = new double*[nBatches]();
2998 AccProbBatchedAverages = new double[nBatches]();
2999 for (int i = 0; i < nBatches; ++i) {
3000 BatchedAverages[i] = new double[nDraw];
3001 AccProbBatchedAverages[i] = 0;
3002 for (int j = 0; j < nDraw; ++j) {
3003 BatchedAverages[i][j] = 0.0;
3004 }
3005 }
3006 std::vector<double> ParStepBranch(nDraw);
3007 std::vector<double> SampleValuesBranch(nSamples);
3008 std::vector<double> SystValuesBranch(nSysts);
3009 int StepNumberBranch = 0;
3010 double AccProbValuesBranch = 0;
3011 // Set the branch addresses for params
3012 for (int j = 0; j < nDraw; ++j) {
3013 Chain->SetBranchStatus(BranchNames[j].Data(), true);
3014 Chain->SetBranchAddress(BranchNames[j].Data(), &ParStepBranch[j]);
3015 }
3016 // Set the branch addresses for samples
3017 for (int j = 0; j < nSamples; ++j) {
3018 Chain->SetBranchStatus(SampleName_v[j].Data(), true);
3019 Chain->SetBranchAddress(SampleName_v[j].Data(), &SampleValuesBranch[j]);
3020 }
3021 // Set the branch addresses for systematics
3022 for (int j = 0; j < nSysts; ++j) {
3023 Chain->SetBranchStatus(SystName_v[j].Data(), true);
3024 Chain->SetBranchAddress(SystName_v[j].Data(), &SystValuesBranch[j]);
3025 }
3026 // Only needed for Geweke right now
3027 Chain->SetBranchStatus("step", true);
3028 Chain->SetBranchAddress("step", &StepNumberBranch);
3029 // Turn on the branches which we want for acc prob
3030 Chain->SetBranchStatus("accProb", true);
3031 Chain->SetBranchAddress("accProb", &AccProbValuesBranch);
3032
3033 // Loop over the entries
3034 //KS: This is really a bottleneck right now, thus revisit with ROOT6 https://pep-root6.github.io/docs/analysis/parallell/root.html
3035 for (int i = 0; i < nEntries; ++i) {
3036 // Fill up the arrays
3037 Chain->GetEntry(i);
3038
3039 if (i % countwidth == 0)
3040 MaCh3Utils::PrintProgressBar(i, nEntries);
3041
3042 // Set the branch addresses for params
3043 for (int j = 0; j < nDraw; ++j) {
3044 ParStep[j][i] = ParStepBranch[j];
3045 }
3046 // Set the branch addresses for samples
3047 for (int j = 0; j < nSamples; ++j) {
3048 SampleValues[i][j] = SampleValuesBranch[j];
3049 }
3050 // Set the branch addresses for systematics
3051 for (int j = 0; j < nSysts; ++j) {
3052 SystValues[i][j] = SystValuesBranch[j];
3053 }
3054
3055 // Set the branch addresses for Acceptance Probability
3056 AccProbValues[i] = AccProbValuesBranch;
3057 StepNumber[i] = StepNumberBranch;
3058
3059 // Find which batch the event belongs in
3060 int BatchNumber = -1;
3061 // I'm so lazy! But it's OK, the major overhead here is GetEntry: saved by ROOT!
3062 for (int j = 0; j < nBatches; ++j) {
3063 if (i < (j+1)*BatchLength) {
3064 BatchNumber = j;
3065 break;
3066 }
3067 }
3068 // Fill up the sum for each j param
3069 for (int j = 0; j < nDraw; ++j) {
3070 ParamSums[j] += ParStep[j][i];
3071 BatchedAverages[BatchNumber][j] += ParStep[j][i];
3072 }
3073
3074 //KS: Could easily add this to above loop but I accProb is different beast so better keep it like this
3075 AccProbBatchedAverages[BatchNumber] += AccProbValues[i];
3076 }
3077 clock.Stop();
3078 MACH3LOG_INFO("Took {:.2f}s to finish caching statistic for Diag MCMC with {} steps", clock.RealTime(), nEntries);
3079
3080 // Make the sums into average
3081 #ifdef MULTITHREAD
3082 #pragma omp parallel for
3083 #endif
3084 for (int i = 0; i < nDraw; ++i) {
3085 ParamSums[i] /= double(nEntries);
3086 for (int j = 0; j < nBatches; ++j) {
3087 // Divide by the total number of events in the batch
3088 BatchedAverages[j][i] /= BatchLength;
3089 if(i == 0) AccProbBatchedAverages[j] /= BatchLength; //KS: we have only one accProb, keep it like this for now
3090 }
3091 }
3092
3093 // And make our sweet output file
3094 if (OutputFile == nullptr) MakeOutputFile();
3095}

◆ PrintInfo()

void MCMCProcessor::PrintInfo ( ) const
inlineprotected

Print info like how many params have been loaded etc.

Definition at line 4131 of file MCMCProcessor.cpp.

4131 {
4132// **************************
4133 // KS: Create a map to store the counts of unique strings
4134 std::unordered_map<std::string, int> paramCounts;
4135 std::vector<std::string> orderedKeys;
4136
4137 for (const std::string& param : ParameterGroup) {
4138 if (paramCounts[param] == 0) {
4139 orderedKeys.push_back(param); // preserve order of first appearance
4140 }
4141 paramCounts[param]++;
4142 }
4143
4144 MACH3LOG_INFO("************************************************");
4145 MACH3LOG_INFO("Scanning output branches...");
4146 MACH3LOG_INFO("# Useful entries in tree: \033[1;32m {} \033[0m ", nDraw);
4147 MACH3LOG_INFO("# Model params: \033[1;32m {} starting at {} \033[0m ", nParam[kXSecPar], ParamTypeStartPos[kXSecPar]);
4148 MACH3LOG_INFO("# With following groups: ");
4149 for (const std::string& key : orderedKeys) {
4150 MACH3LOG_INFO(" # {} params: {}", key, paramCounts[key]);
4151 }
4152 MACH3LOG_INFO("# ND params (legacy): \033[1;32m {} starting at {} \033[0m ", nParam[kNDPar], ParamTypeStartPos[kNDPar]);
4153 MACH3LOG_INFO("# FD params (legacy): \033[1;32m {} starting at {} \033[0m ", nParam[kFDDetPar], ParamTypeStartPos[kFDDetPar]);
4154 MACH3LOG_INFO("************************************************");
4155}

◆ ReadFDFile()

void MCMCProcessor::ReadFDFile ( )
inlineprotected

Read the FD cov file and get the input central values and errors.

Definition at line 2373 of file MCMCProcessor.cpp.

2373 {
2374// ***************
2375 // Do the same for the FD
2376 TFile *FDdetFile = new TFile(CovPos[kFDDetPar].back().c_str(), "open");
2377 if (FDdetFile->IsZombie()) {
2378 MACH3LOG_ERROR("Couldn't find NDdetFile {}", CovPos[kFDDetPar].back());
2379 throw MaCh3Exception(__FILE__ , __LINE__ );
2380 }
2381 FDdetFile->cd();
2382
2383 TMatrixDSym *FDdetMatrix = FDdetFile->Get<TMatrixDSym>("SKJointError_Erec_Total");
2384
2385 for (int i = 0; i < FDdetMatrix->GetNrows(); ++i)
2386 {
2387 //KS: FD parameters start at 1. in contrary to ND280
2388 ParamNom[kFDDetPar].push_back(1.);
2389 ParamCentral[kFDDetPar].push_back(1.);
2390
2391 ParamErrors[kFDDetPar].push_back( std::sqrt((*FDdetMatrix)(i,i)) );
2392 ParamNames[kFDDetPar].push_back( Form("FD Det %i", i) );
2393
2394 //KS: Currently we can only set it via config, change it in future
2395 ParamFlat[kFDDetPar].push_back( false );
2396 }
2397 //KS: The last parameter is p scale
2398 if(FancyPlotNames) ParamNames[kFDDetPar].back() = "Momentum Scale";
2399
2400 FDdetFile->Close();
2401 delete FDdetFile;
2402 delete FDdetMatrix;
2403}

◆ ReadInputCov()

void MCMCProcessor::ReadInputCov ( )
inlineprotected

CW: Read the input Covariance matrix entries. Get stuff like parameter input errors, names, and so on.

Definition at line 2162 of file MCMCProcessor.cpp.

2162 {
2163// **************************
2164 FindInputFiles();
2165 if(CovPos[kXSecPar].back() != "none") ReadModelFile();
2166}
MCMCProcessor::ReadModelFile
void ReadModelFile()
Read the xsec file and get the input central values and errors.
Definition: MCMCProcessor.cpp:2280
MCMCProcessor::FindInputFiles
void FindInputFiles()
Read the output MCMC file and find what inputs were used.
Definition: MCMCProcessor.cpp:2180

◆ ReadInputCovLegacy()

void MCMCProcessor::ReadInputCovLegacy ( )
inlineprotected
Warning
This will no longer be supported in future

Definition at line 2171 of file MCMCProcessor.cpp.

2171 {
2172// **************************
2173 FindInputFilesLegacy();
2174 if(nParam[kNDPar] > 0) ReadNDFile();
2175 if(nParam[kFDDetPar] > 0) ReadFDFile();
2176}
MCMCProcessor::ReadNDFile
void ReadNDFile()
Read the ND cov file and get the input central values and errors.
Definition: MCMCProcessor.cpp:2330
MCMCProcessor::FindInputFilesLegacy
void FindInputFilesLegacy()
Definition: MCMCProcessor.cpp:2239
MCMCProcessor::ReadFDFile
void ReadFDFile()
Read the FD cov file and get the input central values and errors.
Definition: MCMCProcessor.cpp:2373

◆ ReadModelFile()

void MCMCProcessor::ReadModelFile ( )
inlineprotected

Read the xsec file and get the input central values and errors.

Definition at line 2280 of file MCMCProcessor.cpp.

2280 {
2281// ***************
2282 YAML::Node XSecFile = CovConfig[kXSecPar];
2283
2284 auto systematics = XSecFile["Systematics"];
2285 int paramIndex = 0;
2286 for (auto it = systematics.begin(); it != systematics.end(); ++it, ++paramIndex )
2287 {
2288 auto const &param = *it;
2289 // Push back the name
2290 std::string TempString = (param["Systematic"]["Names"]["FancyName"].as<std::string>());
2291
2292 bool rejected = false;
2293 for (unsigned int ik = 0; ik < ExcludedNames.size(); ++ik)
2294 {
2295 if (TempString.rfind(ExcludedNames.at(ik), 0) == 0)
2296 {
2297 rejected = true;
2298 break;
2299 }
2300 }
2301 if(rejected) continue;
2302
2303 ParamNames[kXSecPar].push_back(TempString);
2304 ParamCentral[kXSecPar].push_back(param["Systematic"]["ParameterValues"]["PreFitValue"].as<double>());
2305 ParamNom[kXSecPar].push_back(param["Systematic"]["ParameterValues"]["Generated"].as<double>());
2306 ParamErrors[kXSecPar].push_back(param["Systematic"]["Error"].as<double>() );
2307 ParamFlat[kXSecPar].push_back(GetFromManager<bool>(param["Systematic"]["FlatPrior"], false));
2308
2309 ParameterGroup.push_back(param["Systematic"]["ParameterGroup"].as<std::string>());
2310
2311 nParam[kXSecPar]++;
2312 ParamType.push_back(kXSecPar);
2313 // Params from osc group have branch name equal to fancy name while all others are basically xsec_0 for example
2314 if(ParameterGroup.back() == "Osc") {
2315 BranchNames.push_back(ParamNames[kXSecPar].back());
2316 } else {
2317 BranchNames.push_back("xsec_" + std::to_string(paramIndex));
2318 }
2319
2320 // Check that the branch exists before setting address
2321 if (!Chain->GetBranch(BranchNames.back())) {
2322 MACH3LOG_ERROR("Couldn't find branch '{}'", BranchNames.back());
2323 throw MaCh3Exception(__FILE__, __LINE__);
2324 }
2325 }
2326}
MCMCProcessor::ExcludedNames
std::vector< std::string > ExcludedNames
Definition: MCMCProcessor.h:396

◆ ReadNDFile()

void MCMCProcessor::ReadNDFile ( )
inlineprotected

Read the ND cov file and get the input central values and errors.

Definition at line 2330 of file MCMCProcessor.cpp.

2330 {
2331// ***************
2332 // Do the same for the ND280
2333 TFile *NDdetFile = new TFile(CovPos[kNDPar].back().c_str(), "open");
2334 if (NDdetFile->IsZombie()) {
2335 MACH3LOG_ERROR("Couldn't find NDdetFile {}", CovPos[kNDPar].back());
2336 throw MaCh3Exception(__FILE__ , __LINE__ );
2337 }
2338 NDdetFile->cd();
2339
2340 TMatrixDSym *NDdetMatrix = NDdetFile->Get<TMatrixDSym>("nddet_cov");
2341 TVectorD *NDdetNominal = NDdetFile->Get<TVectorD>("det_weights");
2342 TDirectory *BinningDirectory = NDdetFile->Get<TDirectory>("Binning");
2343
2344 for (int i = 0; i < NDdetNominal->GetNrows(); ++i)
2345 {
2346 ParamNom[kNDPar].push_back( (*NDdetNominal)(i) );
2347 ParamCentral[kNDPar].push_back( (*NDdetNominal)(i) );
2348
2349 ParamErrors[kNDPar].push_back( std::sqrt((*NDdetMatrix)(i,i)) );
2350 ParamNames[kNDPar].push_back( Form("ND Det %i", i) );
2351 //KS: Currently we can only set it via config, change it in future
2352 ParamFlat[kNDPar].push_back( false );
2353 }
2354
2355 TIter next(BinningDirectory->GetListOfKeys());
2356 TKey *key = nullptr;
2357 // Loop through all entries
2358 while ((key = static_cast<TKey*>(next())))
2359 {
2360 std::string name = std::string(key->GetName());
2361 TH2Poly* RefPoly = BinningDirectory->Get<TH2Poly>((name).c_str());
2362 int size = RefPoly->GetNumberOfBins();
2363 NDSamplesBins.push_back(size);
2364 NDSamplesNames.push_back(RefPoly->GetTitle());
2365 }
2366
2367 NDdetFile->Close();
2368 delete NDdetFile;
2369}

◆ ResetHistograms()

void MCMCProcessor::ResetHistograms ( )

Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.

Definition at line 2458 of file MCMCProcessor.cpp.

2458 {
2459// **************************************************
2460 #ifdef MULTITHREAD
2461 #pragma omp parallel for
2462 #endif
2463 for (int i = 0; i < nDraw; ++i)
2464 {
2465 for (int j = 0; j <= i; ++j)
2466 {
2467 // TH2D to hold the Correlation
2468 hpost2D[i][j]->Reset("");
2469 hpost2D[i][j]->Fill(0.0, 0.0, 0.0);
2470 }
2471 }
2472}

◆ ReweightPrior()

void MCMCProcessor::ReweightPrior ( const std::vector< std::string > &  Names,
const std::vector< double > &  NewCentral,
const std::vector< double > &  NewError 
)

Reweight Prior by giving new central value and new error.

Parameters
NamesParameter names for which we do reweighting
NewCentralNew central value for which we reweight
NewErrorNew error used for calculating weight

Definition at line 2737 of file MCMCProcessor.cpp.

2739 {
2740// **************************
2741 MACH3LOG_INFO("Reweighting Prior");
2742
2743 if( (Names.size() != NewCentral.size()) || (NewCentral.size() != NewError.size()))
2744 {
2745 MACH3LOG_ERROR("Size of passed vectors doesn't match in ReweightPrior");
2746 throw MaCh3Exception(__FILE__ , __LINE__ );
2747 }
2748 std::vector<int> Param;
2749 std::vector<double> OldCentral;
2750 std::vector<double> OldError;
2751 std::vector<bool> FlatPrior;
2752
2753 //KS: First we need to find parameter number based on name
2754 for(unsigned int k = 0; k < Names.size(); ++k)
2755 {
2756 //KS: First we need to find parameter number based on name
2757 int ParamNo = GetParamIndexFromName(Names[k]);
2758 if(ParamNo == _UNDEF_)
2759 {
2760 MACH3LOG_WARN("Couldn't find param {}. Can't reweight Prior", Names[k]);
2761 continue;
2762 }
2763
2764 TString Title = "";
2765 double Prior = 1.0, PriorError = 1.0;
2766 GetNthParameter(ParamNo, Prior, PriorError, Title);
2767
2768 Param.push_back(ParamNo);
2769 OldCentral.push_back(Prior);
2770 OldError.push_back(PriorError);
2771
2772 ParameterEnum ParType = ParamType[ParamNo];
2773 int ParamTemp = ParamNo - ParamTypeStartPos[ParType];
2774
2775 FlatPrior.push_back(ParamFlat[ParType][ParamTemp]);
2776 }
2777 std::vector<double> ParameterPos(Names.size());
2778
2779 std::string InputFile = MCMCFile+".root";
2780 std::string OutputFilename = MCMCFile + "_reweighted.root";
2781
2782 //KS: Simply create copy of file and add there new branch
2783 int ret = system(("cp " + InputFile + " " + OutputFilename).c_str());
2784 if (ret != 0)
2785 MACH3LOG_WARN("Error: system call to copy file failed with code {}", ret);
2786
2787 TFile *OutputChain = new TFile(OutputFilename.c_str(), "UPDATE");
2788 OutputChain->cd();
2789 TTree *post = OutputChain->Get<TTree>("posteriors");
2790
2791 double Weight = 1.;
2792
2793 post->SetBranchStatus("*",false);
2794 // Set the branch addresses for params
2795 for (unsigned int j = 0; j < Names.size(); ++j) {
2796 post->SetBranchStatus(BranchNames[Param[j]].Data(), true);
2797 post->SetBranchAddress(BranchNames[Param[j]].Data(), &ParameterPos[j]);
2798 }
2799 TBranch *bpt = post->Branch("Weight", &Weight, "Weight/D");
2800 post->SetBranchStatus("Weight", true);
2801
2802 for (int i = 0; i < nEntries; ++i)
2803 {
2804 post->GetEntry(i);
2805 Weight = 1.;
2806
2807 //KS: Calculate reweight weight. Weights are multiplicative so we can do several reweights at once. FIXME Big limitation is that code only works for uncorrelated parameters :(
2808 for (unsigned int j = 0; j < Names.size(); ++j)
2809 {
2810 double new_chi = (ParameterPos[j] - NewCentral[j])/NewError[j];
2811 double new_prior = std::exp(-0.5 * new_chi * new_chi);
2812
2813 double old_chi = -1;
2814 double old_prior = -1;
2815 if(FlatPrior[j]) {
2816 old_prior = 1.0;
2817 } else {
2818 old_chi = (ParameterPos[j] - OldCentral[j])/OldError[j];
2819 old_prior = std::exp(-0.5 * old_chi * old_chi);
2820 }
2821 Weight *= new_prior/old_prior;
2822 }
2823 bpt->Fill();
2824 }
2825 post->SetBranchStatus("*",true);
2826 OutputChain->cd();
2827 post->Write("posteriors", TObject::kOverwrite);
2828 OutputChain->Close();
2829 delete OutputChain;
2830
2831 OutputFile->cd();
2832}

◆ ScanInput()

void MCMCProcessor::ScanInput ( )
inlineprotected

Scan Input etc.

Definition at line 1920 of file MCMCProcessor.cpp.

1920 {
1921// **************************
1922 // KS: This can reduce time necessary for caching even by half
1923 #ifdef MULTITHREAD
1924 //ROOT::EnableImplicitMT();
1925 #endif
1926
1927 // Open the Chain
1928 Chain = new TChain("posteriors","posteriors");
1929 Chain->Add(MCMCFile.c_str());
1930
1931 nEntries = int(Chain->GetEntries());
1932
1933 //Only is suboptimality we might want to change it, therefore set it high enough so it doesn't affect other functionality
1934 UpperCut = nEntries+1;
1935
1936 // Get the list of branches
1937 TObjArray* brlis = Chain->GetListOfBranches();
1938
1939 // Get the number of branches
1940 nBranches = brlis->GetEntries();
1941
1942 BranchNames.reserve(nBranches);
1943 ParamType.reserve(nBranches);
1944
1945 // Read the input Covariances
1946 ReadInputCov();
1947
1948 // Set all the branches to off
1949 Chain->SetBranchStatus("*", false);
1950
1951 // Loop over the number of branches
1952 // Find the name and how many of each systematic we have
1953 for (int i = 0; i < nBranches; i++)
1954 {
1955 // Get the TBranch and its name
1956 TBranch* br = static_cast<TBranch*>(brlis->At(i));
1957 if(!br){
1958 MACH3LOG_ERROR("Invalid branch at position {}", i);
1959 throw MaCh3Exception(__FILE__,__LINE__);
1960 }
1961 TString bname = br->GetName();
1962
1963 //KS: Exclude parameter types
1964 bool rejected = false;
1965 for(unsigned int ik = 0; ik < ExcludedTypes.size(); ++ik )
1966 {
1967 if(bname.BeginsWith(ExcludedTypes[ik]))
1968 {
1969 rejected = true;
1970 break;
1971 }
1972 }
1973 if(rejected) continue;
1974
1975 // Turn on the branches which we want for parameters
1976 Chain->SetBranchStatus(bname.Data(), true);
1977
1978 if (bname.BeginsWith("ndd_"))
1979 {
1980 BranchNames.push_back(bname);
1981 ParamType.push_back(kNDPar);
1982 nParam[kNDPar]++;
1983 }
1984 else if (bname.BeginsWith("skd_joint_"))
1985 {
1986 BranchNames.push_back(bname);
1987 ParamType.push_back(kFDDetPar);
1988 nParam[kFDDetPar]++;
1989 }
1990
1991 //KS: as a bonus get LogL systematic
1992 if (bname.BeginsWith("LogL_sample_")) {
1993 SampleName_v.push_back(bname);
1994 nSamples++;
1995 }
1996 else if (bname.BeginsWith("LogL_systematic_")) {
1997 SystName_v.push_back(bname);
1998 nSysts++;
1999 }
2000 }
2001 nDraw = int(BranchNames.size());
2002
2003 // Read the input Covariances
2004 ReadInputCovLegacy();
2005
2006 // Check order of parameter types
2007 ScanParameterOrder();
2008
2009 IamVaried.resize(nDraw, true);
2010
2011 // Print useful Info
2012 PrintInfo();
2013
2014 nSteps = Chain->GetMaximum("step");
2015 // Set the step cut to be 20%
2016 int cut = nSteps/5;
2017 SetStepCut(cut);
2018
2019 // Basically allow loading oscillation parameters
2020 LoadAdditionalInfo();
2021}
MCMCProcessor::ScanParameterOrder
void ScanParameterOrder()
Scan order of params from a different groups.
Definition: MCMCProcessor.cpp:2082
MCMCProcessor::PrintInfo
void PrintInfo() const
Print info like how many params have been loaded etc.
Definition: MCMCProcessor.cpp:4131
MCMCProcessor::LoadAdditionalInfo
virtual void LoadAdditionalInfo()
allow loading additional info for example used for oscillation parameters
Definition: MCMCProcessor.h:293
MCMCProcessor::ExcludedTypes
std::vector< std::string > ExcludedTypes
Definition: MCMCProcessor.h:395
MCMCProcessor::ReadInputCovLegacy
void ReadInputCovLegacy()
Definition: MCMCProcessor.cpp:2171
MCMCProcessor::SetStepCut
void SetStepCut(const std::string &Cuts)
Set the step cutting by string.
Definition: MCMCProcessor.cpp:2407
MCMCProcessor::nBranches
int nBranches
Number of branches in a TTree.
Definition: MCMCProcessor.h:381
MCMCProcessor::ReadInputCov
void ReadInputCov()
CW: Read the input Covariance matrix entries. Get stuff like parameter input errors,...
Definition: MCMCProcessor.cpp:2162

◆ ScanParameterOrder()

void MCMCProcessor::ScanParameterOrder ( )
inlineprotected

Scan order of params from a different groups.

Definition at line 2082 of file MCMCProcessor.cpp.

2082 {
2083// *****************************
2084 for(int i = 0; i < kNParameterEnum; i++)
2085 {
2086 for(unsigned int j = 0; j < ParamType.size(); j++)
2087 {
2088 if(ParamType[j] == ParameterEnum(i))
2089 {
2090 //KS: When we find that i-th parameter types start at j, save and move to the next parameter.
2091 ParamTypeStartPos[i] = j;
2092 break;
2093 }
2094 }
2095 }
2096}

◆ SetExcludedNames()

void MCMCProcessor::SetExcludedNames ( std::vector< std::string >  Name)
inline

Definition at line 263 of file MCMCProcessor.h.

263{ExcludedNames = Name; };

◆ SetExcludedTypes()

void MCMCProcessor::SetExcludedTypes ( std::vector< std::string >  Name)
inline

Setter related what parameters we want to exclude from analysis, for example if cross-section parameters look like xsec_, then passing "xsec_" will.

Parameters
BatchesVector with parameters type names we want to exclude

Definition at line 262 of file MCMCProcessor.h.

262{ExcludedTypes = Name; };

◆ SetFancyNames()

void MCMCProcessor::SetFancyNames ( const bool  PlotOrNot)
inline

Definition at line 251 of file MCMCProcessor.h.

251{FancyPlotNames = PlotOrNot; };

◆ SetLegendStyle()

void MCMCProcessor::SetLegendStyle ( TLegend *  Legend,
const double  size 
) const
protected

Configures the style of a TLegend object.

Parameters
LegendPointer to the TLegend object to modify
sizeThe text size to set for the legend

Definition at line 4190 of file MCMCProcessor.cpp.

4190 {
4191// **************************
4192 Legend->SetTextSize(size);
4193 Legend->SetLineColor(0);
4194 Legend->SetLineStyle(0);
4195 Legend->SetFillColor(0);
4196 Legend->SetFillStyle(0);
4197 Legend->SetBorderSize(0);
4198}

◆ SetMargins()

void MCMCProcessor::SetMargins ( std::unique_ptr< TCanvas > &  Canv,
const std::vector< double > &  margins 
)
protected

Set TCanvas margins to specified values.

Definition at line 4165 of file MCMCProcessor.cpp.

4165 {
4166// **************************
4167 if (!Canv) {
4168 MACH3LOG_ERROR("Canv is nullptr");
4169 throw MaCh3Exception(__FILE__, __LINE__);
4170 }
4171 if (margins.size() != 4) {
4172 MACH3LOG_ERROR("Margin vector must have exactly 4 elements");
4173 throw MaCh3Exception(__FILE__, __LINE__);
4174 }
4175 Canv->SetTopMargin(margins[0]);
4176 Canv->SetBottomMargin(margins[1]);
4177 Canv->SetLeftMargin(margins[2]);
4178 Canv->SetRightMargin(margins[3]);
4179}

◆ SetnBatches()

void MCMCProcessor::SetnBatches ( const int  Batches)
inline

Set value of Nbatches used for batched mean, this need to be done earlier as batches are made when reading tree.

Parameters
BatchesNumber of batches, default is 20

Definition at line 267 of file MCMCProcessor.h.

267{nBatches = Batches; };

◆ SetnLags()

void MCMCProcessor::SetnLags ( const int  nLags)
inline

Definition at line 268 of file MCMCProcessor.h.

268{AutoCorrLag = nLags; };

◆ SetOutputSuffix()

void MCMCProcessor::SetOutputSuffix ( const std::string  Suffix)
inline

Sett output suffix, this way jobs using the same file will have different names.

Definition at line 271 of file MCMCProcessor.h.

271{OutputSuffix = Suffix; };

◆ SetPlotBinValue()

void MCMCProcessor::SetPlotBinValue ( const bool  PlotOrNot)
inline

Definition at line 250 of file MCMCProcessor.h.

250{plotBinValue = PlotOrNot; };

◆ SetPlotErrorForFlatPrior()

void MCMCProcessor::SetPlotErrorForFlatPrior ( const bool  PlotOrNot)
inline

Set whether you want to plot error for parameters which have flat prior.

Definition at line 249 of file MCMCProcessor.h.

249{PlotFlatPrior = PlotOrNot; };

◆ SetPlotRelativeToPrior()

void MCMCProcessor::SetPlotRelativeToPrior ( const bool  PlotOrNot)
inline

You can set relative to prior or relative to generated. It is advised to use relate to prior.

Parameters
PlotOrNotbool controlling plotRelativeToPrior argument

Definition at line 246 of file MCMCProcessor.h.

246{plotRelativeToPrior = PlotOrNot; };

◆ SetPost2DPlotThreshold()

void MCMCProcessor::SetPost2DPlotThreshold ( const double  Threshold)
inline

Code will only plot 2D posteriors if Correlation are larger than defined threshold.

Parameters
ThresholdThis threshold is compared with correlation value

Definition at line 256 of file MCMCProcessor.h.

256{Post2DPlotThreshold = Threshold; };

◆ SetPosterior1DCut()

void MCMCProcessor::SetPosterior1DCut ( const std::string  Cut)
inline

Allow to set addtional cuts based on ROOT TBrowser cut, for to only affect one mass ordering.

Definition at line 273 of file MCMCProcessor.h.

273{Posterior1DCut = Cut; };

◆ SetPrintToPDF()

void MCMCProcessor::SetPrintToPDF ( const bool  PlotOrNot)
inline

Definition at line 247 of file MCMCProcessor.h.

247{printToPDF = PlotOrNot; };

◆ SetSmoothing()

void MCMCProcessor::SetSmoothing ( const bool  PlotOrNot)
inline

Set whether want to use smoothing for histograms using ROOT algorithm.

Definition at line 253 of file MCMCProcessor.h.

253{ApplySmoothing = PlotOrNot; };

◆ SetStepCut() [1/2]

void MCMCProcessor::SetStepCut ( const int  Cuts)

Set the step cutting by int.

Parameters
Cutsinteger telling cut value

Definition at line 2415 of file MCMCProcessor.cpp.

2415 {
2416// ***************
2417 std::stringstream TempStream;
2418 TempStream << "step > " << Cuts;
2419 StepCut = TempStream.str();
2420 BurnInCut = Cuts;
2421}

◆ SetStepCut() [2/2]

void MCMCProcessor::SetStepCut ( const std::string &  Cuts)

Set the step cutting by string.

Parameters
Cutsstring telling cut value

Definition at line 2407 of file MCMCProcessor.cpp.

2407 {
2408// ***************
2409 StepCut = Cuts;
2410 BurnInCut = std::stoi( Cuts );
2411}

◆ SetTLineStyle()

void MCMCProcessor::SetTLineStyle ( TLine *  Line,
const Color_t  Colour,
const Width_t  Width,
const ELineStyle  Style 
) const
protected

Configures a TLine object with the specified style parameters.

Parameters
LinePointer to the TLine object to modify. Must not be nullptr.
ColourThe color to set for the line.
WidthThe width of the line.
StyleThe line style (e.g., solid, dashed, etc.).

Definition at line 4182 of file MCMCProcessor.cpp.

4182 {
4183// **************************
4184 Line->SetLineColor(Colour);
4185 Line->SetLineWidth(Width);
4186 Line->SetLineStyle(Style);
4187}

◆ SetupOutput()

void MCMCProcessor::SetupOutput ( )
inlineprotected

Prepare all objects used for output.

Definition at line 2025 of file MCMCProcessor.cpp.

2025 {
2026// ****************************
2027 // Make sure we can read files located anywhere and strip the .root ending
2028 MCMCFile = MCMCFile.substr(0, MCMCFile.find(".root"));
2029
2030 // Check if the output file is ready
2031 if (OutputFile == nullptr) MakeOutputFile();
2032
2033 CanvasName = MCMCFile + OutputSuffix + ".pdf[";
2034 if(printToPDF) Posterior->Print(CanvasName);
2035
2036 // Once the pdf file is open no longer need to bracket
2037 CanvasName.ReplaceAll("[","");
2038
2039 // We fit with this Gaussian
2040 Gauss = std::make_unique<TF1>("Gauss", "[0]/sqrt(2.0*3.14159)/[2]*TMath::Exp(-0.5*pow(x-[1],2)/[2]/[2])", -5, 5);
2041 Gauss->SetLineWidth(2);
2042 Gauss->SetLineColor(kOrange-5);
2043
2044 // Declare the TVectors
2045 Covariance = new TMatrixDSym(nDraw);
2046 Correlation = new TMatrixDSym(nDraw);
2047 Central_Value = new TVectorD(nDraw);
2048 Means = new TVectorD(nDraw);
2049 Errors = new TVectorD(nDraw);
2050 Means_Gauss = new TVectorD(nDraw);
2051 Errors_Gauss = new TVectorD(nDraw);
2052 Means_HPD = new TVectorD(nDraw);
2053 Errors_HPD = new TVectorD(nDraw);
2054 Errors_HPD_Positive = new TVectorD(nDraw);
2055 Errors_HPD_Negative = new TVectorD(nDraw);
2056
2057 // Initialise to something silly
2058 #ifdef MULTITHREAD
2059 #pragma omp parallel for
2060 #endif
2061 for (int i = 0; i < nDraw; ++i)
2062 {
2063 (*Central_Value)(i) = _UNDEF_;
2064 (*Means)(i) = _UNDEF_;
2065 (*Errors)(i) = _UNDEF_;
2066 (*Means_Gauss)(i) = _UNDEF_;
2067 (*Errors_Gauss)(i) = _UNDEF_;
2068 (*Means_HPD)(i) = _UNDEF_;
2069 (*Errors_HPD)(i) = _UNDEF_;
2070 (*Errors_HPD_Positive)(i) = _UNDEF_;
2071 (*Errors_HPD_Negative)(i) = _UNDEF_;
2072 for (int j = 0; j < nDraw; ++j) {
2073 (*Covariance)(i, j) = _UNDEF_;
2074 (*Correlation)(i, j) = _UNDEF_;
2075 }
2076 }
2077 hpost.resize(nDraw);
2078}

◆ SetUseFFTAutoCorrelation()

void MCMCProcessor::SetUseFFTAutoCorrelation ( const bool  useFFT)
inline

Toggle using the FFT-based autocorrelation calculator.

Definition at line 258 of file MCMCProcessor.h.

258{useFFTAutoCorrelation = useFFT; };

◆ ThinMCMC()

void MCMCProcessor::ThinMCMC ( const int  ThinningCut)
inline

Thin MCMC Chain, to save space and maintain low autocorrelations.

Parameters
ThinningCutevery which entry you want to thin
AverageIf true will perform MCMC averaging instead of thinning

Definition at line 186 of file MCMCProcessor.h.

186{ ThinningMCMC(MCMCFile+".root", ThinningCut); };
ThinningMCMC
void ThinningMCMC(const std::string &FilePath, const int ThinningCut)
Thin MCMC Chain, to save space and maintain low autocorrelations.
Definition: StatisticalUtils.cpp:487

Member Data Documentation

◆ AccProbBatchedAverages

double* MCMCProcessor::AccProbBatchedAverages
protected

Holds all accProb in batches.

Definition at line 530 of file MCMCProcessor.h.

◆ AccProbValues

double* MCMCProcessor::AccProbValues
protected

Holds all accProb.

Definition at line 528 of file MCMCProcessor.h.

◆ ApplySmoothing

bool MCMCProcessor::ApplySmoothing
protected

Apply smoothing for 2D histos using root algorithm.

Definition at line 444 of file MCMCProcessor.h.

◆ AutoCorrLag

int MCMCProcessor::AutoCorrLag
protected

LagL used in AutoCorrelation.

Definition at line 515 of file MCMCProcessor.h.

◆ BatchedAverages

double** MCMCProcessor::BatchedAverages
protected

Values of batched average for every param and batch.

Definition at line 520 of file MCMCProcessor.h.

◆ BranchNames

std::vector<TString> MCMCProcessor::BranchNames
protected

Definition at line 394 of file MCMCProcessor.h.

◆ BurnInCut

int MCMCProcessor::BurnInCut
protected

Value of burn in cut.

Definition at line 379 of file MCMCProcessor.h.

◆ CacheMCMC

bool MCMCProcessor::CacheMCMC
protected

MCMC Chain has been cached.

Definition at line 507 of file MCMCProcessor.h.

◆ CanvasName

TString MCMCProcessor::CanvasName
protected

Name of canvas which help to save to the sample pdf.

Definition at line 427 of file MCMCProcessor.h.

◆ Central_Value

TVectorD* MCMCProcessor::Central_Value
protected

Vector with central value for each parameter.

Definition at line 464 of file MCMCProcessor.h.

◆ Chain

TChain* MCMCProcessor::Chain
protected

Main chain storing all steps etc.

Definition at line 371 of file MCMCProcessor.h.

◆ Correlation

TMatrixDSym* MCMCProcessor::Correlation
protected

Posterior Correlation Matrix.

Definition at line 485 of file MCMCProcessor.h.

◆ Covariance

TMatrixDSym* MCMCProcessor::Covariance
protected

Posterior Covariance Matrix.

Definition at line 483 of file MCMCProcessor.h.

◆ CovConfig

std::vector<YAML::Node> MCMCProcessor::CovConfig
protected

Covariance matrix config.

Definition at line 368 of file MCMCProcessor.h.

◆ CovPos

std::vector<std::vector<std::string> > MCMCProcessor::CovPos
protected

Covariance matrix name position.

Definition at line 366 of file MCMCProcessor.h.

◆ doDiagMCMC

bool MCMCProcessor::doDiagMCMC
protected

Doing MCMC Diagnostic.

Definition at line 509 of file MCMCProcessor.h.

◆ DrawRange

double MCMCProcessor::DrawRange
protected

Drawrange for SetMaximum.

Definition at line 504 of file MCMCProcessor.h.

◆ Errors

TVectorD* MCMCProcessor::Errors
protected

Vector with errors values using RMS.

Definition at line 468 of file MCMCProcessor.h.

◆ Errors_Gauss

TVectorD* MCMCProcessor::Errors_Gauss
protected

Vector with error values using Gaussian fit.

Definition at line 472 of file MCMCProcessor.h.

◆ Errors_HPD

TVectorD* MCMCProcessor::Errors_HPD
protected

Vector with error values using Highest Posterior Density.

Definition at line 476 of file MCMCProcessor.h.

◆ Errors_HPD_Negative

TVectorD* MCMCProcessor::Errors_HPD_Negative
protected

Vector with negative error (left hand side) values using Highest Posterior Density.

Definition at line 480 of file MCMCProcessor.h.

◆ Errors_HPD_Positive

TVectorD* MCMCProcessor::Errors_HPD_Positive
protected

Vector with positive error (right hand side) values using Highest Posterior Density.

Definition at line 478 of file MCMCProcessor.h.

◆ ExcludedNames

std::vector<std::string> MCMCProcessor::ExcludedNames
protected

Definition at line 396 of file MCMCProcessor.h.

◆ ExcludedTypes

std::vector<std::string> MCMCProcessor::ExcludedTypes
protected

Definition at line 395 of file MCMCProcessor.h.

◆ FancyPlotNames

bool MCMCProcessor::FancyPlotNames
protected

Whether we want fancy plot names or not.

Definition at line 440 of file MCMCProcessor.h.

◆ Gauss

std::unique_ptr<TF1> MCMCProcessor::Gauss
protected

Gaussian fitter.

Definition at line 454 of file MCMCProcessor.h.

◆ hpost

std::vector<TH1D*> MCMCProcessor::hpost
protected

Holds 1D Posterior Distributions.

Definition at line 488 of file MCMCProcessor.h.

◆ hpost2D

std::vector<std::vector<TH2D*> > MCMCProcessor::hpost2D
protected

Holds 2D Posterior Distributions.

Definition at line 490 of file MCMCProcessor.h.

◆ hviolin

std::unique_ptr<TH2D> MCMCProcessor::hviolin
protected

Holds violin plot for all dials.

Definition at line 492 of file MCMCProcessor.h.

◆ hviolin_prior

std::unique_ptr<TH2D> MCMCProcessor::hviolin_prior
protected

Holds prior violin plot for all dials,.

Definition at line 494 of file MCMCProcessor.h.

◆ IamVaried

std::vector<bool> MCMCProcessor::IamVaried
protected

Is the ith parameter varied.

Definition at line 399 of file MCMCProcessor.h.

◆ MadePostfit

bool MCMCProcessor::MadePostfit
protected

Sanity check if Postfit is already done to not make several times.

Definition at line 436 of file MCMCProcessor.h.

◆ MCMCFile

std::string MCMCProcessor::MCMCFile
protected

Name of MCMC file.

Definition at line 362 of file MCMCProcessor.h.

◆ Means

TVectorD* MCMCProcessor::Means
protected

Vector with mean values using Arithmetic Mean.

Definition at line 466 of file MCMCProcessor.h.

◆ Means_Gauss

TVectorD* MCMCProcessor::Means_Gauss
protected

Vector with mean values using Gaussian fit.

Definition at line 470 of file MCMCProcessor.h.

◆ Means_HPD

TVectorD* MCMCProcessor::Means_HPD
protected

Vector with mean values using Highest Posterior Density.

Definition at line 474 of file MCMCProcessor.h.

◆ nBatches

int MCMCProcessor::nBatches
protected

Number of batches for Batched Mean.

Definition at line 513 of file MCMCProcessor.h.

◆ nBins

int MCMCProcessor::nBins
protected

Number of bins.

Definition at line 502 of file MCMCProcessor.h.

◆ nBranches

int MCMCProcessor::nBranches
protected

Number of branches in a TTree.

Definition at line 381 of file MCMCProcessor.h.

◆ nDraw

int MCMCProcessor::nDraw
protected

Number of all parameters used in the analysis.

Definition at line 391 of file MCMCProcessor.h.

◆ NDSamplesBins

std::vector<int> MCMCProcessor::NDSamplesBins
protected

Definition at line 450 of file MCMCProcessor.h.

◆ NDSamplesNames

std::vector<std::string> MCMCProcessor::NDSamplesNames
protected

Definition at line 451 of file MCMCProcessor.h.

◆ nEntries

int MCMCProcessor::nEntries
protected

KS: For merged chains number of entries will be different from nSteps.

Definition at line 383 of file MCMCProcessor.h.

◆ nParam

std::vector<int> MCMCProcessor::nParam
protected

Number of parameters per type.

Definition at line 410 of file MCMCProcessor.h.

◆ nSamples

int MCMCProcessor::nSamples
protected

Number of sample PDF objects.

Definition at line 387 of file MCMCProcessor.h.

◆ nSteps

int MCMCProcessor::nSteps
protected

KS: For merged chains number of entries will be different from nSteps.

Definition at line 385 of file MCMCProcessor.h.

◆ nSysts

int MCMCProcessor::nSysts
protected

Number of covariance objects.

Definition at line 389 of file MCMCProcessor.h.

◆ OutputFile

TFile* MCMCProcessor::OutputFile
protected

The output file.

Definition at line 457 of file MCMCProcessor.h.

◆ OutputName

std::string MCMCProcessor::OutputName
protected

Name of output files.

Definition at line 425 of file MCMCProcessor.h.

◆ OutputSuffix

std::string MCMCProcessor::OutputSuffix
protected

Output file suffix useful when running over same file with different settings.

Definition at line 364 of file MCMCProcessor.h.

◆ ParamCentral

std::vector<std::vector<double> > MCMCProcessor::ParamCentral
protected

Parameters central values which we are going to analyse.

Definition at line 403 of file MCMCProcessor.h.

◆ ParamErrors

std::vector<std::vector<double> > MCMCProcessor::ParamErrors
protected

Uncertainty on a single parameter.

Definition at line 406 of file MCMCProcessor.h.

◆ ParameterGroup

std::vector<std::string> MCMCProcessor::ParameterGroup
protected

Definition at line 417 of file MCMCProcessor.h.

◆ ParamFlat

std::vector<std::vector<bool> > MCMCProcessor::ParamFlat
protected

Whether Param has flat prior or not.

Definition at line 408 of file MCMCProcessor.h.

◆ ParamNames

std::vector<std::vector<TString> > MCMCProcessor::ParamNames
protected

Name of parameters which we are going to analyse.

Definition at line 401 of file MCMCProcessor.h.

◆ ParamNom

std::vector<std::vector<double> > MCMCProcessor::ParamNom
protected

Definition at line 404 of file MCMCProcessor.h.

◆ ParamSums

double* MCMCProcessor::ParamSums
protected

Total parameter sum for each param.

Definition at line 518 of file MCMCProcessor.h.

◆ ParamType

std::vector<ParameterEnum> MCMCProcessor::ParamType
protected

Make an enum for which class this parameter belongs to so we don't have to keep string comparing.

Definition at line 412 of file MCMCProcessor.h.

◆ ParamTypeStartPos

std::vector<int> MCMCProcessor::ParamTypeStartPos
protected

KS: in MCMC output there is order of parameters so for example first goes xsec then nd det etc. Idea is that this parameter will keep track of it so code is flexible

Definition at line 415 of file MCMCProcessor.h.

◆ ParStep

double** MCMCProcessor::ParStep
protected

Array holding values for all parameters.

Definition at line 497 of file MCMCProcessor.h.

◆ plotBinValue

bool MCMCProcessor::plotBinValue
protected

If true it will print value on each bin of covariance matrix.

Definition at line 442 of file MCMCProcessor.h.

◆ PlotFlatPrior

bool MCMCProcessor::PlotFlatPrior
protected

Whether we plot flat prior or not, we usually provide error even for flat prior params.

Definition at line 430 of file MCMCProcessor.h.

◆ plotRelativeToPrior

bool MCMCProcessor::plotRelativeToPrior
protected

Whether we plot relative to prior or nominal, in most cases is prior.

Definition at line 434 of file MCMCProcessor.h.

◆ Post2DPlotThreshold

double MCMCProcessor::Post2DPlotThreshold
protected

KS: Set Threshold when to plot 2D posterior as by default we get a LOT of plots.

Definition at line 446 of file MCMCProcessor.h.

◆ Posterior

std::unique_ptr<TCanvas> MCMCProcessor::Posterior
protected

Fancy canvas used for our beautiful plots.

Definition at line 460 of file MCMCProcessor.h.

◆ Posterior1DCut

std::string MCMCProcessor::Posterior1DCut
protected

Cut used when making 1D Posterior distribution.

Definition at line 375 of file MCMCProcessor.h.

◆ printToPDF

bool MCMCProcessor::printToPDF
protected

Will plot all plot to PDF not only to root file.

Definition at line 438 of file MCMCProcessor.h.

◆ SampleName_v

std::vector<TString> MCMCProcessor::SampleName_v
protected

Vector of each systematic.

Definition at line 420 of file MCMCProcessor.h.

◆ SampleValues

double** MCMCProcessor::SampleValues
protected

Holds the sample values.

Definition at line 523 of file MCMCProcessor.h.

◆ StepCut

std::string MCMCProcessor::StepCut
protected

BurnIn Cuts.

Definition at line 373 of file MCMCProcessor.h.

◆ StepNumber

int* MCMCProcessor::StepNumber
protected

Step number for step, important if chains were merged.

Definition at line 499 of file MCMCProcessor.h.

◆ SystName_v

std::vector<TString> MCMCProcessor::SystName_v
protected

Vector of each sample PDF object.

Definition at line 422 of file MCMCProcessor.h.

◆ SystValues

double** MCMCProcessor::SystValues
protected

Holds the systs values.

Definition at line 525 of file MCMCProcessor.h.

◆ UpperCut

int MCMCProcessor::UpperCut
protected

KS: Used only for SubOptimality.

Definition at line 377 of file MCMCProcessor.h.

◆ useFFTAutoCorrelation

bool MCMCProcessor::useFFTAutoCorrelation
protected

MJR: Use FFT-based autocorrelation algorithm (save time & resources)?

Definition at line 448 of file MCMCProcessor.h.

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.6