This class extends MCMC and allow specialised for Oscillation parameters analysis which require specialised hardcoding. More...

#include <Fitters/OscProcessor.h>

Inheritance diagram for OscProcessor:

[legend]

Collaboration diagram for OscProcessor:

[legend]

Public Member Functions
	OscProcessor (const std::string &InputFile)
	Constructs an OscProcessor object with the specified input file and options. More...

virtual	~OscProcessor ()
	Destroys the OscProcessor object. More...

void	PerformJarlskogAnalysis ()
	Perform Several Jarlskog Plotting. More...

void	MakePiePlot ()
	Make fancy Pie plot for delta CP. More...

Public Member Functions inherited from MCMCProcessor
	MCMCProcessor (const std::string &InputFile)
	Constructs an MCMCProcessor object with the specified input file and options. More...

virtual	~MCMCProcessor ()
	Destroys the MCMCProcessor object. More...

void	Initialise ()
	Scan chain, what parameters we have and load information from covariance matrices. More...

void	MakePostfit (const std::map< std::string, std::pair< double, double >> &Edges={})
	Make 1D projection for each parameter and prepare structure. More...

void	MakeCovariance ()
	Calculate covariance by making 2D projection of each combination of parameters. More...

void	CacheSteps ()
	KS:By caching each step we use multithreading. More...

void	MakeCovariance_MP (const bool Mute=false)
	Calculate covariance by making 2D projection of each combination of parameters using multithreading. More...

void	MakeSubOptimality (const int NIntervals=10)
	Make and Draw SubOptimality [28]. More...

void	Reset2DPosteriors ()
	Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut. More...

void	DrawPostfit ()
	Draw the post-fit comparisons. More...

void	MakeViolin ()
	Make and Draw Violin. More...

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. More...

void	DrawCovariance ()
	Draw the post-fit covariances. More...

void	MakeCovarianceYAML (const std::string &OutputYAMLFile, const std::string &MeansMethod) const
	Make YAML file from post-fit covariance. More...

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. More...

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. More...

void	CheckCredibleIntervalsOrder (const std::vector< double > &CredibleIntervals, const std::vector< Color_t > &CredibleIntervalsColours) const
	Checks the order and size consistency of the `CredibleIntervals` and `CredibleIntervalsColours` vectors. More...

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. More...

void	GetPolarPlot (const std::vector< std::string > &ParNames)
	Make funny polar plot. More...

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. More...

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. More...

void	SavageDickeyPlot (std::unique_ptr< TH1D > &PriorHist, std::unique_ptr< TH1D > &PosteriorHist, const std::string &Title, const double EvaluationPoint) const
	Produce Savage Dickey plot. More...

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. More...

void	SmearChain (const std::vector< std::string > &Names, const std::vector< double > &Error, const bool &SaveBranch) const
	Smear chain contours. More...

void	ParameterEvolution (const std::vector< std::string > &Names, const std::vector< int > &NIntervals)
	Make .gif of parameter evolution. More...

void	ThinMCMC (const int ThinningCut) const
	Thin MCMC Chain, to save space and maintain low autocorrelations. More...

void	DiagMCMC ()
	KS: Perform MCMC diagnostic including Autocorrelation, Trace etc. More...

int	GetNParams () const
	Get total number of used parameters. More...

int	GetNXSec () const

int	GetNND () const

int	GetNFD () const

YAML::Node	GetCovConfig (const int i) const
	Get Yaml config obtained from a Chain. More...

int	GetGroup (const std::string &name) const
	Number of params from a given group, for example flux. More...

TH1D *	GetHpost (const int i) const
	Get 1D posterior for a given parameter. More...

TH2D *	GetHpost2D (const int i, const int j) const
	Get 2D posterior for a given parameter combination. More...

TH2D *	GetViolin () const
	Get Violin plot for all parameters with posterior values. More...

TH2D *	GetViolinPrior () const
	Get Violin plot for all parameters with prior values. More...

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) More...

void	GetCovariance (TMatrixDSym &Cov, TMatrixDSym &Corr)
	Get the post-fit covariances and correlations. More...

void	GetPostfit_Ind (TVectorD &Central, TVectorD &Errors, TVectorD *&Peaks, ParameterEnum kParam)
	Or the individual post-fits. More...

const std::vector< TString > &	GetBranchNames () const
	Get the vector of branch names from root file. More...

void	GetNthParameter (const int param, double &Prior, double &PriorError, TString &Title) const
	Get properties of parameter by passing it number. More...

int	GetParamIndexFromName (const std::string &Name) const
	Get parameter number based on name. More...

Long64_t	GetnEntries ()
	Get Number of entries that Chain has, for merged chains will not be the same Nsteps. More...

Long64_t	GetnSteps ()
	Get Number of Steps that Chain has, for merged chains will not be the same nEntries. More...

void	SetNBins (const int NewBins)
	Modify number of bins used for 1D and 2D Histograms. More...

void	SetEntries (const int NewEntries)
	Set number of entries to make potentially MCMC Processing faster. More...

void	SetStepCut (const std::string &Cuts)
	Set the step cutting by string. More...

void	SetStepCut (const int Cuts)
	Set the step cutting by int. More...

void	CheckStepCut () const
	Check if step cut isn't larger than highest values of step in a chain. More...

void	SetPlotRelativeToPrior (const bool PlotOrNot)
	You can set relative to prior or relative to generated. It is advised to use relate to prior. More...

void	SetPrintToPDF (const bool PlotOrNot)
	Whether to dump all plots into PDF. More...

void	SetPlotErrorForFlatPrior (const bool PlotOrNot)
	Set whether you want to plot error for parameters which have flat prior. More...

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. More...

void	SetPost2DPlotThreshold (const double Threshold)
	Code will only plot 2D posteriors if Correlation are larger than defined threshold. More...

void	SetUseFFTAutoCorrelation (const bool useFFT)
	Toggle using the FFT-based autocorrelation calculator. More...

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 param_, then passing "param_" will. More...

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

void	SetExcludedGroups (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. More...

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. More...

void	SetPosterior1DCut (const std::string Cut)
	Allow to set addtional cuts based on ROOT TBrowser cut, for to only affect one mass ordering. More...

Protected Member Functions
void	LoadAdditionalInfo () override
	Read the Osc cov file and get the input central values and errors Here we allow Jarlskog Shenanigans. More...

void	MakeJarlskogPlot (const std::unique_ptr< TH1D > &jarl, const std::unique_ptr< TH1D > &jarl_flatsindcp, const std::unique_ptr< TH1D > &jarl_NH, const std::unique_ptr< TH1D > &jarl_NH_flatsindcp, const std::unique_ptr< TH1D > &jarl_IH, const std::unique_ptr< TH1D > &jarl_IH_flatsindcp)
	Perform Jarlskog Plotting. More...

double	CalcJarlskog (const double s2th13, const double s2th23, const double s2th12, const double dcp) const
	Calculate Jarlskog Invariant using oscillation parameters. More...

double	SamplePriorForParam (const int paramIndex, const std::unique_ptr< TRandom3 > &randGen, const std::vector< double > &FlatBounds) const
	Draw Prior value. More...

Protected Member Functions inherited from MCMCProcessor
std::unique_ptr< TH1D >	MakePrefit ()
	Prepare prefit histogram for parameter overlay plot. More...

void	MakeOutputFile ()
	prepare output root file and canvas to which we will save EVERYTHING More...

void	DrawCorrelations1D ()
	Draw 1D correlations which might be more helpful than looking at huge 2D Corr matrix. More...

void	DrawCorrelationsGroup (const std::unique_ptr< TH2D > &CorrMatrix) const
	Produces correlation matrix but instead of giving name for each param it only give name for param group. More...

void	ReadInputCov ()
	CW: Read the input Covariance matrix entries. Get stuff like parameter input errors, names, and so on. More...

void	ReadInputCovLegacy ()

void	FindInputFiles ()
	Read the output MCMC file and find what inputs were used. More...

void	FindInputFilesLegacy ()

void	ReadModelFile ()
	Read the xsec file and get the input central values and errors. More...

void	ReadNDFile ()
	Read the ND cov file and get the input central values and errors. More...

void	ReadFDFile ()
	Read the FD cov file and get the input central values and errors. More...

void	PrintInfo () const
	Print info like how many params have been loaded etc. More...

void	ScanInput ()
	Scan Input etc. More...

void	ScanParameterOrder ()
	Scan order of params from a different groups. More...

void	SetupOutput ()
	Prepare all objects used for output. More...

void	PrepareDiagMCMC ()
	CW: Prepare branches etc. for DiagMCMC. More...

std::vector< double >	GetParameterSums ()
	Computes the average of each parameter across all MCMC entries. Useful for autocorrelation. More...

void	ParamTraces ()
	CW: Draw trace plots of the parameters i.e. parameter vs step. More...

void	AutoCorrelation ()
	KS: Calculate autocorrelations supports both OpenMP and CUDA :) More...

void	AutoCorrelation_FFT ()
	MJR: Autocorrelation function using FFT algorithm for extra speed. More...

void	CalculateESS (const int nLags, const std::vector< std::vector< double >> &LagL)
	KS: calc Effective Sample Size. More...

void	BatchedAnalysis ()
	Get the batched means variance estimation and variable indicating if number of batches is sensible [4] [29]. More...

void	BatchedMeans ()
	CW: Batched means, literally read from an array and chuck into TH1D. More...

void	GewekeDiagnostic ()
	Geweke Diagnostic based on the methods described by Fang (2014) and Karlsbakk (2011). [8] [21]. More...

void	AcceptanceProbabilities ()
	Acceptance Probability. More...

void	PowerSpectrumAnalysis ()
	RC: Perform spectral analysis of MCMC [7]. More...

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. More...

void	SetMargins (std::unique_ptr< TCanvas > &Canv, const std::vector< double > &margins)
	Set TCanvas margins to specified values. More...

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. More...

void	SetLegendStyle (TLegend *Legend, const double size) const
	Configures the style of a TLegend object. More...

Protected Attributes
bool	PlotJarlskog
	Will plot Jarlskog Invariant using information in the chain. More...

bool	OscEnabled
	Will plot Jarlskog Invariant using information in the chain. More...

std::string	Sin2Theta13Name
	Name of the parameter representing \(\sin^2\theta_{13}\). More...

std::string	Sin2Theta12Name
	Name of the parameter representing \(\sin^2\theta_{12}\). More...

std::string	Sin2Theta23Name
	Name of the parameter representing \(\sin^2\theta_{23}\). More...

std::string	DeltaCPName
	Name of the parameter representing \(\delta_{\mathrm{CP}}\) (the CP-violating phase). More...

std::string	DeltaM2_23Name
	Name of the parameter representing \(\Delta m^2_{32}\) (mass-squared difference). More...

int	Sin2Theta13Index
	Index of \(\sin^2\theta_{13}\) in the parameter list. More...

int	Sin2Theta12Index
	Index of \(\sin^2\theta_{12}\) in the parameter list. More...

int	Sin2Theta23Index
	Index of \(\sin^2\theta_{23}\) in the parameter list. More...

int	DeltaCPIndex
	Index of \(\delta_{\mathrm{CP}}\) in the parameter list. More...

int	DeltaM2_23Index
	Index of \(\Delta m^2_{32}\) in the parameter list. More...

Protected Attributes inherited from MCMCProcessor
std::string	MCMCFile
	Name of MCMC file. More...

std::string	OutputSuffix
	Output file suffix useful when running over same file with different settings. More...

std::vector< std::vector< std::string > >	CovPos
	Covariance matrix file name position. More...

std::vector< std::string >	CovNamePos
	Covariance matrix name position. More...

std::vector< YAML::Node >	CovConfig
	Covariance matrix config. More...

TChain *	Chain
	Main chain storing all steps etc. More...

std::string	StepCut
	BurnIn Cuts. More...

std::string	Posterior1DCut
	Cut used when making 1D Posterior distribution. More...

unsigned int	UpperCut
	KS: Used only for SubOptimality. More...

unsigned int	BurnInCut
	Value of burn in cut. More...

int	nBranches
	Number of branches in a TTree. More...

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

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

int	nSampleHandlers
	Number of sample PDF objects. More...

int	nParameterHandlers
	Number of covariance objects. More...

int	nDraw
	Number of all parameters used in the analysis. More...

std::vector< TString >	BranchNames

std::vector< std::string >	ExcludedTypes

std::vector< std::string >	ExcludedNames

std::vector< std::string >	ExcludedGroups

std::vector< bool >	IamVaried
	Is the ith parameter varied. More...

std::vector< std::vector< TString > >	ParamNames
	Name of parameters which we are going to analyse. More...

std::vector< std::vector< double > >	ParamCentral
	Parameters central values which we are going to analyse. More...

std::vector< std::vector< double > >	ParamErrors
	Uncertainty on a single parameter. More...

std::vector< std::vector< bool > >	ParamFlat
	Whether Param has flat prior or not. More...

std::vector< int >	nParam
	Number of parameters per type. More...

std::vector< ParameterEnum >	ParamType
	Make an enum for which class this parameter belongs to so we don't have to keep string comparing. More...

std::vector< int >	ParamTypeStartPos

std::vector< std::string >	ParameterGroup

std::vector< TString >	SampleName_v
	Vector of each systematic. More...

std::vector< TString >	SystName_v
	Vector of each sample PDF object. More...

std::string	OutputName
	Name of output files. More...

TString	CanvasName
	Name of canvas which help to save to the sample pdf. More...

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

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

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

bool	printToPDF
	Will plot all plot to PDF not only to root file. More...

bool	FancyPlotNames
	Whether we want fancy plot names or not. More...

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

bool	ApplySmoothing
	Apply smoothing for 2D histos using root algorithm. More...

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

bool	useFFTAutoCorrelation
	MJR: Use FFT-based autocorrelation algorithm (save time & resources)? More...

std::vector< int >	NDSamplesBins

std::vector< std::string >	NDSamplesNames

std::unique_ptr< TF1 >	Gauss
	Gaussian fitter. More...

TFile *	OutputFile
	The output file. More...

std::unique_ptr< TCanvas >	Posterior
	Fancy canvas used for our beautiful plots. More...

TVectorD *	Central_Value
	Vector with central value for each parameter. More...

TVectorD *	Means
	Vector with mean values using Arithmetic Mean. More...

TVectorD *	Errors
	Vector with errors values using RMS. More...

TVectorD *	Means_Gauss
	Vector with mean values using Gaussian fit. More...

TVectorD *	Errors_Gauss
	Vector with error values using Gaussian fit. More...

TVectorD *	Means_HPD
	Vector with mean values using Highest Posterior Density. More...

TVectorD *	Errors_HPD
	Vector with error values using Highest Posterior Density. More...

TVectorD *	Errors_HPD_Positive
	Vector with positive error (right hand side) values using Highest Posterior Density. More...

TVectorD *	Errors_HPD_Negative
	Vector with negative error (left hand side) values using Highest Posterior Density. More...

TMatrixDSym *	Covariance
	Posterior Covariance Matrix. More...

TMatrixDSym *	Correlation
	Posterior Correlation Matrix. More...

std::vector< TH1D * >	hpost
	Holds 1D Posterior Distributions. More...

std::vector< std::vector< TH2D * > >	hpost2D
	Holds 2D Posterior Distributions. More...

std::unique_ptr< TH2D >	hviolin
	Holds violin plot for all dials. More...

std::unique_ptr< TH2D >	hviolin_prior
	Holds prior violin plot for all dials,. More...

M3::float_t **	ParStep
	Array holding values for all parameters. More...

unsigned int *	StepNumber
	Step number for step, important if chains were merged. More...

int	nBins
	Number of bins. More...

double	DrawRange
	Drawrange for SetMaximum. More...

bool	CacheMCMC
	MCMC Chain has been cached. More...

bool	doDiagMCMC
	Doing MCMC Diagnostic. More...

int	nBatches
	Number of batches for Batched Mean. More...

int	AutoCorrLag
	LagL used in AutoCorrelation. More...

double **	BatchedAverages
	Values of batched average for every param and batch. More...

double **	SampleValues
	Holds the sample values. More...

double **	SystValues
	Holds the systs values. More...

double *	AccProbValues
	Holds all accProb. More...

double *	AccProbBatchedAverages
	Holds all accProb in batches. More...

bool	ReweightPosterior
	Whether to apply reweighting weight or not. More...

std::string	ReweightName
	Name of branch used for chain reweighting. More...

double *	WeightValue
	Stores value of weight for each step. More...

Detailed Description

This class extends MCMC and allow specialised for Oscillation parameters analysis which require specialised hardcoding.

Author: Clarence Wret; Kamil Skwarczynski

Definition at line 10 of file OscProcessor.h.

Constructor & Destructor Documentation

◆ OscProcessor()

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

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

Parameters

InputFile The path to the input file containing MCMC data.

Todo:: Here where we should add all unitarity triangles, fancy Jarlskog studies and other hacky things that only make sense for oscitations

Definition at line 14 of file OscProcessor.cpp.

                                                      : MCMCProcessor(InputFile) {
 // ****************************
   //KS: WARNING this only work when you project from Chain, will nor work when you try SetBranchAddress etc. Turn it on only if you know how to use it
   PlotJarlskog = false;
  
   Sin2Theta13Index = M3::_BAD_INT_;
   Sin2Theta12Index = M3::_BAD_INT_;
   Sin2Theta23Index = M3::_BAD_INT_;
   DeltaCPIndex = M3::_BAD_INT_;
   DeltaM2_23Index = M3::_BAD_INT_;
 }

◆ ~OscProcessor()

OscProcessor::~OscProcessor ( )

virtual

Destroys the OscProcessor object.

Definition at line 29 of file OscProcessor.cpp.

                             {
 // ****************************
 }

Member Function Documentation

◆ CalcJarlskog()

double OscProcessor::CalcJarlskog	(	const double	s2th13,
		const double	s2th23,
		const double	s2th12,
		const double	dcp
	)		const

protected

Calculate Jarlskog Invariant using oscillation parameters.

Parameters

s2th13	Value of \( \sin^2\theta_{13} \)
s2th23	Value of \( \sin^2\theta_{23} \)
s2th12	Value of \( \sin^2\theta_{12} \)
dcp	CP-violating phase \( \delta_{\text{CP}} \) (in radians)

Returns: The value of the Jarlskog invariant \( J_{\text{CP}} \) [19]

Definition at line 98 of file OscProcessor.cpp.

                                                                                                                        {
 // ***************
   const double s13  = std::sqrt(s2th13);
   const double s23  = std::sqrt(s2th23);
   const double s12  = std::sqrt(s2th12);
   const double sdcp = std::sin(dcp);
   const double c13  = std::sqrt(1.-s2th13);
   const double c12  = std::sqrt(1.-s2th12);
   const double c23  = std::sqrt(1.-s2th23);
  
   const double j = s13*c13*c13*s12*c12*s23*c23*sdcp;
  
   return j;
 }

◆ LoadAdditionalInfo()

void OscProcessor::LoadAdditionalInfo ( )

overrideprotectedvirtual

Read the Osc cov file and get the input central values and errors Here we allow Jarlskog Shenanigans.

Todo:: remove this hardcoding (e.g., use a map or enum-to-name function)

Todo:: we should actually calculate central value and prior error but leave it for now...

Reimplemented from MCMCProcessor.

Definition at line 35 of file OscProcessor.cpp.

                                       {
 // ***************
   // KS: Check if OscParams were enabled, in future we will also get
   for(size_t i = 0; i < ParameterGroup.size(); i++) {
     if(ParameterGroup[i] == "Osc"){
       OscEnabled = true;
       break;
     }
   }
  
   if(OscEnabled)
   {
     for (int i = 0; i < nDraw; ++i)
     {
       //Those keep which parameter type we run currently and relative number
       const int ParamEnum = ParamType[i];
       const int ParamNo = i - ParamTypeStartPos[ParameterEnum(ParamEnum)];
       const std::string CurrentName = ParamNames[ParamEnum][ParamNo].Data();
  
       if (CurrentName == "sin2th_13") {
         Sin2Theta13Index = i;
         Sin2Theta13Name  = CurrentName;
       } else if (CurrentName == "sin2th_12") {
         Sin2Theta12Index = i;
         Sin2Theta12Name  = CurrentName;
       } else if (CurrentName == "sin2th_23") {
         Sin2Theta23Index = i;
         Sin2Theta23Name  = CurrentName;
       } else if (CurrentName == "delta_cp") {
         DeltaCPIndex     = i;
         DeltaCPName      = CurrentName;
       } else if (CurrentName == "delm2_23") {
         DeltaM2_23Index  = i;
         DeltaM2_23Name   = CurrentName;
       }
     }
   } else{
     MACH3LOG_WARN("Didn't find oscillation parameters");
   }
  
   if(PlotJarlskog && OscEnabled)
   {
     Chain->SetAlias("J_cp", "TMath::Sqrt(sin2th_13)*TMath::Sqrt(1.-sin2th_13)*TMath::Sqrt(1.-sin2th_13)*TMath::Sqrt(sin2th_12)*TMath::Sqrt(1.-sin2th_12)*TMath::Sqrt(sin2th_23)*TMath::Sqrt(1.-sin2th_23)*TMath::Sin(delta_cp)");
     BranchNames.push_back("J_cp");
     ParamType.push_back(kXSecPar);
     nParam[kXSecPar]++;
     nDraw++;
  
     ParamCentral[kXSecPar].push_back( 0. );
     ParamErrors[kXSecPar].push_back( 1. );
     // Push back the name
     ParamNames[kXSecPar].push_back("J_cp");
     ParamFlat[kXSecPar].push_back( false );
   } else if(PlotJarlskog && !OscEnabled) {
     MACH3LOG_ERROR("Trying to enable Jarlskog without oscillations");
     throw MaCh3Exception(__FILE__,__LINE__);
   }
 }

◆ MakeJarlskogPlot()

void OscProcessor::MakeJarlskogPlot	(	const std::unique_ptr< TH1D > &	jarl,
		const std::unique_ptr< TH1D > &	jarl_flatsindcp,
		const std::unique_ptr< TH1D > &	jarl_NH,
		const std::unique_ptr< TH1D > &	jarl_NH_flatsindcp,
		const std::unique_ptr< TH1D > &	jarl_IH,
		const std::unique_ptr< TH1D > &	jarl_IH_flatsindcp
	)

protected

Perform Jarlskog Plotting.

Author: Kevin Wood

Note: based on drawJarl_dcpPriorComparison.C

Definition at line 385 of file OscProcessor.cpp.

                                                                                    {
 // ***************
   MACH3LOG_INFO("Starting {}", __func__);
   int originalErrorLevel = gErrorIgnoreLevel;
   gErrorIgnoreLevel = kFatal;
  
   // 1-->NH, 0-->both, -1-->IH
   for(int hierarchy = -1; hierarchy <= 1; hierarchy++)
   {
     std::unique_ptr<TH1D> j_hist;
     std::unique_ptr<TH1D> j_hist_sdcp;
     if(hierarchy == 1) {
       j_hist = M3::Clone(jarl_NH.get(), "");
       j_hist_sdcp = M3::Clone(jarl_NH_flatsindcp.get(), "");
       j_hist->SetTitle(";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
     } else if(hierarchy == 0) {
       j_hist = M3::Clone(jarl.get(), "");
       j_hist_sdcp = M3::Clone(jarl_flatsindcp.get(), "");
       j_hist->SetTitle(";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
     } else if(hierarchy == -1) {
       j_hist = M3::Clone(jarl_IH.get(), "");
       j_hist_sdcp = M3::Clone(jarl_IH_flatsindcp.get(), "");
       j_hist->SetTitle(";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
     } else {
       MACH3LOG_ERROR("Invalid hierarchy option. 1 for NH, 0 for both, -1 for IH");
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
  
     j_hist->Rebin(7);
     j_hist_sdcp->Rebin(7);
  
     j_hist->SetLineColor(kAzure-2);
     j_hist_sdcp->SetLineColor(kOrange+1);
     j_hist->SetLineWidth(2);
     j_hist_sdcp->SetLineWidth(2);
  
     auto StyleAxis = [](TH1* h) {
       auto xAxis = h->GetXaxis();
       auto yAxis = h->GetYaxis();
  
       xAxis->SetLabelSize(0.04);
       xAxis->SetLabelFont(132);
       xAxis->SetTitleSize(0.04);
       xAxis->SetTitleOffset(0.80);
       xAxis->SetTitleFont(132);
       xAxis->SetNdivisions(505);
       xAxis->SetTickSize(0.04);
  
       yAxis->SetLabelSize(0.04);
       yAxis->SetLabelFont(132);
       yAxis->SetTitleSize(0.04);
       yAxis->SetTitleOffset(1.2);
       yAxis->SetTitleFont(132);
       yAxis->SetNdivisions(505);
       yAxis->SetTickSize(0.04);
     };
  
     StyleAxis(j_hist.get());
  
     j_hist->GetXaxis()->SetRangeUser(-0.04,0.04);
     j_hist->Scale(1./j_hist->Integral());
     j_hist_sdcp->Scale(1./j_hist_sdcp->Integral());
  
     std::unique_ptr<TH1D> j_hist_copy = M3::Clone(j_hist.get(), "j_hist_copy");
     std::unique_ptr<TH1D> j_hist_1sig = M3::Clone(j_hist.get(), "j_hist_1sig");
     std::unique_ptr<TH1D> j_hist_2sig = M3::Clone(j_hist.get(), "j_hist_2sig");
     std::unique_ptr<TH1D> j_hist_3sig = M3::Clone(j_hist.get(), "j_hist_3sig");
  
     //upper and lower edges
     double j_bf = j_hist_copy->GetXaxis()->GetBinCenter(j_hist_copy->GetMaximumBin());
     double j_1sig_low = 9999999.;
     double j_1sig_up  = -9999999.;
     double j_2sig_low = 9999999.;;
     double j_2sig_up  = -9999999.;
     double j_3sig_low = 9999999.;;
     double j_3sig_up  = -9999999.;
  
  
     std::unique_ptr<TH1D> j_hist_sdcp_copy = M3::Clone(j_hist_sdcp.get(), "j_hist_sdcp_copy");
     std::unique_ptr<TH1D> j_hist_sdcp_1sig = M3::Clone(j_hist_sdcp.get(), "j_hist_sdcp_1sig");
     std::unique_ptr<TH1D> j_hist_sdcp_2sig = M3::Clone(j_hist_sdcp.get(), "j_hist_sdcp_2sig");
     std::unique_ptr<TH1D> j_hist_sdcp_3sig = M3::Clone(j_hist_sdcp.get(), "j_hist_sdcp_3sig");
  
     //upper and lower edges
     double j_sdcp_1sig_low = 9999999.;
     double j_sdcp_1sig_up  = -9999999.;
     double j_sdcp_2sig_low = 9999999.;;
     double j_sdcp_2sig_up  = -9999999.;
     double j_sdcp_3sig_low = 9999999.;;
     double j_sdcp_3sig_up  = -9999999.;
  
     double contlevel1 = 0.68;
     double contlevel2 = 0.90;
     double contlevel4 = 0.99;
     double contlevel5 = 0.9973;
     double integral, tsum = 0.;
  
     integral = j_hist_copy->Integral();
  
     while((tsum/integral)<contlevel5) {
       double tmax = j_hist_copy->GetMaximum();
       int bin = j_hist_copy->GetMaximumBin();
       double xval = j_hist_copy->GetXaxis()->GetBinCenter(bin);
       double xwidth = j_hist_copy->GetXaxis()->GetBinWidth(bin);
       if((tsum/integral)<contlevel1) {
         j_hist_copy->SetBinContent(bin,-1.0);
         j_hist_1sig->SetBinContent(bin,0.);
         j_hist_2sig->SetBinContent(bin,0.);
         j_hist_3sig->SetBinContent(bin,0.);
         if(xval<j_1sig_low && xval<j_bf) j_1sig_low = xval - xwidth/2.;
         if(xval>j_1sig_up && xval>j_bf) j_1sig_up = xval + xwidth/2.;
       }
       if((tsum/integral)<contlevel2  && (tsum / integral > contlevel1) ) {
         j_hist_copy->SetBinContent(bin,-5.0);
         j_hist_2sig->SetBinContent(bin,0.);
         j_hist_3sig->SetBinContent(bin,0.);
         if(xval<j_2sig_low && xval<j_bf) j_2sig_low = xval - xwidth/2.;
         if(xval>j_2sig_up && xval>j_bf) j_2sig_up = xval + xwidth/2.;
       }
       if((tsum/integral)<contlevel4  && (tsum / integral > contlevel1) ) {
         j_hist_copy->SetBinContent(bin,-9.0);
         j_hist_3sig->SetBinContent(bin,0.);
         if(xval < j_3sig_low && xval <j_bf) j_3sig_low = xval - xwidth/2.;
         if(xval > j_3sig_up && xval > j_bf) j_3sig_up = xval + xwidth/2.;
       }
       tsum+=tmax;
     }
  
     integral = j_hist_sdcp_copy->Integral();
     tsum = 0.;
  
     while((tsum/integral)<contlevel5) {
       double tmax = j_hist_sdcp_copy->GetMaximum();
       int bin = j_hist_sdcp_copy->GetMaximumBin();
       double xval = j_hist_sdcp_copy->GetXaxis()->GetBinCenter(bin);
       double xwidth = j_hist_sdcp_copy->GetXaxis()->GetBinWidth(bin);
       if((tsum/integral)<contlevel1) {
         j_hist_sdcp_copy->SetBinContent(bin,-1.0);
         j_hist_sdcp_1sig->SetBinContent(bin,0.);
         j_hist_sdcp_2sig->SetBinContent(bin,0.);
         j_hist_sdcp_3sig->SetBinContent(bin,0.);
         if(xval<j_sdcp_1sig_low && xval<j_bf) j_sdcp_1sig_low = xval - xwidth/2.;
         if(xval>j_sdcp_1sig_up && xval>j_bf) j_sdcp_1sig_up = xval + xwidth/2.;
       }
       if((tsum/integral)<contlevel2  && (tsum / integral > contlevel1) ) {
         j_hist_sdcp_copy->SetBinContent(bin,-5.0);
         j_hist_sdcp_2sig->SetBinContent(bin,0.);
         j_hist_sdcp_3sig->SetBinContent(bin,0.);
         if(xval<j_sdcp_2sig_low && xval<j_bf) j_sdcp_2sig_low = xval - xwidth/2.;
         if(xval>j_sdcp_2sig_up && xval>j_bf) j_sdcp_2sig_up = xval + xwidth/2.;
       }
       if((tsum/integral)<contlevel4  && (tsum / integral > contlevel1) ) {
         j_hist_sdcp_copy->SetBinContent(bin,-9.0);
         j_hist_sdcp_3sig->SetBinContent(bin,0.);
         if(xval<j_sdcp_3sig_low && xval<j_bf) j_sdcp_3sig_low = xval - xwidth/2.;
         if(xval>j_sdcp_3sig_up && xval>j_bf) j_sdcp_3sig_up = xval + xwidth/2.;
       }
       tsum+=tmax;
     }
  
     j_hist_1sig->SetLineStyle(9);
     j_hist_sdcp_1sig->SetLineStyle(9);
     j_hist_2sig->SetLineStyle(7);
     j_hist_sdcp_2sig->SetLineStyle(7);
     j_hist_3sig->SetLineStyle(2);
     j_hist_sdcp_3sig->SetLineStyle(2);
  
     auto ldash = std::make_unique<TH1D>("ldash", "solid",  10, -0.04, 0.04);
     auto sdash = std::make_unique<TH1D>("sdash", "dashed", 10, -0.04, 0.04);
     auto fdash = std::make_unique<TH1D>("fdash", "fdashed",10, -0.04, 0.04);
     ldash->SetLineColor(kBlack);
     sdash->SetLineColor(kBlack);
     fdash->SetLineColor(kBlack);
     ldash->SetLineWidth(2);
     sdash->SetLineWidth(2);
     fdash->SetLineWidth(2);
     ldash->SetLineStyle(9);
     sdash->SetLineStyle(7);
     fdash->SetLineStyle(2);
  
     double vertUp = 0.5 * j_hist->GetMaximum();
     auto jline_1sig_low = std::make_unique<TLine>(j_1sig_low, 0., j_1sig_low, vertUp);
     auto jline_2sig_low = std::make_unique<TLine>(j_2sig_low, 0., j_2sig_low, vertUp);
     auto jline_3sig_low = std::make_unique<TLine>(j_3sig_low, 0., j_3sig_low, vertUp);
  
     auto jline_1sig_up = std::make_unique<TLine>(j_1sig_up, 0., j_1sig_up,vertUp);
     auto jline_2sig_up = std::make_unique<TLine>(j_2sig_up, 0., j_2sig_up,vertUp);
     auto jline_3sig_up = std::make_unique<TLine>(j_3sig_up, 0., j_3sig_up,vertUp);
  
     auto jline_sdcp_1sig_low = std::make_unique<TLine>(j_sdcp_1sig_low, 0., j_sdcp_1sig_low, vertUp);
     auto jline_sdcp_2sig_low = std::make_unique<TLine>(j_sdcp_2sig_low, 0., j_sdcp_2sig_low, vertUp);
     auto jline_sdcp_3sig_low = std::make_unique<TLine>(j_sdcp_3sig_low, 0., j_sdcp_3sig_low, vertUp);
  
     auto jline_sdcp_1sig_up = std::make_unique<TLine>(j_sdcp_1sig_up, 0., j_sdcp_1sig_up, vertUp);
     auto jline_sdcp_2sig_up = std::make_unique<TLine>(j_sdcp_2sig_up, 0., j_sdcp_2sig_up, vertUp);
     auto jline_sdcp_3sig_up = std::make_unique<TLine>(j_sdcp_3sig_up, 0., j_sdcp_3sig_up, vertUp);
  
     double arrowLength = 0.003;
     double arrowHeight = vertUp;
  
     auto MakeArrow = [&](double x, Color_t color, Width_t width) -> std::unique_ptr<TArrow> {
       auto arrow = std::make_unique<TArrow>(x, arrowHeight, x - arrowLength, arrowHeight, 0.02, ">");
       arrow->SetLineColor(color);
       arrow->SetLineWidth(width);
       return arrow;
     };
  
     auto j_arrow_1sig_up = MakeArrow(j_1sig_up, j_hist_1sig->GetLineColor(), j_hist_1sig->GetLineWidth());
     auto j_arrow_2sig_up = MakeArrow(j_2sig_up, j_hist_2sig->GetLineColor(), j_hist_2sig->GetLineWidth());
     auto j_arrow_3sig_up = MakeArrow(j_3sig_up, j_hist_3sig->GetLineColor(), j_hist_3sig->GetLineWidth());
  
     auto j_sdcp_arrow_1sig_up = MakeArrow(j_sdcp_1sig_up, j_hist_sdcp_1sig->GetLineColor(), j_hist_sdcp_1sig->GetLineWidth());
     auto j_sdcp_arrow_2sig_up = MakeArrow(j_sdcp_2sig_up, j_hist_sdcp_2sig->GetLineColor(), j_hist_sdcp_2sig->GetLineWidth());
     auto j_sdcp_arrow_3sig_up = MakeArrow(j_sdcp_3sig_up, j_hist_sdcp_3sig->GetLineColor(), j_hist_sdcp_3sig->GetLineWidth());
  
     MACH3LOG_DEBUG("j_1sig_low = {:.4f}, j_2sig_low = {:.4f}, j_3sig_low = {:.4f}", j_1sig_low, j_2sig_low, j_3sig_low);
     MACH3LOG_DEBUG("j_1sig_up = {:.4f}, j_2sig_up = {:.4f}, j_3sig_up = {:.4f}", j_1sig_up, j_2sig_up, j_3sig_up);
  
     auto CopyLineStyle = [](const TH1D* src, TLine* dst) {
       dst->SetLineColor(src->GetLineColor());
       dst->SetLineStyle(src->GetLineStyle());
       dst->SetLineWidth(src->GetLineWidth());
     };
  
     CopyLineStyle(j_hist_1sig.get(), jline_1sig_low.get());
     CopyLineStyle(j_hist_1sig.get(), jline_1sig_up.get());
     CopyLineStyle(j_hist_2sig.get(), jline_2sig_low.get());
     CopyLineStyle(j_hist_2sig.get(), jline_2sig_up.get());
     CopyLineStyle(j_hist_3sig.get(), jline_3sig_low.get());
     CopyLineStyle(j_hist_3sig.get(), jline_3sig_up.get());
  
     CopyLineStyle(j_hist_sdcp_1sig.get(), jline_sdcp_1sig_low.get());
     CopyLineStyle(j_hist_sdcp_1sig.get(), jline_sdcp_1sig_up.get());
     CopyLineStyle(j_hist_sdcp_2sig.get(), jline_sdcp_2sig_low.get());
     CopyLineStyle(j_hist_sdcp_2sig.get(), jline_sdcp_2sig_up.get());
     CopyLineStyle(j_hist_sdcp_3sig.get(), jline_sdcp_3sig_low.get());
     CopyLineStyle(j_hist_sdcp_3sig.get(), jline_sdcp_3sig_up.get());
  
     auto leg = std::make_unique<TLegend>(0.45, 0.60, 0.75, 0.90);
     leg->SetTextSize(0.05);
     leg->SetFillStyle(0);
     leg->SetNColumns(1);
     leg->SetTextFont(132);
     leg->SetBorderSize(0);
  
     leg->AddEntry(j_hist.get(), "Prior flat in #delta_{CP}", "l");
     leg->AddEntry(j_hist_sdcp.get(), "Prior flat in sin#delta_{CP}", "l");
     leg->AddEntry(ldash.get(), "68% CI", "l");
     leg->AddEntry(sdash.get(), "90% CI", "l");
     leg->AddEntry(fdash.get(), "99% CI", "l");
  
     j_hist->GetYaxis()->SetRangeUser(0., j_hist->GetMaximum()*1.15);
     j_hist->Draw("h");
     j_hist_sdcp->Draw("same h");
  
     jline_sdcp_1sig_up->Draw("same");
     jline_sdcp_2sig_up->Draw("same");
     jline_sdcp_3sig_up->Draw("same");
     jline_1sig_up->Draw("same");
     jline_2sig_up->Draw("same");
     jline_3sig_up->Draw("same");
  
     j_arrow_1sig_up->Draw();
     j_arrow_2sig_up->Draw();
     j_arrow_3sig_up->Draw();
     j_sdcp_arrow_1sig_up->Draw();
     j_sdcp_arrow_2sig_up->Draw();
     j_sdcp_arrow_3sig_up->Draw();
     leg->Draw("same");
  
     auto ttext = std::make_unique<TText>();
     ttext->SetNDC(); // Use normalized device coordinates
     ttext->SetTextSize(0.03); // Adjust size as needed
     ttext->SetTextAlign(13); // Align left-top
  
     if (hierarchy == 1) ttext->DrawText(0.15, 0.85, "Normal Ordering");
     else if (hierarchy == 0) ttext->DrawText(0.15, 0.85, "Both Orderings");
     else if (hierarchy == -1) ttext->DrawText(0.15, 0.85, "Inverted Ordering");
  
     gPad->RedrawAxis();
     Posterior->Update();
     gPad->Update();
  
     Posterior->Print(CanvasName);
  
     if(hierarchy == 1)  Posterior->Write("jarl1D_NH_comp");
     else if(hierarchy == 0)  Posterior->Write("jarl1D_both_comp");
     else if(hierarchy == -1)  Posterior->Write("jarl1D_IH_comp");
   }
  
   gErrorIgnoreLevel = originalErrorLevel;
 }

◆ MakePiePlot()

void OscProcessor::MakePiePlot ( )

Make fancy Pie plot for delta CP.

Author: Asher Kaboth

Note: based on makePiePlot.C

Definition at line 684 of file OscProcessor.cpp.

                                {
 // ***************
   if(!OscEnabled || DeltaCPIndex == M3::_BAD_INT_)
   {
     MACH3LOG_WARN("Will not {}, as oscillation parameters are missing", __func__);
     return;
   }
   MACH3LOG_INFO("Starting {}", __func__);
  
   // get best fit for delta CP
   const double best_fit = (*Means_HPD)(DeltaCPIndex);
  
   const double sigma_p = (*Errors_HPD_Positive)(DeltaCPIndex);
   const double sigma_n = (*Errors_HPD_Negative)(DeltaCPIndex);
   // make sure result is between -pi and pi
   auto wrap_pi = [](double x) {
     while (x > TMath::Pi())  x -= 2*TMath::Pi();
     while (x < -TMath::Pi()) x += 2*TMath::Pi();
     return x;
   };
  
   std::array<double, 6> bds;
   bds[0] = wrap_pi(best_fit - 3.0 * sigma_n); // -3σ
   bds[1] = wrap_pi(best_fit - 2.0 * sigma_n); // -2σ
   bds[2] = wrap_pi(best_fit - 1.0 * sigma_n); // -1σ
   bds[3] = wrap_pi(best_fit + 1.0 * sigma_p); // +1σ
   bds[4] = wrap_pi(best_fit + 2.0 * sigma_p); // +2σ
   bds[5] = wrap_pi(best_fit + 3.0 * sigma_p); // +3σ
  
   constexpr double radius = 0.4;
   constexpr double rad_to_deg = 180.0 / TMath::Pi();
  
   // ROOT expects TEllipse angles in degrees, counterclockwise from the x-axis.
   // If phimax < phimin, ROOT draws counterclockwise across the full circle, causing overlaps.
   // This ensures threesigA slice stays within the intended range.
   auto normalize_angle = [](double rad) {
     // If rad is negative, add 2*pi to wrap into [0, 2*pi)
     if (rad < 0) rad += 2.0 * TMath::Pi();
     return rad;
   };
  
   TEllipse onesig   (0.5, 0.5, radius, radius, bds[2] * rad_to_deg, bds[4] * rad_to_deg);
   TEllipse twosigA  (0.5, 0.5, radius, radius, bds[1] * rad_to_deg, bds[2] * rad_to_deg);
   TEllipse twosigB  (0.5, 0.5, radius, radius, bds[3] * rad_to_deg, bds[4] * rad_to_deg);
  
   // three sigma slices
   TEllipse threesigA(0.5, 0.5, radius, radius, bds[0] * rad_to_deg, normalize_angle(bds[1]) * rad_to_deg);
   TEllipse threesigB(0.5, 0.5, radius, radius, bds[4] * rad_to_deg, bds[5] * rad_to_deg);
  
   // Remaining slices
   TEllipse rest(0.5, 0.5, radius, radius, bds[5]*rad_to_deg, bds[0]*rad_to_deg);
   TEllipse restA(0.5, 0.5, radius, radius, bds[5]*rad_to_deg, 180.0);
   TEllipse restB(0.5, 0.5, radius, radius, -180.0, bds[0]*rad_to_deg);
  
   onesig.SetFillColor(13);
   twosigA.SetFillColor(12);
   twosigB.SetFillColor(12);
   threesigA.SetFillColor(11);
   threesigB.SetFillColor(11);
   TLine line1(0.5 - radius, 0.5, 0.5 + radius, 0.5);
   line1.SetLineWidth(3);
  
   TLine line2(0.5, 0.5 - radius, 0.5, 0.5 + radius);
   line2.SetLineWidth(3);
  
   TArrow bf(0.5, 0.5, 0.5 + radius * cos(best_fit),0.5 + radius * sin(best_fit),0.04, "|>");
   bf.SetLineWidth(3);
   bf.SetLineColor(kRed);
   bf.SetFillColor(kRed);
  
   TCanvas canvas("canvas", "canvas", 0, 0, 1000, 1000);
   onesig.Draw();
   twosigA.Draw();
   twosigB.Draw();
   threesigA.Draw();
   threesigB.Draw();
  
   // Check if the rest wraps around the circle
   if (bds[5] > 0) {
     // Single rest slice
     rest.Draw();
   } else {
     // Split rest into two slices
     restA.Draw();
     restB.Draw();
   }
  
   line1.Draw();
   line2.Draw();
   bf.Draw();
  
   TLegend leg(0.0, 0.8, 0.23, 0.95);
   leg.AddEntry(&bf,        "Best Fit", "L");
   leg.AddEntry(&onesig,    "1#sigma",  "F");
   leg.AddEntry(&twosigA,   "2#sigma",  "F");
   leg.AddEntry(&threesigA, "3#sigma",  "F");
   leg.Draw();
  
   // KS: Simple lambda to avoid copy-pasting
   auto draw_text = [](auto& txt, Color_t color = kBlack) {
     txt.SetTextAlign(22);
     txt.SetTextColor(color);
     txt.SetTextFont(43);
     txt.SetTextSize(40);
     txt.SetTextAngle(0);
     txt.Draw();
   };
  
   //KS: If best fit point is somehow very close text we simply not plot it
   // Define a threshold for "too close"
   constexpr double too_close_threshold = 0.1;
  
   // Position of tbf
   const double tbf_x = 0.5 + (radius + 0.02) * cos(best_fit);
   const double tbf_y = 0.5 + (radius + 0.02) * sin(best_fit);
  
   // Function to calculate distance between two points
   auto distance = [](double x1, double y1, double x2, double y2) {
     return std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
   };
  
   // Check and draw right 0
   constexpr double t0_x = 0.5 + radius + 0.02;
   constexpr double t0_y = 0.5;
   TText t0(t0_x, t0_y, "0");
   if (distance(tbf_x, tbf_y, t0_x, t0_y) > too_close_threshold) {
     draw_text(t0);
   }
  
   // Check and draw left pi
   constexpr double tp_x = 0.5 - radius - 0.02;
   constexpr double tp_y = 0.5;
   TLatex tp(tp_x, tp_y, "#pi");
   if (distance(tbf_x, tbf_y, tp_x, tp_y) > too_close_threshold) {
     draw_text(tp);
   }
  
   // Check and draw top pi/2
   constexpr double tp2_x = 0.5;
   constexpr double tp2_y = 0.5 + radius + 0.04;
   TLatex tp2(tp2_x, tp2_y, "#frac{#pi}{2}");
   if (distance(tbf_x, tbf_y, tp2_x, tp2_y) > too_close_threshold) {
     draw_text(tp2);
   }
  
   // Check and draw bottom -pi/2
   constexpr double tmp2_x = 0.5;
   constexpr double tmp2_y = 0.5 - radius - 0.04;
   TLatex tmp2(tmp2_x, tmp2_y, "-#frac{#pi}{2}");
   if (distance(tbf_x, tbf_y, tmp2_x, tmp2_y) > too_close_threshold) {
     draw_text(tmp2);
   }
  
   TLatex tbf(0.5 + (radius + 0.02) * cos(best_fit),
              0.5 + (radius + 0.02) * sin(best_fit),
              fmt::format("{:.2f}", best_fit).c_str());
   draw_text(tbf, kRed);
  
   canvas.Print(CanvasName);
 }

◆ PerformJarlskogAnalysis()

void OscProcessor::PerformJarlskogAnalysis ( )

Perform Several Jarlskog Plotting.

Author: Kevin Wood

Note: based on makeJarlskog.C

Definition at line 137 of file OscProcessor.cpp.

                                            {
 // ***************
   if(!OscEnabled ||
     Sin2Theta13Index == M3::_BAD_INT_ ||
     Sin2Theta12Index == M3::_BAD_INT_ ||
     Sin2Theta23Index == M3::_BAD_INT_ ||
     DeltaCPIndex == M3::_BAD_INT_||
     DeltaM2_23Index == M3::_BAD_INT_)
   {
     MACH3LOG_WARN("Will not {}, as oscillation parameters are missing", __func__);
     return;
   }
   MACH3LOG_INFO("Starting {}", __func__);
  
   bool DoReweight = false;
  
   double s2th13, s2th23, s2th12, dcp, dm2 = M3::_BAD_DOUBLE_;
   double weight = 1.0;
   std::pair<double, double> Sin13_NewPrior;
  
   // Now read the MCMC file
   TFile *TempFile = M3::Open((MCMCFile + ".root"), "open", __FILE__, __LINE__);
  
   // Get the settings for the MCMC
   TMacro *Config = TempFile->Get<TMacro>("Reweight_Config");
  
   if (Config != nullptr) {
     MACH3LOG_INFO("Found Reweight_Config in chain");
     
     // Print the reweight configuration for user info
     YAML::Node Settings = TMacroToYAML(*Config); 
     // Simple check: only enable DoReweight if it's a 1D sin2th_13 Gaussian reweight since Savage Dickey process later on generates values from the Gaussian
     if(CheckNodeExists(Settings, "ReweightMCMC")) {
       YAML::Node firstReweight = Settings["ReweightMCMC"].begin()->second;
       int dimension = GetFromManager<int>(firstReweight["ReweightDim"], 1);
       std::string reweightType = GetFromManager<std::string>(firstReweight["ReweightType"], "");
       auto paramNames = GetFromManager<std::vector<std::string>>(firstReweight["ReweightVar"], {});
       if (dimension == 1 && reweightType == "Gaussian" && paramNames.size() == 1){
         Sin13_NewPrior = Get<std::pair<double, double>>(firstReweight["ReweightPrior"],__FILE__,__LINE__);
         DoReweight = true;
       } else {
         MACH3LOG_INFO("No valid reweighting configuration (1D Gaussian on sin2th_13 only) found for Jarlskog analysis");
       }
     } else {
       MACH3LOG_INFO("No reweighting configuration found for Jarlskog analysis");
     }
 }
  
   TempFile->Close();
   delete TempFile;
  
   TDirectory *JarlskogDir = OutputFile->mkdir("Jarlskog");
   JarlskogDir->cd();
  
   unsigned int step = 0;
   Chain->SetBranchStatus("*", false);
  
   Chain->SetBranchStatus(Sin2Theta13Name.c_str(), true);
   Chain->SetBranchAddress(Sin2Theta13Name.c_str(), &s2th13);
  
   Chain->SetBranchStatus(Sin2Theta23Name.c_str(), true);
   Chain->SetBranchAddress(Sin2Theta23Name.c_str(), &s2th23);
  
   Chain->SetBranchStatus(Sin2Theta12Name.c_str(), true);
   Chain->SetBranchAddress(Sin2Theta12Name.c_str(), &s2th12);
  
   Chain->SetBranchStatus(DeltaCPName.c_str(), true);
   Chain->SetBranchAddress(DeltaCPName.c_str(), &dcp);
  
   Chain->SetBranchStatus(DeltaM2_23Name.c_str(), true);
   Chain->SetBranchAddress(DeltaM2_23Name.c_str(), &dm2);
  
   Chain->SetBranchStatus("step", true);
   Chain->SetBranchAddress("step", &step);
  
   if(DoReweight) {
     Chain->SetBranchStatus("Weight", true);
     Chain->SetBranchAddress("Weight", &weight);
   } else {
     MACH3LOG_WARN("Not applying reweighting weight");
     weight = 1.0;
   }
  
   // Original histograms
   auto jarl = std::make_unique<TH1D>("jarl", "jarl", 1000, -0.05, 0.05);
   jarl->SetDirectory(nullptr);
   auto jarl_th23 = std::make_unique<TH2D>("jarl_th23", "jarl_th23", 500, -0.05, 0.05, 500, 0.3, 0.7);
   jarl_th23->SetDirectory(nullptr);
   auto jarl_dcp = std::make_unique<TH2D>("jarl_dcp", "jarl_dcp", 500, -0.05, 0.05, 500, -1. * TMath::Pi(), TMath::Pi());
   jarl_dcp->SetDirectory(nullptr);
  
   jarl->SetTitle("Jarlskog Invariant;J #equiv s_{13}c_{13}^{2}s_{12}c_{12}s_{23}c_{23}sin#delta;Posterior probability");
   jarl_th23->SetTitle("Jarlskog Invariant;J #equiv s_{13}c_{13}^{2}s_{12}c_{12}s_{23}c_{23}sin#delta;Posterior probability");
  
   // Clones
   auto jarl_IH  = M3::Clone(jarl.get(), "jarl_IH");
   auto jarl_NH  = M3::Clone(jarl.get(), "jarl_NH");
  
   auto jarl_th23_IH  = M3::Clone(jarl_th23.get(), "jarl_th23_IH");
   auto jarl_th23_NH  = M3::Clone(jarl_th23.get(), "jarl_th23_NH");
  
   auto jarl_dcp_IH  = M3::Clone(jarl_dcp.get(), "jarl_dcp_IH");
   auto jarl_dcp_NH  = M3::Clone(jarl_dcp.get(), "jarl_dcp_NH");
  
   auto jarl_flatsindcp          = M3::Clone(jarl.get(), "jarl_flatsindcp");
   auto jarl_IH_flatsindcp       = M3::Clone(jarl.get(), "jarl_IH_flatsindcp");
   auto jarl_NH_flatsindcp       = M3::Clone(jarl.get(), "jarl_NH_flatsindcp");
  
   auto jarl_th23_flatsindcp     = M3::Clone(jarl_th23.get(), "jarl_th23_flatsindcp");
   auto jarl_th23_IH_flatsindcp  = M3::Clone(jarl_th23.get(), "jarl_th23_IH_flatsindcp");
   auto jarl_th23_NH_flatsindcp  = M3::Clone(jarl_th23.get(), "jarl_th23_NH_flatsindcp");
  
   auto jarl_prior               = M3::Clone(jarl.get(), "jarl_prior");
   auto jarl_prior_flatsindcp    = M3::Clone(jarl.get(), "jarl_prior_flatsindcp");
   std::unique_ptr<TH1D> jarl_wRC_prior, jarl_wRC_prior_flatsindcp, jarl_wRC_prior_t2kth23;
   // Only use this if chain has reweigh weight [mostly coming from Reactor Constrains]
   if(DoReweight){
     jarl_wRC_prior            = M3::Clone(jarl.get(), "jarl_wRC_prior");
     jarl_wRC_prior_flatsindcp = M3::Clone(jarl.get(), "jarl_wRC_prior_flatsindcp");
     jarl_wRC_prior_t2kth23    = M3::Clone(jarl.get(), "jarl_wRC_prior_flatsindcp");
   }
  
   // to apply a prior that is flat in sin(dcp) intead of dcp
   auto prior3 = std::make_unique<TF1>("prior3", "TMath::Abs(TMath::Cos(x))");
  
   // T2K prior is flat (and uncorrelated) in dcp, sin^2(th13), sin^2(th23)
   auto randGen = std::make_unique<TRandom3>(0);
   const Long64_t countwidth = nEntries/5;
  
   for(int i = 0; i < nEntries; ++i) {
     if (i % countwidth == 0) {
       MaCh3Utils::PrintProgressBar(i, nEntries);
       MaCh3Utils::EstimateDataTransferRate(Chain, i);
     } else {
       Chain->GetEntry(i);
     }
  
     if(step < BurnInCut) continue; // burn-in cut
  
     const double j = CalcJarlskog(s2th13, s2th23, s2th12, dcp);
     const double prior_weight = prior3->Eval(dcp);
  
     jarl->Fill(j, weight);
     jarl_th23->Fill(j, s2th23, weight);
     jarl_dcp->Fill(j, dcp, weight);
  
     jarl_flatsindcp->Fill(j, prior_weight*weight);
     jarl_th23_flatsindcp->Fill(j, s2th23, prior_weight*weight);
  
     const double prior_s2th13 = SamplePriorForParam(Sin2Theta13Index, randGen, {0.,1.});
     const double prior_s2th23 = SamplePriorForParam(Sin2Theta23Index, randGen, {0.,1.});
     const double prior_s2th12 = SamplePriorForParam(Sin2Theta12Index, randGen, {0.,1.});
     const double prior_dcp = SamplePriorForParam(DeltaCPIndex, randGen, {-1.*TMath::Pi(),TMath::Pi()});
     // KS: This is hardcoded but we always assume flat in delta CP so probably fine
     const double prior_sindcp = randGen->Uniform(-1., 1.);
  
     const double prior_s13          = std::sqrt(prior_s2th13);
     const double prior_s23          = std::sqrt(prior_s2th23);
     const double prior_s12          = std::sqrt(prior_s2th12);
     const double prior_sdcp         = std::sin(prior_dcp);
     const double prior_c13          = std::sqrt(1.-prior_s2th13);
     const double prior_c12          = std::sqrt(1.-prior_s2th12);
     const double prior_c23          = std::sqrt(1.-prior_s2th23);
     const double prior_j            = prior_s13*prior_c13*prior_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sdcp;
     const double prior_flatsindcp_j = prior_s13*prior_c13*prior_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sindcp;
  
     jarl_prior->Fill(prior_j);
     jarl_prior_flatsindcp->Fill(prior_flatsindcp_j);
  
     if(DoReweight) {
       const double prior_wRC_s2th13       = randGen->Gaus(Sin13_NewPrior.first, Sin13_NewPrior.second);
       const double prior_wRC_s13          = std::sqrt(prior_wRC_s2th13);
       const double prior_wRC_c13          = std::sqrt(1.-prior_wRC_s2th13);
       const double prior_wRC_j            = prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sdcp;
       const double prior_wRC_flatsindcp_j = prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sindcp;
       const double s23                    = std::sqrt(s2th23);
       const double c23                    = std::sqrt(1.-s2th23);
  
       jarl_wRC_prior->Fill(prior_wRC_j);
       jarl_wRC_prior_flatsindcp->Fill(prior_wRC_flatsindcp_j);
       jarl_wRC_prior_t2kth23->Fill(prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*s23*c23*prior_sdcp);
     }
  
     if(dm2 > 0.) {
       jarl_NH->Fill(j, weight);
       jarl_th23_NH->Fill(j, s2th23, weight);
       jarl_dcp_NH->Fill(j, dcp, weight);
       jarl_NH_flatsindcp->Fill(j, prior_weight*weight);
       jarl_th23_NH_flatsindcp->Fill(j, s2th23, prior_weight*weight);
     }
     else if(dm2 < 0.) {
       jarl_IH->Fill(j, weight);
       jarl_th23_IH->Fill(j, s2th23, weight);
       jarl_dcp_IH->Fill(j, dcp, weight);
       jarl_IH_flatsindcp->Fill(j, prior_weight*weight);
       jarl_th23_IH_flatsindcp->Fill(j, s2th23, prior_weight*weight);
     }
   }
  
   jarl->Write("jarlskog_both");
   jarl_NH->Write("jarlskog_NH");
   jarl_IH->Write("jarlskog_IH");
   jarl_th23->Write("jarlskog_th23_both");
   jarl_th23_NH->Write("jarlskog_th23_NH");
   jarl_th23_IH->Write("jarlskog_th23_IH");
  
   jarl_dcp->Write("jarlskog_dcp_both");
   jarl_dcp_NH->Write("jarlskog_dcp_NH");
   jarl_dcp_IH->Write("jarlskog_dcp_IH");
  
  
   jarl_flatsindcp->Write("jarlskog_both_flatsindcp");
   jarl_NH_flatsindcp->Write("jarlskog_NH_flatsindcp");
   jarl_IH_flatsindcp->Write("jarlskog_IH_flatsindcp");
   jarl_th23_flatsindcp->Write("jarlskog_th23_both_flatsindcp");
   jarl_th23_NH_flatsindcp->Write("jarlskog_th23_NH_flatsindcp");
   jarl_th23_IH_flatsindcp->Write("jarlskog_th23_IH_flatsindcp");
  
   jarl_prior->Write("jarl_prior");
   jarl_prior_flatsindcp->Write("jarl_prior_flatsindcp");
   if(DoReweight) {
     jarl_wRC_prior->Write("jarl_wRC_prior");
     jarl_wRC_prior_flatsindcp->Write("jarl_wRC_prior_flatsindcp");
     jarl_wRC_prior_t2kth23->Write("jarl_wRC_prior_t2kth23");
   }
  
   MakeJarlskogPlot(jarl, jarl_flatsindcp,
                    jarl_NH, jarl_NH_flatsindcp,
                    jarl_IH, jarl_IH_flatsindcp);
  
   // Perform Savage Dickey analysis
   if(DoReweight) {
     SavageDickeyPlot(jarl, jarl_wRC_prior, "Jarlskog flat #delta_{CP}", 0);
     SavageDickeyPlot(jarl_flatsindcp, jarl_wRC_prior_flatsindcp, "Jarlskog flat sin#delta_{CP}", 0);
   } else {
     SavageDickeyPlot(jarl, jarl_prior, "Jarlskog flat #delta_{CP}", 0);
     SavageDickeyPlot(jarl_flatsindcp, jarl_prior_flatsindcp, "Jarlskog flat sin#delta_{CP}", 0);
   }
  
   JarlskogDir->Close();
   delete JarlskogDir;
  
   Chain->SetBranchStatus("*", true);
   OutputFile->cd();
 }

◆ SamplePriorForParam()

double OscProcessor::SamplePriorForParam	(	const int	paramIndex,
		const std::unique_ptr< TRandom3 > &	randGen,
		const std::vector< double > &	FlatBounds
	)		const

protected

Draw Prior value.

Definition at line 114 of file OscProcessor.cpp.

                                                                                                                                                 {
 // ***************
   TString Title = "";
   double Prior = 1.0, PriorError = 1.0;
   bool FlatPrior = false;
  
   // Get info for this parameter
   GetNthParameter(paramIndex, Prior, PriorError, Title);
  
   ParameterEnum ParType = ParamType[paramIndex];
   int ParamTemp = paramIndex - ParamTypeStartPos[ParType];
   FlatPrior = ParamFlat[ParType][ParamTemp];
  
   if (FlatPrior) {
     return randGen->Uniform(FlatBounds[0], FlatBounds[1]);
   } else {
     // Gaussian prior centered at Prior with width PriorError
     return randGen->Gaus(Prior, PriorError);
   }
 }

Member Data Documentation

◆ DeltaCPIndex

int OscProcessor::DeltaCPIndex

protected

Index of \(\delta_{\mathrm{CP}}\) in the parameter list.

Definition at line 80 of file OscProcessor.h.

◆ DeltaCPName

std::string OscProcessor::DeltaCPName

protected

Name of the parameter representing \(\delta_{\mathrm{CP}}\) (the CP-violating phase).

Definition at line 69 of file OscProcessor.h.

◆ DeltaM2_23Index

int OscProcessor::DeltaM2_23Index

protected

Index of \(\Delta m^2_{32}\) in the parameter list.

Definition at line 82 of file OscProcessor.h.

◆ DeltaM2_23Name

std::string OscProcessor::DeltaM2_23Name

protected

Name of the parameter representing \(\Delta m^2_{32}\) (mass-squared difference).

Definition at line 71 of file OscProcessor.h.

◆ OscEnabled

bool OscProcessor::OscEnabled

protected

Will plot Jarlskog Invariant using information in the chain.

Definition at line 60 of file OscProcessor.h.

◆ PlotJarlskog

bool OscProcessor::PlotJarlskog

protected

Will plot Jarlskog Invariant using information in the chain.

Definition at line 57 of file OscProcessor.h.

◆ Sin2Theta12Index

int OscProcessor::Sin2Theta12Index

protected

Index of \(\sin^2\theta_{12}\) in the parameter list.

Definition at line 76 of file OscProcessor.h.

◆ Sin2Theta12Name

std::string OscProcessor::Sin2Theta12Name

protected

Name of the parameter representing \(\sin^2\theta_{12}\).

Definition at line 65 of file OscProcessor.h.

◆ Sin2Theta13Index

int OscProcessor::Sin2Theta13Index

protected

Index of \(\sin^2\theta_{13}\) in the parameter list.

Definition at line 74 of file OscProcessor.h.

◆ Sin2Theta13Name

std::string OscProcessor::Sin2Theta13Name

protected

Name of the parameter representing \(\sin^2\theta_{13}\).

Definition at line 63 of file OscProcessor.h.

◆ Sin2Theta23Index

int OscProcessor::Sin2Theta23Index

protected

Index of \(\sin^2\theta_{23}\) in the parameter list.

Definition at line 78 of file OscProcessor.h.

◆ Sin2Theta23Name

std::string OscProcessor::Sin2Theta23Name

protected

Name of the parameter representing \(\sin^2\theta_{23}\).

Definition at line 67 of file OscProcessor.h.

The documentation for this class was generated from the following files:

Fitters/OscProcessor.h
Fitters/OscProcessor.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ OscProcessor()

◆ ~OscProcessor()

Member Function Documentation

◆ CalcJarlskog()

◆ LoadAdditionalInfo()

◆ MakeJarlskogPlot()

◆ MakePiePlot()

◆ PerformJarlskogAnalysis()

◆ SamplePriorForParam()

Member Data Documentation

◆ DeltaCPIndex

◆ DeltaCPName

◆ DeltaM2_23Index

◆ DeltaM2_23Name

◆ OscEnabled

◆ PlotJarlskog

◆ Sin2Theta12Index

◆ Sin2Theta12Name

◆ Sin2Theta13Index

◆ Sin2Theta13Name

◆ Sin2Theta23Index

◆ Sin2Theta23Name