SampleSummary Class Reference

Class to calculate pvalue produce posterior predictive and many fancy Bayesian stuff [10]. More...

#include <Fitters/SampleSummary.h>

Collaboration diagram for SampleSummary:

Public Member Functions
	SampleSummary (const int n_Samples, const std::string &Filename, SampleHandlerBase *const sample, const int nSteps)
	Constructor. More...
	~SampleSummary ()
	Destructor. More...
void	AddData (std::vector< TH2Poly * > &DataHist)
	KS: Add data histograms. More...
void	AddNominal (std::vector< TH2Poly * > &NominalHist, std::vector< TH2Poly * > &W2Nom)
	KS: Add prior histograms. More...
void	AddThrow (std::vector< TH2Poly * > &MCHist, std::vector< TH2Poly * > &W2Hist, const double LLHPenalty=0.0, const double Weight=1.0, const int DrawNumber=0)
	KS: Add histograms with throws. More...
void	AddThrowByMode (std::vector< std::vector< TH2Poly * >> &SampleVector_ByMode)
	KS: Add histograms for each mode. More...
void	Write ()
	KS: Write results into root file. More...
void	SetLikelihood (const TestStatistic TestStat)
	KS: Set likelihood type. More...
void	SetNModelParams (const int nPars)
	Set number of model params used for BIC. More...

Private Member Functions
void	MakePredictive ()
	Finalise the distributions from the thrown samples. More...
void	PrepareOutput ()
	KS: Prepare output tree and necessary variables. More...
void	CalcLLH (TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
	Helper functions to calculate likelihoods using TH2Poly, will modify MC hist title to include LLH. More...
void	CalcLLH (TH1D *const &Data, TH1D *const &MC, TH1D *const &W2)
	Helper functions to calculate likelihoods using TH1D, will modify MC hist title to include LLH. More...
double	GetLLH (TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
	Helper functions to calculate likelihoods using TH2Poly. More...
double	GetLLH (TH1D *const &Data, TH1D *const &MC, TH1D *const &W2)
	Helper functions to calculate likelihoods using TH1D. More...
void	PlotBetaParameters ()
	KS: In Barlow Beeston we have Beta Parameters which scale generated MC. More...
void	StudyKinematicCorrelations ()
	KS: Study how correlated are sample or kinematic bins. More...
void	MakeCutLLH ()
	Make the cut LLH histogram. More...
void	MakeCutLLH1D (TH1D *Histogram, double llh_ref=-999)
void	MakeCutEventRate (TH1D *Histogram, const double DataRate)
	Make the 1D Event Rate Hist. More...
void	MakeChi2Hists ()
	Make the fluctuated histograms (2D and 1D) for the chi2s Essentially taking the MCMC draws and calculating their LLH to the Posterior predictive distribution And additionally taking the data histogram and calculating the LLH to the predictive distribution Additionally we calculate the chi2 of the draws (fluctuated) of the MC with the prior/posterior predictive and plot it vs the chi2 from the draws of MCMC and the data. More...
bool	CheckSamples (const int Length)
	Check the length of samples agrees. More...
TH1D *	ProjectHist (TH2D *Histogram, const bool ProjectX)
	Helper to project TH2D onto axis. More...
TH1D *	ProjectPoly (TH2Poly *Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist=false)
	Helper to project TH2Poly onto axis. More...
void	MakeFluctuatedHistogram (TH1D *FluctHist, TH1D *PolyHist)
	Make Poisson fluctuation of TH1D hist. More...
void	MakeFluctuatedHistogram (TH2Poly *FluctHist, TH2Poly *PolyHist)
	Make Poisson fluctuation of TH2Poly hist. More...
void	StudyInformationCriterion (M3::kInfCrit Criterion)
	Information Criterion. More...
void	StudyBIC ()
	Study Bayesian Information Criterion (BIC) [11]. More...
void	StudyDIC ()
	KS: Get the Deviance Information Criterion (DIC) [29] [31]. More...
void	StudyWAIC ()
	KS: Get the Watanabe-Akaike information criterion (WAIC) [11] [15]. More...

Private Attributes
std::unique_ptr< TRandom3 >	rnd
	Random number generator. More...
bool	first_pass
	KS: Hacky flag to let us know if this is first toy. More...
bool	StandardFluctuation
	KS: We have two methods for Poissonian fluctuation. More...
std::vector< std::vector< TH2Poly * > >	MCVector
	Vector of vectors which holds the loaded MC histograms. More...
std::vector< std::vector< TH2Poly * > >	W2MCVector
	Vector of vectors which holds the loaded W2 histograms. More...
std::vector< std::vector< std::vector< TH2Poly * > > >	MCVectorByMode
	Vector of vectors which holds the loaded MC histograms for each mode. More...
std::vector< double >	LLHPenaltyVector
	Vector to hold the penalty term. More...
std::vector< double >	WeightVector
	Vector holding weight. More...
int	nSamples
	Number of samples. More...
std::vector< std::string >	SampleNames
	name for each sample More...
std::vector< std::vector< std::unique_ptr< TH1D > > >	PosteriorHist
	The posterior predictive for the whole selection: this gets built after adding in the toys. Now an array of Th1ds, 1 for each poly bin, for each sample. More...
std::vector< std::vector< std::unique_ptr< TH1D > > >	w2Hist
	The posterior predictive for the whole selection: this gets built after adding in the toys. Now an array of Th1ds, 1 for each poly bin, for each sample for W2. More...
std::vector< TH2D * >	ViolinHists_ProjectX
	Posterior predictive but for X projection but as a violin plot. More...
std::vector< TH2D * >	ViolinHists_ProjectY
	Posterior predictive but for Y projection but as a violin plot. More...
std::vector< TH2Poly * >	DataHist
	The data histogram for the selection. More...
std::vector< TH1D * >	DataHist_ProjectX
	The data histogram for the selection X projection. More...
std::vector< TH1D * >	DataHist_ProjectY
	The data histogram for the selection Y projection. More...
std::vector< TH2Poly * >	NominalHist
	The nominal histogram for the selection. More...
std::vector< TH2Poly * >	W2NomHist
	Pointer to the w2 histograms (for nominal values). More...
std::vector< TH2Poly * >	W2MeanHist
	Pointer to the w2 histograms (for mean values). More...
std::vector< TH2Poly * >	W2ModeHist
	Pointer to the w2 histograms (for mode values). More...
std::unique_ptr< TH1D >	lnLHist
	The histogram containing the lnL for each throw. More...
std::unique_ptr< TH1D >	lnLHist_drawfluc
	The lnLhist for the draw vs MC fluctuated. More...
std::unique_ptr< TH1D >	lnLHist_drawflucdraw
	The lnLhist for the draw vs draw fluctuated. More...
std::unique_ptr< TH1D >	lnLHist_drawdata
	The lnLhist for the draw vs data. More...
std::unique_ptr< TH2D >	lnLDrawHist
	The 2D lnLhist, showing (draw vs data) and (draw vs fluct), anything above y=x axis is the p-value. More...
std::unique_ptr< TH2D >	lnLFlucHist
	The 2D lnLHist, showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-value. More...
std::unique_ptr< TH2D >	lnLDrawHistRate
	The 2D lnLhist, showing (draw vs data) and (draw vs fluct), using rate, anything above y=x axis is the p-value. More...
std::unique_ptr< TH2D >	lnLFlucHist_ProjectX
	The 2D lnLHist but for ProjectionX histogram (pmu), showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-value. More...
std::vector< TH1D * >	lnLHist_Sample_DrawData
	The histogram containing the lnL (draw vs data) for each throw for each sample. More...
std::vector< TH1D * >	lnLHist_Sample_DrawflucDraw
	The histogram containing the lnL (draw vs draw fluct) for each throw for each sample. More...
std::vector< TH1D * >	lnLHist_Sample_PredflucDraw
	The histogram containing the lnL (draw vs pred fluct) for each throw for each sample. More...
std::vector< TH2Poly * >	lnLHist_Mean
	The LLH distribution in pmu cosmu for using the mean in each bin. More...
std::vector< TH2Poly * >	lnLHist_Mode
	The LLH distribution in pmu cosmu for using the mode in each bin. More...
std::vector< TH1D * >	lnLHist_Mean_ProjectX
	The LLH distribution in pmu using the mean in each bin. More...
std::vector< TH2Poly * >	MeanHist
	The posterior predictive distribution in pmu cosmu using the mean. More...
std::vector< TH2Poly * >	MeanHistCorrected
	The posterior predictive distribution in pmu cosmu using the mean after applying Barlow-Beeston Correction. More...
std::vector< TH2Poly * >	ModeHist
	The posterior predictive distribution in pmu cosmu using the mode. More...
std::vector< TH1D * >	lnLHist_Mean1D
	Holds the bin-by-bin LLH for the mean posterior predictive vs the data. More...
std::vector< TH1D * >	lnLHist_Mode1D
	Holds the bin-by-bin LLH for the mode posterior predictive vs the data. More...
std::unique_ptr< TH1D >	RandomHist
	Holds the history of which entries have been drawn in the MCMC file. More...
std::vector< std::vector< std::unique_ptr< TH1D > > >	BetaHist
	Distribution of beta parameters in Barlow Beeston formalisms. More...
bool	DoBetaParam
	Are we making Beta Histograms. More...
unsigned int	nChainSteps
	Number of throws by user. More...
bool	isPriorPredictive
	bool whether we have Prior or Posterior Predictive More...
bool	doShapeOnly
	bool whether to normalise each toy to have shape based p-value and pos pred distribution More...
unsigned int	nThrows
	Number of throws. More...
std::vector< int >	maxBins
	Max Number of Bins per each sample. More...
double	llh_total
	Total LLH for the posterior predictive distribution. More...
std::string	OutputName
	Output filename. More...
TFile *	Outputfile
	Output filename. More...
std::vector< TDirectory * >	Dir
	Directory for each sample. More...
TTree *	OutputTree
	TTree which we save useful information to. More...
std::vector< double >	llh_data_draw
	Data vs Draw. More...
std::vector< double >	llh_drawfluc_draw
	Fluctuated Draw vs Draw. More...
std::vector< double >	llh_predfluc_draw
	Fluctuated Predictive vs Draw. More...
std::vector< double >	llh_rate_data_draw
	Data vs Draw using rate only. More...
std::vector< double >	llh_rate_predfluc_draw
	Fluctuated Predictive vs Draw using rate only. More...
std::vector< double >	llh_data_drawfluc
	Data vs Fluctuated Draw. More...
std::vector< double >	llh_data_predfluc
	Data vs Fluctuated Predictive. More...
std::vector< double >	llh_draw_pred
	Draw vs Predictive. More...
std::vector< double >	llh_drawfluc_pred
	Fluctuated Draw vs Predictive. More...
std::vector< double >	llh_predfluc_pred
	Fluctuated Predictive vs Predictive. More...
std::vector< double >	llh_drawfluc_predfluc
	Fluctuated Draw vs Fluctuated Predictive. More...
std::vector< double >	llh_datafluc_draw
	Fluctuated Data vs Draw. More...
std::vector< double >	llh_data_draw_ProjectX
	Projection X (most likely muon momentum) of LLH. More...
std::vector< double >	llh_drawfluc_draw_ProjectX
double	llh_penalty
	LLH penalty for each throw. More...
double	total_llh_data_draw
	Data vs Draw. More...
double	total_llh_drawfluc_draw
	Fluctuated Draw vs Draw. More...
double	total_llh_predfluc_draw
	Fluctuated Predictive vs Draw. More...
double	total_llh_rate_data_draw
	Rate Data vs Draw. More...
double	total_llh_rate_predfluc_draw
	Fluctuated Predictive vs Draw using Rate. More...
double	total_llh_data_predfluc
	Data vs Fluctuated Predictive. More...
double	total_llh_data_drawfluc
	Data vs Fluctuated Draw. More...
double	total_llh_draw_pred
	Draw vs Predictive. More...
double	total_llh_drawfluc_pred
	Fluctuated Draw vs Predictive. More...
double	total_llh_drawfluc_predfluc
	Fluctuated Draw vs Fluctuated Predictive. More...
double	total_llh_datafluc_draw
	Fluctuated Data vs Draw. More...
double	total_llh_predfluc_pred
	Fluctuated Predictive vs Predictive. More...
double	total_llh_data_draw_ProjectX
	Data vs Draw for projection X (most likely muon momentum) More...
double	total_llh_drawfluc_draw_ProjectX
	Fluctuated Draw vs Draw for projection X (most likely muon momentum) More...
bool	DoByModePlots
	By mode variables. More...
std::vector< std::vector< TH2Poly * > >	MeanHist_ByMode
	The posterior predictive distribution in pmu cosmu using the mean. More...
TH1D ****	PosteriorHist_ByMode
	Histogram which corresponds to each bin in the sample's th2poly. More...
SampleHandlerBase *	SampleHandler
	Pointer to SampleHandler object, mostly used to get sample names, binning etc. More...
MaCh3Modes *	Modes
	MaCh3 Modes. More...
TestStatistic	likelihood
	Type of likelihood for example Poisson, Barlow-Beeston or Ice Cube. More...
int	nModelParams
	Number of parameters. More...
int	Debug
	Tells Debug level to save additional histograms. More...

Detailed Description

Class to calculate pvalue produce posterior predictive and many fancy Bayesian stuff [10].

For more information, visit the Wiki.

Author

Clarence Wret

Kamil Skwarczynski

Definition at line 22 of file SampleSummary.h.

Constructor & Destructor Documentation

SampleSummary()

SampleSummary::SampleSummary	(	const int	n_Samples,
		const std::string &	Filename,
		SampleHandlerBase *const	sample,
		const int	nSteps
	)

Constructor.

Parameters

n_Samples	total number of samples
Filename	name of output file
sample	pointer to sample PDF object
nChainSteps	number of steps in a chain, 0 indicate prior predictive was used

Definition at line 9 of file SampleSummary.cpp.

                                                                                                                              {
// *******************
  MACH3LOG_DEBUG("Making sample summary class...");
  #ifdef MULTITHREAD
  MACH3LOG_DEBUG("With OpenMP and {} threads", omp_get_max_threads());
  #endif
  
  StandardFluctuation = true;
  
  if(StandardFluctuation) MACH3LOG_INFO("Using standard method of statistical fluctuation");
  else MACH3LOG_INFO("Using alternative method of statistical fluctuation, which is much slower");
  
  //KS: If true it will print posterior predictive for every beta parameter it is quite useful but make root big number of plots
  DoBetaParam = true;
  if(DoBetaParam) MACH3LOG_INFO("I will calculate #beta parameters from Barlow-Beeston formalism");
 
  //If true code will normalise each histogram, this way you can calculate shape only error. etc. pvalue will be completely wrong unfortunately
  doShapeOnly = false;
 
  nChainSteps = nSteps;
  //KS: nChainSteps == 0 means we run PriorPredcitive
  if(nChainSteps == 0) isPriorPredictive = true;
  else isPriorPredictive = false;
  
  OutputName = Filename;
  nSamples = n_Samples;
  SampleHandler = sample;
 
  //Get mach3 modes from manager
  Modes = SampleHandler->GetMaCh3Modes();
 
  nThrows = 0;
  first_pass = true;
  Outputfile = nullptr;
  OutputTree = nullptr;
  rnd = std::make_unique<TRandom3>();
 
  DataHist.resize(nSamples);
  DataHist_ProjectX.resize(nSamples);
  DataHist_ProjectY.resize(nSamples);
  NominalHist.resize(nSamples);
  PosteriorHist.resize(nSamples);
  W2NomHist.resize(nSamples);
  w2Hist.resize(nSamples);
 
  ViolinHists_ProjectX.resize(nSamples);
  ViolinHists_ProjectY.resize(nSamples);
    
  if(DoBetaParam) BetaHist.resize(nSamples);
 
  maxBins.resize(nSamples);
  
  lnLHist_Mean.resize(nSamples);
  lnLHist_Mode.resize(nSamples);
  lnLHist_Mean_ProjectX.resize(nSamples);
  MeanHist.resize(nSamples);;
  if(DoBetaParam) MeanHistCorrected.resize(nSamples);;
  ModeHist.resize(nSamples);
  W2MeanHist.resize(nSamples);
  W2ModeHist.resize(nSamples);
  lnLHist_Mean1D.resize(nSamples);
  lnLHist_Mode1D.resize(nSamples);
  lnLHist_Sample_DrawData.resize(nSamples);
  lnLHist_Sample_DrawflucDraw.resize(nSamples);
  lnLHist_Sample_PredflucDraw.resize(nSamples);
    
  //KS: When a histogram is created with an axis lower limit greater or equal to its upper limit ROOT will automatically adjust histogram range
  // https://root.cern.ch/doc/master/classTH1.html#auto-bin
  lnLHist = std::make_unique<TH1D>("lnLHist_predpredfluc", "lnLHist_predpredfluc", 100, 1, -1);
  lnLHist->SetDirectory(nullptr);
  lnLHist->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Pred)");
  lnLHist->GetYaxis()->SetTitle("Counts");
 
  lnLHist_drawdata = std::make_unique<TH1D>("lnLHist_drawdata", "lnLHist_drawdata", 100, 1, -1);
  lnLHist_drawdata->SetDirectory(nullptr);
  lnLHist_drawdata->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
  lnLHist_drawdata->GetYaxis()->SetTitle("Counts");
 
  lnLHist_drawfluc = std::make_unique<TH1D>("lnLHist_drawpredfluc", "lnLHist_drawpredfluc", 100, 1, -1);
  lnLHist_drawfluc->SetDirectory(nullptr);
  lnLHist_drawfluc->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
  lnLHist_drawfluc->GetYaxis()->SetTitle("Counts");
 
  lnLHist_drawflucdraw = std::make_unique<TH1D>("lnLHist_drawflucdraw", "lnLHist_drawflucdraw", 100, 1, -1);
  lnLHist_drawflucdraw->SetDirectory(nullptr);
  lnLHist_drawflucdraw->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
  lnLHist_drawflucdraw->GetYaxis()->SetTitle("Counts");
 
  lnLDrawHist = std::make_unique<TH2D>("lnLDrawHist", "lnLDrawHist", 50, 1, -1, 50, 1, -1);
  lnLDrawHist->SetDirectory(nullptr);
  lnLDrawHist->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
  lnLDrawHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
 
  lnLFlucHist = std::make_unique<TH2D>("lnLFlucHist", "lnLFlucHist", 50, 1, -1, 50, 1, -1);
  lnLFlucHist->SetDirectory(nullptr);
  lnLFlucHist->GetXaxis()->SetTitle("-2LLH_{Draw Fluc, Draw}");
  lnLFlucHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
 
  lnLDrawHistRate = std::make_unique<TH2D>("lnLDrawHistRate", "lnLDrawHistRate", 50, 1, -1, 50, 1, -1);
  lnLDrawHistRate->SetDirectory(nullptr);
  lnLDrawHistRate->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
  lnLDrawHistRate->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
 
  //KS: This is silly as it assumes all samples uses same kinematics
  lnLFlucHist_ProjectX = std::make_unique<TH2D>("lnLFlucHist_ProjectX", "lnLFlucHist_ProjectX", 50, 1, -1, 50, 1, -1);
  lnLFlucHist_ProjectX->SetDirectory(nullptr);
  lnLFlucHist_ProjectX->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
  lnLFlucHist_ProjectX->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
  
  // Holds the hist of random number draws, only works for posterior predictive
  if(!isPriorPredictive)
  {
    RandomHist = std::make_unique<TH1D>("RandomHist", "RandomHist", 100, 0, nChainSteps);
    RandomHist->SetDirectory(nullptr);
    RandomHist->GetXaxis()->SetTitle("Step");
    const double binwidth = nChainSteps/RandomHist->GetNbinsX();
    std::stringstream ss;
    ss << "Draws/" << binwidth;
    RandomHist->GetYaxis()->SetTitle(ss.str().c_str());
    RandomHist->SetLineWidth(2);
  }
  else RandomHist = nullptr;
 
  for (int i = 0; i < nSamples; ++i)
  {
    DataHist[i] = nullptr;
    DataHist_ProjectX[i] = nullptr;
    DataHist_ProjectY[i] = nullptr;
    NominalHist[i] = nullptr;
    
    MeanHist[i] = nullptr;
    if(DoBetaParam) MeanHistCorrected[i] = nullptr;
    W2MeanHist[i] = nullptr;
    W2ModeHist[i] = nullptr;
    lnLHist_Mean[i] = nullptr;
    lnLHist_Mode[i] = nullptr;
    lnLHist_Mean_ProjectX[i] = nullptr;
    lnLHist_Mean1D[i] = nullptr;
    lnLHist_Mode1D[i] = nullptr;
    lnLHist_Sample_DrawData[i] = nullptr;
    lnLHist_Sample_DrawflucDraw[i] = nullptr;
    lnLHist_Sample_PredflucDraw[i] = nullptr;
  }//end loop over samples
 
  DoByModePlots = false;
  PosteriorHist_ByMode = nullptr;
 
  nModelParams = 0;
 
  Debug = 0;
}

~SampleSummary()

SampleSummary::~SampleSummary

(

)

Destructor.

Definition at line 163 of file SampleSummary.cpp.

                              {
// *******************
  Outputfile->cd();
 
  //ROOT is weird and once you write TFile claim ownership of histograms. Best is to first delete histograms and then close file
  Outputfile->Close();
  delete Outputfile;
 
  if(DoByModePlots)
  {
    for (int i = 0; i < nSamples; ++i)
    {
      if(DataHist[i] == nullptr) continue;
      for (int j = 0; j < Modes->GetNModes()+1; j++)
      {
        for (int k = 1; k <= maxBins[i]; ++k)
        {
          if(PosteriorHist_ByMode[i][j][k] != nullptr) delete PosteriorHist_ByMode[i][j][k];
        }
        delete[] PosteriorHist_ByMode[i][j];
        if(MeanHist_ByMode[i][j] != nullptr) delete MeanHist_ByMode[i][j];
      }
      delete[] PosteriorHist_ByMode[i];
    }
    delete[] PosteriorHist_ByMode;
  }
 
  for (int i = 0; i < nSamples; ++i)
  {
    if(DataHist[i] == nullptr) continue;
    if(DataHist[i] != nullptr) delete DataHist[i];
    if(NominalHist[i] != nullptr) delete NominalHist[i];
    if(MeanHist[i] != nullptr) delete MeanHist[i];
    if(ModeHist[i] != nullptr) delete ModeHist[i];
    if(DoBetaParam && MeanHistCorrected[i] != nullptr) delete MeanHistCorrected[i];
    if(W2MeanHist[i] != nullptr) delete W2MeanHist[i];
    if(W2ModeHist[i] != nullptr) delete W2ModeHist[i];
    
    if(ViolinHists_ProjectX[i] != nullptr) delete ViolinHists_ProjectX[i];
    if(ViolinHists_ProjectY[i] != nullptr) delete ViolinHists_ProjectY[i];
    
    if(lnLHist_Mean[i] != nullptr) delete lnLHist_Mean[i];
    if(lnLHist_Mode[i] != nullptr) delete lnLHist_Mode[i];
    if(lnLHist_Mean_ProjectX[i] != nullptr) delete lnLHist_Mean_ProjectX[i];
    if(lnLHist_Mean1D[i] != nullptr) delete lnLHist_Mean1D[i];
    if(lnLHist_Mode1D[i] != nullptr) delete lnLHist_Mode1D[i];
    if(lnLHist_Sample_DrawData[i] != nullptr) delete lnLHist_Sample_DrawData[i];
    if(lnLHist_Sample_DrawflucDraw[i] != nullptr) delete lnLHist_Sample_DrawflucDraw[i];
    if(lnLHist_Sample_PredflucDraw[i] != nullptr) delete lnLHist_Sample_PredflucDraw[i];
  }
}

Member Function Documentation

AddData()

void SampleSummary::AddData

(

std::vector< TH2Poly * > &

DataHist

)

KS: Add data histograms.

Parameters

DataHist

Histogram with data even rates for each sample

Definition at line 233 of file SampleSummary.cpp.

                                                     {
// *******************
  const int Length = int(Data.size());
  // Check length of samples are OK
  if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
  for (int i = 0; i < Length; ++i) {
    if (Data[i] == nullptr) {
      DataHist[i] = nullptr;
      DataHist_ProjectX[i] = nullptr;
      DataHist_ProjectY[i] = nullptr;
      maxBins[i] = 0;
    } else {
      std::string classname = std::string(DataHist[i]->Class_Name());
      if(classname == "TH2Poly")
      {
        DataHist[i] = static_cast<TH2Poly*>(Data[i]->Clone());
        if(doShapeOnly) NormaliseTH2Poly(DataHist[i]);
        DataHist_ProjectX[i] = ProjectPoly(DataHist[i], true, i);
        DataHist_ProjectY[i] = ProjectPoly(DataHist[i], false, i);
        maxBins[i] = DataHist[i]->GetNumberOfBins();
      } else {
        MACH3LOG_ERROR("Somehow sample {} doesn't use TH2Poly", SampleHandler->GetSampleName(i));
        MACH3LOG_ERROR("Right now I only support TH2Poly but I am ambitious piece of code and surely will have more support in the future");
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
    }
  }
}

AddNominal()

void SampleSummary::AddNominal	(	std::vector< TH2Poly * > &	NominalHist,
		std::vector< TH2Poly * > &	W2Nom
	)

KS: Add prior histograms.

Definition at line 264 of file SampleSummary.cpp.

                                                                                       {
// *******************
  const int Length = int(Nominal.size());
  if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
  
  //KS: ROOT is super annoying and you cannot use clone with openMP, hence we have another loop below
  for (int i = 0; i < Length; ++i) 
  {
    if (Nominal[i] == nullptr) {
      NominalHist[i] = nullptr;
      W2NomHist[i] = nullptr;
      lnLHist_Mean[i] = nullptr;
      lnLHist_Mode[i] = nullptr;
      lnLHist_Mean_ProjectX[i] = nullptr;
      MeanHist[i] = nullptr;
      if(DoBetaParam) MeanHistCorrected[i] = nullptr;
      ModeHist[i] = nullptr;
      W2MeanHist[i] = nullptr;
      W2ModeHist[i] = nullptr;
      lnLHist_Sample_DrawData[i] = nullptr;
      lnLHist_Sample_DrawflucDraw[i] = nullptr;
      lnLHist_Sample_PredflucDraw[i] = nullptr;
    // If not nullptr it indicates the selection was turned on, so initialise the privates
    } else {
      NominalHist[i] = static_cast<TH2Poly*>(Nominal[i]->Clone());
      if(doShapeOnly) NormaliseTH2Poly(NominalHist[i]);
      W2NomHist[i] = static_cast<TH2Poly*>(NomW2[i]->Clone());
      
      lnLHist_Mean[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      lnLHist_Mean[i]->SetDirectory(nullptr);
      lnLHist_Mode[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      lnLHist_Mode[i]->SetDirectory(nullptr);
      lnLHist_Mean_ProjectX[i] = static_cast<TH1D*>(DataHist_ProjectX[i]->Clone());
      MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      if(DoBetaParam) MeanHistCorrected[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      W2MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
      W2ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
    }
  }
  
  // Loop over the length of nominal and set the initial distributions up
  //KS: Don't multithread, mostly due to fact that we initialise histograms
  for (int i = 0; i < Length; ++i) {
  // If NULL it indicates the selection was turned off, so initialise all the hists to NULL
    if (Nominal[i] != nullptr)
    {
      std::string name = std::string(NominalHist[i]->GetName());
      name = name.substr(0, name.find("_nom"));
 
      PosteriorHist[i].resize(maxBins[i]+1);
      w2Hist[i].resize(maxBins[i]+1);
 
      if(DoBetaParam) BetaHist[i].resize(maxBins[i]+1);
 
      for (int j = 0; j <= maxBins[i]; ++j)
      {
        PosteriorHist[i][j] = nullptr;
      }
      lnLHist_Mean[i]->SetNameTitle((name+"_MeanlnL").c_str(), (name+"_MeanlnL").c_str());
      lnLHist_Mean[i]->Reset("");
      lnLHist_Mean[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
      
      lnLHist_Mode[i]->SetNameTitle((name+"_ModelnL").c_str(), (name+"_ModelnL").c_str());
      lnLHist_Mode[i]->Reset("");
      lnLHist_Mode[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
 
      lnLHist_Mean_ProjectX[i]->SetNameTitle((name+"_MeanlnL_ProjectX").c_str(), (name+"_MeanlnL_ProjectX").c_str()); 
      lnLHist_Mean_ProjectX[i]->Reset("");
      lnLHist_Mean_ProjectX[i]->GetYaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
 
      MeanHist[i]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
      MeanHist[i]->Reset("");
      MeanHist[i]->GetZaxis()->SetTitle("Mean");
      
      if(DoBetaParam)
      {
        MeanHistCorrected[i]->SetNameTitle((name+"_mean_corrected").c_str(), (name+"_mean_corrected").c_str());
        MeanHistCorrected[i]->Reset("");
        MeanHistCorrected[i]->GetZaxis()->SetTitle("Mean");
      }
      std::vector<double> xbins;
      std::vector<double> ybins;
 
      SampleHandler->SetupBinning(M3::int_t(i), xbins, ybins);
      
      //KS: Y axis is number of events to get estimate of maximal number we use integral
      const int MaxBinning = doShapeOnly ? 1 : int(NoOverflowIntegral(NominalHist[i])/4);
      ViolinHists_ProjectX[i] = new TH2D((name+"_Violin_ProjectX").c_str(), (name+"_Violin_ProjectX").c_str(), int(xbins.size()-1), &xbins[0] , 400, 0, MaxBinning);
      ViolinHists_ProjectX[i]->GetYaxis()->SetTitle("Events");
      ViolinHists_ProjectX[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetXaxis()->GetTitle()).c_str() );
      ViolinHists_ProjectX[i]->SetDirectory(nullptr);
 
      ViolinHists_ProjectY[i] = new TH2D((name+"_Violin_ProjectY").c_str(), (name+"_Violin_ProjectY").c_str(), int(ybins.size()-1), &ybins[0] , 400, 0, MaxBinning);
      ViolinHists_ProjectY[i]->GetYaxis()->SetTitle("Events");
      ViolinHists_ProjectY[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetYaxis()->GetTitle()).c_str());
      ViolinHists_ProjectY[i]->SetDirectory(nullptr);
 
      ModeHist[i]->SetNameTitle((name+"_mode").c_str(), (name+"_mode").c_str());
      ModeHist[i]->Reset("");
      ModeHist[i]->GetZaxis()->SetTitle("Mode");
 
      W2MeanHist[i]->SetNameTitle((name+"_w2mean").c_str(), (name+"_w2mean").c_str());
      W2MeanHist[i]->Reset("");
      W2MeanHist[i]->GetZaxis()->SetTitle("W2Mean");
 
      W2ModeHist[i]->SetNameTitle((name+"_w2mode").c_str(), (name+"_w2mode").c_str());
      W2ModeHist[i]->Reset("");
      W2ModeHist[i]->GetZaxis()->SetTitle("W2Mode");
 
      // Declare the lnL histograms
      lnLHist_Mean1D[i] = new TH1D((name+"_MeanlnL1D").c_str(),(name+"_MeanlnL1D").c_str(), 50, 1, -1);
      lnLHist_Mean1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
      lnLHist_Mean1D[i]->GetYaxis()->SetTitle("Counts");
 
      lnLHist_Mode1D[i] = new TH1D((name+"_ModelnL1D").c_str(),(name+"_ModelnL1D").c_str(), 50, 1, -1);
      lnLHist_Mode1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
      lnLHist_Mode1D[i]->GetYaxis()->SetTitle("Counts");
 
      lnLHist_Sample_DrawData[i] = new TH1D((name+"_lnLdrawdata").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
      lnLHist_Sample_DrawData[i]->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
      lnLHist_Sample_DrawData[i]->GetYaxis()->SetTitle("Counts");
      
      lnLHist_Sample_DrawflucDraw[i] = new TH1D((name+"_lnLdrawfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
      lnLHist_Sample_DrawflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
      lnLHist_Sample_DrawflucDraw[i]->GetYaxis()->SetTitle("Counts");
      
      lnLHist_Sample_PredflucDraw[i] = new TH1D((name+"_lnLpredfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
      lnLHist_Sample_PredflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
      lnLHist_Sample_PredflucDraw[i]->GetYaxis()->SetTitle("Counts");
    }
  }
  //KS: Separate loop for thread safe reasons
  for (int i = 0; i < Length; ++i)
  {
    //KS: We copy histograms so delete original
    delete Nominal[i];
    delete NomW2[i];
  }
}

AddThrow()

void SampleSummary::AddThrow	(	std::vector< TH2Poly * > &	MCHist,
		std::vector< TH2Poly * > &	W2Hist,
		const double	LLHPenalty = 0.0,
		const double	Weight = 1.0,
		const int	DrawNumber = 0
	)

KS: Add histograms with throws.

Definition at line 408 of file SampleSummary.cpp.

                                                                                                                                                              {
// *******************
  nThrows++;
  //KS: Only make sense for PosteriorPredictive
  if( !isPriorPredictive )RandomHist->Fill(DrawNumber);
 
  const int size = int(SampleVector.size());
  if (!CheckSamples(size)) throw MaCh3Exception(__FILE__ , __LINE__ );
 
  // Push back the throw
  MCVector.push_back(SampleVector);
  LLHPenaltyVector.push_back(LLHPenalty);
  WeightVector.push_back(Weight);
  W2MCVector.push_back(W2Vec);
 
  // Initialise the posterior hist
  if (first_pass)
  {
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      const int nXBins = 500;
      //Initialise TH1D which corresponds to each bin in the sample's th2poly
      std::string name = std::string(SampleVector[SampleNum]->GetName());
      for (int i = 1; i <= maxBins[SampleNum]; ++i)
      {
        //Get PolyBin
        TH2PolyBin* bin = static_cast<TH2PolyBin*>(SampleVector[SampleNum]->GetBins()->At(i-1));
 
        // Just make a little fancy name
        std::stringstream ss2;
        ss2 << name << "_";
        ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
        ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
 
        PosteriorHist[SampleNum][i] = std::make_unique<TH1D>(ss2.str().c_str(), ss2.str().c_str(),nXBins, 1, -1);
        PosteriorHist[SampleNum][i]->SetDirectory(nullptr);
        w2Hist[SampleNum][i] = std::make_unique<TH1D>(("w2_"+ss2.str()).c_str(), ("w2_"+ss2.str()).c_str(),nXBins, 1, -1);
        w2Hist[SampleNum][i]->SetDirectory(nullptr);
        if(DoBetaParam)
        {
          std::string betaName = "#beta_param_";
          BetaHist[SampleNum][i] = std::make_unique<TH1D>((betaName + ss2.str()).c_str(), (betaName + ss2.str()).c_str(), 70, 1, -1);
          BetaHist[SampleNum][i]->SetDirectory(nullptr);
          BetaHist[SampleNum][i]->GetXaxis()->SetTitle("#beta parameter value");
          BetaHist[SampleNum][i]->GetYaxis()->SetTitle("Counts");
        }
      } //end loop over bins
    }//end loop over samples
  }
  first_pass = false;
 
  // Loop over the samples
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
  {
    if (SampleVector[SampleNum] == nullptr) continue;
    if(doShapeOnly) NormaliseTH2Poly(SampleVector[SampleNum]);
    // Loop over the distribution and fill the prior/posterior predictive
    for (int i = 1; i <= maxBins[SampleNum]; ++i) {
      const double Content = SampleVector[SampleNum]->GetBinContent(i);
      PosteriorHist[SampleNum][i]->Fill(Content, Weight);
      const double w2 = W2Vec[SampleNum]->GetBinContent(i);
      w2Hist[SampleNum][i]->Fill(w2, Weight);
      if(DoBetaParam)
      {
        const double data = DataHist[SampleNum]->GetBinContent(i);
        const double BetaParam = GetBetaParameter(data, Content, w2, likelihood);
        BetaHist[SampleNum][i]->Fill(BetaParam, Weight);
      }
    } // end bin loop
  } // end samples loop
} // end AddThrow

AddThrowByMode()

void SampleSummary::AddThrowByMode

(

std::vector< std::vector< TH2Poly * >> &

SampleVector_ByMode

)

KS: Add histograms for each mode.

Definition at line 486 of file SampleSummary.cpp.

                                                                                      {
// *******************
  MCVectorByMode.push_back(SampleVector_ByMode);
 
  //KS: This means this is first time
  if(!DoByModePlots)
  {
    MACH3LOG_INFO("Turning reaction breadkwon mode, brum brum");
    PosteriorHist_ByMode = new TH1D***[nSamples];
    MeanHist_ByMode.resize(nSamples);
    for (int SampleNum = 0;  SampleNum < nSamples; SampleNum++)
    {
      if (DataHist[SampleNum] == nullptr) continue;
 
      PosteriorHist_ByMode[SampleNum] = new TH1D**[Modes->GetNModes()+1];
      MeanHist_ByMode[SampleNum].resize(Modes->GetNModes()+1);
      for (int j = 0; j < Modes->GetNModes()+1; j++)
      {
        PosteriorHist_ByMode[SampleNum][j] = new TH1D*[maxBins[SampleNum]+1];
        constexpr int nXBins = 500;
 
        std::string name = std::string(NominalHist[SampleNum]->GetName());
        name = name.substr(0, name.find("_nom"));
        name = name + "_"+Modes->GetMaCh3ModeName(j);
 
        for (int i = 1; i <= maxBins[SampleNum]; i++)
        {
          //Get PolyBin
          TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i-1));
 
          // Just make a little fancy name
          std::stringstream ss2;
          ss2 << name << "_";
          ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
          ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
 
          //Initialise TH1D which corresponds to each bin in the sample's th2poly
          PosteriorHist_ByMode[SampleNum][j][i] = new TH1D((name+ss2.str()).c_str(),(name+ss2.str()).c_str(),nXBins, 1, -1);
        }
        MeanHist_ByMode[SampleNum][j] = static_cast<TH2Poly*>(NominalHist[SampleNum]->Clone());
        MeanHist_ByMode[SampleNum][j]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
        MeanHist_ByMode[SampleNum][j]->Reset("");
        MeanHist_ByMode[SampleNum][j]->GetZaxis()->SetTitle("Mean");
      }
    }
  }
  DoByModePlots = true;    
  // Loop over the sameples
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int SampleNum = 0;  SampleNum < nSamples; SampleNum++)
  {
    if (DataHist[SampleNum] == nullptr) continue;
    
    for (int j = 0; j < Modes->GetNModes()+1; j++)
    {
      if(doShapeOnly) NormaliseTH2Poly(SampleVector_ByMode[SampleNum][j]);
      // Loop over the distribution and fill the prior/posterior predictive
      for (int i = 1; i <= maxBins[SampleNum]; ++i)
      {
        const double Content = SampleVector_ByMode[SampleNum][j]->GetBinContent(i);
        const int Entries = int(PosteriorHist_ByMode[SampleNum][j][i]->GetEntries());
        PosteriorHist_ByMode[SampleNum][j][i]->Fill(Content, WeightVector[Entries]);
      }
    }
  }
} // end AddThrowByMode

CalcLLH() [1/2]

void SampleSummary::CalcLLH ( TH1D *const &  Data, TH1D *const &  MC, TH1D *const &  W2  )

inlineprivate

Helper functions to calculate likelihoods using TH1D, will modify MC hist title to include LLH.

Parameters

Data	histogram with data distribution for a single sample
MC	histogram with MC distribution for a single sample
W2	histogram with W2 distribution for a single sample

Definition at line 1561 of file SampleSummary.cpp.

                                                                                                {
// ****************
  const double llh = GetLLH(DatHist, MCHist, W2Hist);
  std::stringstream ss;
  ss << "_2LLH=" << llh;
  MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
  MACH3LOG_INFO("{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), DatHist->Integral(), MCHist->Integral(), llh);
}

CalcLLH() [2/2]

void SampleSummary::CalcLLH ( TH2Poly *const &  Data, TH2Poly *const &  MC, TH2Poly *const &  W2  )

inlineprivate

Helper functions to calculate likelihoods using TH2Poly, will modify MC hist title to include LLH.

Parameters

Data	histogram with data distribution for a single sample
MC	histogram with MC distribution for a single sample
W2	histogram with W2 distribution for a single sample

Definition at line 1572 of file SampleSummary.cpp.

                                                                                                         {
// ****************
  const double llh = GetLLH(DatHist, MCHist, W2Hist);
  std::stringstream ss;
  ss << "_2LLH=" << llh;
  MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
  MACH3LOG_INFO("{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), NoOverflowIntegral(DatHist), NoOverflowIntegral(MCHist), llh);
}

CheckSamples()

bool SampleSummary::CheckSamples ( const int  Length )

inlineprivate

Check the length of samples agrees.

Definition at line 217 of file SampleSummary.cpp.

                                                 {
// *******************
  bool ok = (nSamples == Length);
  if (!ok) {
    MACH3LOG_ERROR("Size of SampleVector input != number of defined samples");
    MACH3LOG_ERROR("Size of SampleVector:  {}", Length);
    MACH3LOG_ERROR("Size of defined samples:  {}", nSamples);
    MACH3LOG_ERROR("Something has gone wrong with making the Samples");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  return ok;
}

GetLLH() [1/2]

double SampleSummary::GetLLH ( TH1D *const &  Data, TH1D *const &  MC, TH1D *const &  W2  )

inlineprivate

Helper functions to calculate likelihoods using TH1D.

Parameters

Data	histogram with data distribution for a single sample
MC	histogram with MC distribution for a single sample
W2	histogram with W2 distribution for a single sample

Definition at line 1597 of file SampleSummary.cpp.

                                                                                                 {
// ****************
  double llh = 0.0;
  for (int i = 1; i <= DatHist->GetXaxis()->GetNbins(); ++i)
  {
    const double data = DatHist->GetBinContent(i);
    const double mc = MCHist->GetBinContent(i);
    const double w2 = W2Hist->GetBinContent(i);
    llh += SampleHandler->GetTestStatLLH(data, mc, w2);
  }
  //KS: do times 2 because banff reports chi2
  return 2*llh;
}

GetLLH() [2/2]

double SampleSummary::GetLLH ( TH2Poly *const &  Data, TH2Poly *const &  MC, TH2Poly *const &  W2  )

inlineprivate

Helper functions to calculate likelihoods using TH2Poly.

Parameters

Data	histogram with data distribution for a single sample
MC	histogram with MC distribution for a single sample
W2	histogram with W2 distribution for a single sample

Definition at line 1582 of file SampleSummary.cpp.

                                                                                                          {
// ****************
  double llh = 0.0;
  for (int i = 1; i < DatHist->GetNumberOfBins()+1; ++i)
  {
    const double data = DatHist->GetBinContent(i);
    const double mc = MCHist->GetBinContent(i);
    const double w2 = W2Hist->GetBinContent(i);
    llh += SampleHandler->GetTestStatLLH(data, mc, w2);
  }
  //KS: do times 2 because banff reports chi2
  return 2*llh;
}

MakeChi2Hists()

void SampleSummary::MakeChi2Hists ( )

inlineprivate

Make the fluctuated histograms (2D and 1D) for the chi2s Essentially taking the MCMC draws and calculating their LLH to the Posterior predictive distribution And additionally taking the data histogram and calculating the LLH to the predictive distribution Additionally we calculate the chi2 of the draws (fluctuated) of the MC with the prior/posterior predictive and plot it vs the chi2 from the draws of MCMC and the data.

Definition at line 1134 of file SampleSummary.cpp.

                                  {
// *******************
  MACH3LOG_INFO("Making the chi2 histograms...");
  // Have this to signify if we're doing the first pass
  first_pass = true;
 
  double AveragePenalty = 0;
 
  // Vectors to hold exact LLH
  std::vector<double> LLH_PredFluc_V(nThrows);
  std::vector<double> LLH_DataDraw_V(nThrows);
  std::vector<double> LLH_DrawFlucDraw_V(nThrows);
    
  // Loop over the draws
  // Should look into multi-threading this. Would require temporary THxx structures from arrays
  //KS: Update above would be ideal but currently we loop over samples (see loop below) which isn't as efficient as loop over throws but much much easier to implement
  for (unsigned int i = 0; i < nThrows; ++i)
  {
    if (i % (nThrows/10) == 0) {
      MaCh3Utils::PrintProgressBar(i, nThrows);
    }
 
    // Set the total LLH to zero to initialise
    double total_llh_data_draw_temp = 0.0;
    double total_llh_drawfluc_draw_temp = 0.0;
    double total_llh_predfluc_draw_temp = 0.0;
 
    double total_llh_rate_data_draw_temp = 0.0;
    double total_llh_rate_predfluc_draw_temp = 0.0;
 
    double total_llh_data_drawfluc_temp = 0.0;
    double total_llh_data_predfluc_temp = 0.0;
    double total_llh_draw_pred_temp = 0.0;
    double total_llh_drawfluc_pred_temp = 0.0;
    double total_llh_drawfluc_predfluc_temp = 0.0;
    double total_llh_predfluc_pred_temp = 0.0;
    double total_llh_datafluc_draw_temp = 0.0;
 
    double total_llh_data_draw_ProjectX_temp = 0.0;
    double total_llh_drawfluc_draw_ProjectX_temp = 0.0;
 
    // Save the double that gets written to file
    llh_penalty = LLHPenaltyVector[i];
    AveragePenalty += llh_penalty;
 
    // Make the Poisson fluctuated hist
    std::vector<TH2Poly*> FluctHist(nSamples);
    // Also Poisson fluctuate the drawn MCMC hist
    std::vector<TH2Poly*> FluctDrawHist(nSamples);
    // Finally Poisson fluctuate the data histogram
    std::vector<TH2Poly*> DataFlucHist(nSamples);
 
    // Finally Poisson fluctuate the data histogram
    std::vector<TH1D*> FluctDrawHistProjectX(nSamples);
    std::vector<TH1D*> DrawHistProjectX(nSamples);
    std::vector<TH1D*> DrawHistProjectY(nSamples);
    std::vector<TH1D*> DrawW2HistProjectX(nSamples);
 
    //KS: We have to clone histograms here to avoid cloning in MP loop, we have to make sure binning matches, content doesn't have to
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      FluctHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
      FluctDrawHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
      DataFlucHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
 
      FluctDrawHistProjectX[SampleNum] = static_cast<TH1D*>(DataHist_ProjectX[SampleNum]->Clone());
 
      // Get the ith draw for the jth sample
      TH2Poly *DrawHist = MCVector[i][SampleNum];
      TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
 
      //ProjectPoly calls new TH1D under the hood, never define new ROOT object under MP...
      DrawHistProjectX[SampleNum] = ProjectPoly(DrawHist, true, SampleNum);
      DrawW2HistProjectX[SampleNum] = ProjectPoly(DrawW2Hist, true, SampleNum);
      DrawHistProjectY[SampleNum] = ProjectPoly(DrawHist, false, SampleNum);
    }
    #ifdef MULTITHREAD 
    //KS: might be most obscure OMP reduction I have ever made...
    #pragma omp parallel for reduction(+:total_llh_data_draw_temp, total_llh_drawfluc_draw_temp, total_llh_predfluc_draw_temp, total_llh_rate_data_draw_temp, total_llh_rate_predfluc_draw_temp, total_llh_data_drawfluc_temp, total_llh_data_predfluc_temp, total_llh_draw_pred_temp, total_llh_drawfluc_pred_temp, total_llh_drawfluc_predfluc_temp, total_llh_predfluc_pred_temp, total_llh_datafluc_draw_temp, total_llh_data_draw_ProjectX_temp, total_llh_drawfluc_draw_ProjectX_temp)
    #endif
    // Loop over the samples
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      // Get the ith draw for the jth sample
      TH2Poly *DrawHist = MCVector[i][SampleNum];
      TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
      // Skip empty samples
      if (DrawHist == nullptr) continue;
 
      // Add LLH penalties from the systematics to the LLH that use the drawn histogram
      // Data vs Draw
      llh_data_draw[SampleNum] = llh_penalty;
      // Fluctuated Draw vs Draw
      llh_drawfluc_draw[SampleNum] = llh_penalty;
      // Fluctuated Predicitve vs Draw
      llh_predfluc_draw[SampleNum] = llh_penalty;
      
      // Data vs Draw using rate
      llh_rate_data_draw[SampleNum] = llh_penalty;
      // Fluctuated Predicitve vs Draw using rate
      llh_rate_predfluc_draw[SampleNum] = llh_penalty;
 
       // Data vs Fluctuated Draw
      llh_data_drawfluc[SampleNum] = llh_penalty;
      // Draw vs Predictive
      llh_draw_pred[SampleNum] = llh_penalty;
      // Fluctuated Draw vs Predictive
      llh_drawfluc_pred[SampleNum] = llh_penalty;
      // Fluctuated Draw vs Fluctuated Predictive
      llh_drawfluc_predfluc[SampleNum] = llh_penalty;
      // Fluctuated Data vs Draw
      llh_datafluc_draw[SampleNum] = llh_penalty;
      
      //Some LLH for 1D projections
      llh_data_draw_ProjectX[SampleNum] = llh_penalty;
      llh_drawfluc_draw_ProjectX[SampleNum] = llh_penalty;
      
      //Other get 0 penalty term
      // Fluctuated Predictive vs Predictive
      llh_predfluc_pred[SampleNum] = 0.0;
      // Data vs Fluctuated Predictive
      llh_data_predfluc[SampleNum] = 0.0;
 
      // Make the Poisson fluctuated hist
      MakeFluctuatedHistogram(FluctHist[SampleNum], MeanHist[SampleNum]);
      // Also Poisson fluctuate the drawn MCMC hist
      MakeFluctuatedHistogram(FluctDrawHist[SampleNum], DrawHist);
      // Finally Poisson fluctuate the data histogram
      MakeFluctuatedHistogram(DataFlucHist[SampleNum], DataHist[SampleNum]);
 
      // Likelihood between the drawn histogram and the data
      const double DataDrawLLH = GetLLH(DataHist[SampleNum], DrawHist, DrawW2Hist);
      llh_data_draw[SampleNum] += DataDrawLLH;
      total_llh_data_draw_temp += DataDrawLLH;
      
      // Likelihood between drawn histogram and fluctuated drawn histogram
      const double DrawFlucDrawLLH = GetLLH(FluctDrawHist[SampleNum], DrawHist, DrawW2Hist);
      llh_drawfluc_draw[SampleNum] += DrawFlucDrawLLH;
      total_llh_drawfluc_draw_temp += DrawFlucDrawLLH;
      
      // Likelihood between drawn histogram and fluctuated posterior predictive distribution
      const double PredFlucDrawLLH = GetLLH(FluctHist[SampleNum], DrawHist, DrawW2Hist);
      llh_predfluc_draw[SampleNum] += PredFlucDrawLLH;
      total_llh_predfluc_draw_temp += PredFlucDrawLLH;
 
//Rate Based p-value
      // Likelihood between the drawn histogram and the data
      const double RateDataDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(DataHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
      llh_rate_data_draw[SampleNum] += RateDataDrawLLH;
      total_llh_rate_data_draw_temp += RateDataDrawLLH;
 
      // Likelihood between drawn histogram and fluctuated posterior predictive distribution using rate
      const double RatePredFlucDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(FluctHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
      llh_rate_predfluc_draw[SampleNum] += RatePredFlucDrawLLH;
      total_llh_rate_predfluc_draw_temp += RatePredFlucDrawLLH;
 
//    All LLH below are for validation reason but not used for final P-Value
      // Likelihood between the fluctuated drawn histogram and the data
      const double DataDrawFlucLLH = GetLLH(DataHist[SampleNum], FluctDrawHist[SampleNum], DrawW2Hist);
      llh_data_drawfluc[SampleNum] += DataDrawFlucLLH;
      total_llh_data_drawfluc_temp += DataDrawFlucLLH;
 
      // Likelihood between the drawn histogram and the data
      const double DataPredFlucLLH = GetLLH(DataHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
      llh_data_predfluc[SampleNum] += DataPredFlucLLH;
      total_llh_data_predfluc_temp += DataPredFlucLLH;
 
      // Likelihood between the drawn hist and the Posterior Predictive
      const double DrawPredLLH = GetLLH(DrawHist, MeanHist[SampleNum], W2MeanHist[SampleNum]);
      llh_draw_pred[SampleNum] += DrawPredLLH;
      total_llh_draw_pred_temp += DrawPredLLH;
 
      // Likelihood between fluctuated drawn and predictive
      const double DrawFlucPredLLH = GetLLH(FluctDrawHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
      llh_drawfluc_pred[SampleNum]  += DrawFlucPredLLH;
      total_llh_drawfluc_pred_temp  += DrawFlucPredLLH;
 
      // Likelihood between drawn histogram and fluctuated drawn histogram
      const double DrawFlucPredFlucLLH = GetLLH(FluctDrawHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
      llh_drawfluc_predfluc[SampleNum] += DrawFlucPredFlucLLH;
      total_llh_drawfluc_predfluc_temp += DrawFlucPredFlucLLH;
 
      // Likelihood between the fluctuated drawn histogram and the posterior predictive
      const double PredFlucPredLLH = GetLLH(FluctHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
      llh_predfluc_pred[SampleNum] += PredFlucPredLLH;
      total_llh_predfluc_pred_temp += PredFlucPredLLH;
 
      // Likelihood between fluctuated data histogram and drawn histogram 
      const double DataFlucDrawLLH = GetLLH(DataFlucHist[SampleNum], DrawHist, DrawW2Hist);
      llh_datafluc_draw[SampleNum] += DataFlucDrawLLH;
      total_llh_datafluc_draw_temp += DataFlucDrawLLH;
 
      lnLHist_Sample_DrawData[SampleNum]->Fill(DataDrawLLH);
      lnLHist_Sample_DrawflucDraw[SampleNum]->Fill(DrawFlucDrawLLH);
      lnLHist_Sample_PredflucDraw[SampleNum]->Fill(PredFlucDrawLLH);
      
//    At the end we leave LLH for 1D projections
      MakeFluctuatedHistogram(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum]);
      
      // Likelihood between the drawn histogram and the data for muon momentum
      const double DataDrawLLH_ProjectX = GetLLH(DataHist_ProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
      llh_data_draw_ProjectX[SampleNum] += DataDrawLLH_ProjectX;
      total_llh_data_draw_ProjectX_temp += DataDrawLLH_ProjectX;
      
      const double DrawFlucDrawLLH_ProjectX = GetLLH(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
      llh_drawfluc_draw_ProjectX[SampleNum] += DrawFlucDrawLLH_ProjectX;
      total_llh_drawfluc_draw_ProjectX_temp += DrawFlucDrawLLH_ProjectX;
      
      //KS: This might seem complicated but we make X and Y projection for each sample. Then we add this to the Violin plot making nice Gaussian in each kineatmical bin of x and y axis
      FastViolinFill(ViolinHists_ProjectX[SampleNum], DrawHistProjectX[SampleNum]);
      FastViolinFill(ViolinHists_ProjectY[SampleNum], DrawHistProjectY[SampleNum]);
    } // End loop over samples (still looping throws)
 
    // Delete the temporary histograms
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      delete FluctHist[SampleNum];
      delete FluctDrawHist[SampleNum];
      delete DataFlucHist[SampleNum];
      delete FluctDrawHistProjectX[SampleNum];
      delete DrawHistProjectX[SampleNum];
      delete DrawHistProjectY[SampleNum];
      delete DrawW2HistProjectX[SampleNum];
    }
 
    total_llh_data_draw = total_llh_data_draw_temp;
    total_llh_drawfluc_draw = total_llh_drawfluc_draw_temp;
    total_llh_predfluc_draw = total_llh_predfluc_draw_temp;
 
    total_llh_rate_data_draw = total_llh_rate_data_draw_temp;
    total_llh_rate_predfluc_draw = total_llh_rate_predfluc_draw_temp;
 
    total_llh_data_drawfluc = total_llh_data_drawfluc_temp;
    total_llh_data_predfluc = total_llh_data_predfluc_temp;
    total_llh_draw_pred = total_llh_draw_pred_temp;
    total_llh_drawfluc_pred = total_llh_drawfluc_pred_temp;
    total_llh_drawfluc_predfluc = total_llh_drawfluc_predfluc_temp;
    total_llh_predfluc_pred = total_llh_predfluc_pred_temp;
    total_llh_datafluc_draw = total_llh_datafluc_draw_temp;
 
    total_llh_data_draw_ProjectX = total_llh_data_draw_ProjectX_temp;
    total_llh_drawfluc_draw_ProjectX = total_llh_drawfluc_draw_ProjectX_temp;
 
    // Add LLH penalties from the systematics to the LLH that use the drawn histogram
    total_llh_data_draw     += llh_penalty;
    total_llh_drawfluc_draw += llh_penalty;
    total_llh_predfluc_draw += llh_penalty;
    //Rate based
    total_llh_rate_data_draw += llh_penalty;
    total_llh_rate_predfluc_draw += llh_penalty;
 
    total_llh_data_drawfluc += llh_penalty;
    total_llh_draw_pred     += llh_penalty;
    total_llh_drawfluc_pred += llh_penalty;
    total_llh_drawfluc_predfluc += llh_penalty;
 
    total_llh_data_draw_ProjectX += llh_penalty;
    total_llh_drawfluc_draw_ProjectX += llh_penalty;
 
    lnLHist->Fill(total_llh_predfluc_pred);
    lnLHist_drawdata->Fill(total_llh_data_draw);
    lnLHist_drawfluc->Fill(total_llh_predfluc_draw);
    lnLHist_drawflucdraw->Fill(total_llh_drawfluc_draw);
 
    lnLDrawHist->Fill(total_llh_predfluc_draw, total_llh_data_draw);
    lnLFlucHist->Fill(total_llh_drawfluc_draw, total_llh_data_draw);
 
    lnLDrawHistRate->Fill(total_llh_rate_predfluc_draw, total_llh_rate_data_draw);
 
    lnLFlucHist_ProjectX->Fill(total_llh_drawfluc_draw_ProjectX, total_llh_data_draw_ProjectX);
    
    // Also save to arrays to make sure we have the utmost super accuracy
    LLH_PredFluc_V[i] = total_llh_predfluc_draw;
    LLH_DataDraw_V[i] = total_llh_data_draw;
    LLH_DrawFlucDraw_V[i] = total_llh_drawfluc_draw;
 
    // Write to the output tree
    OutputTree->Fill();
  } // End loop over throws
 
  AveragePenalty = AveragePenalty/double(nThrows);
  MACH3LOG_INFO("Average LLH penalty over toys is {:.2f}", AveragePenalty);
  // Calculate exact p-value instead of binned
  unsigned int Accept_PredFluc = 0;
  unsigned int Accept_DrawFluc = 0;
  for (unsigned int i = 0; i < nThrows; ++i)
  {
    if (LLH_DataDraw_V[i] > LLH_DrawFlucDraw_V[i]) Accept_DrawFluc++;
    if (LLH_DataDraw_V[i] > LLH_PredFluc_V[i]) Accept_PredFluc++;
  }
  const double pvalue_DrawFluc = double(Accept_DrawFluc)/double(nThrows);
  const double pvalue_PredFluc = double(Accept_PredFluc)/double(nThrows);
 
  MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Draw: {:.2f}", pvalue_DrawFluc);
  MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Prediction: {:.2f}", pvalue_PredFluc);
}

MakeCutEventRate()

void SampleSummary::MakeCutEventRate ( TH1D *  Histogram, const double  DataRate  )

inlineprivate

Make the 1D Event Rate Hist.

Definition at line 1511 of file SampleSummary.cpp.

                                                                           {
// ****************
  // For the event rate histogram add a TLine to the data rate
  auto TempLine = std::make_unique<TLine>(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum());
  TempLine->SetLineColor(kRed);
  TempLine->SetLineWidth(2);
  // Also fit a Gaussian because why not?
  TF1 *Fitter = new TF1("Fit", "gaus", Histogram->GetBinLowEdge(1), Histogram->GetBinLowEdge(Histogram->GetNbinsX()+1));
  Histogram->Fit(Fitter, "RQ");
  Fitter->SetLineColor(kRed-5);
  // Calculate a p-value
  double Above = 0.0;
  for (int z = 0; z < Histogram->GetNbinsX(); ++z) {
    const double xvalue = Histogram->GetBinCenter(z+1);
    if (xvalue >= DataRate) {
      Above += Histogram->GetBinContent(z+1);
    }
  }
  const double pvalue = Above/Histogram->Integral();
  auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
  Legend->SetFillColor(0);
  Legend->SetFillStyle(0);
  Legend->SetLineWidth(0);
  Legend->SetLineColor(0);
  Legend->AddEntry(TempLine.get(), Form("Data, %.0f, p-value=%.2f", DataRate, pvalue), "l");
  Legend->AddEntry(Histogram, Form("MC, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
  Legend->AddEntry(Fitter, Form("Gauss, #mu=%.1f#pm%.1f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
  std::string TempTitle = std::string(Histogram->GetName());
  TempTitle += "_canv";
  TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
  TempCanvas->SetGridx();
  TempCanvas->SetGridy();
  TempCanvas->SetRightMargin(0.03);
  TempCanvas->SetBottomMargin(0.08);
  TempCanvas->SetLeftMargin(0.10);
  TempCanvas->SetTopMargin(0.06);
  TempCanvas->cd();
  Histogram->Draw();
  TempLine->Draw("same");
  Fitter->Draw("same");
  Legend->Draw("same");
  TempCanvas->Write();
  Histogram->Write();
 
  delete TempCanvas;
  delete Fitter;
}

MakeCutLLH()

void SampleSummary::MakeCutLLH ( )

inlineprivate

Make the cut LLH histogram.

Definition at line 1433 of file SampleSummary.cpp.

                               {
// *******************
  Outputfile->cd();
  MakeCutLLH1D(lnLHist.get());
  MakeCutLLH1D(lnLHist_drawfluc.get());
  MakeCutLLH1D(lnLHist_drawdata.get());
  MakeCutLLH1D(lnLHist_drawflucdraw.get());
  
  Get2DBayesianpValue(lnLDrawHist.get());
  Get2DBayesianpValue(lnLFlucHist.get());
  Get2DBayesianpValue(lnLDrawHistRate.get());
  Get2DBayesianpValue(lnLFlucHist_ProjectX.get());
}

MakeCutLLH1D()

void SampleSummary::MakeCutLLH1D ( TH1D *  Histogram, double  llh_ref = -999  )

inlineprivate

Definition at line 1449 of file SampleSummary.cpp.

                                                                {
// ****************
  const double TotalIntegral = Histogram->Integral();
  double Above = 0.0;
  // Get the LLH reference from total llh or some reference histogram
  double llh_reference = 0.0;
  if (llh_ref >= 0) {
    llh_reference = llh_ref;
  } else {
    llh_reference = llh_total;
  }
  for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) {
    const double xvalue = Histogram->GetBinCenter(i+1);
    if (xvalue >= llh_reference) {
      Above += Histogram->GetBinContent(i+1);
    }
  }
  const double pvalue = Above/TotalIntegral;
  std::stringstream ss;
  ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue;
  Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str());
 
  // Write a TCanvas and make a line and a filled histogram
  auto TempLine = std::make_unique<TLine>(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum());
  TempLine->SetLineColor(kBlack);
  TempLine->SetLineWidth(2);
 
  // Make the fill histogram
  TH1D *TempHistogram = static_cast<TH1D*>(Histogram->Clone());
  TempHistogram->SetFillStyle(1001);
  TempHistogram->SetFillColor(kRed);
  for (int i = 0; i < TempHistogram->GetNbinsX(); ++i) {
    if (TempHistogram->GetBinCenter(i+1) < llh_reference) {
      TempHistogram->SetBinContent(i+1, 0.0);
    }
  }
 
  auto Legend = std::make_unique<TLegend>(0.6, 0.6, 0.9, 0.9);
  Legend->SetFillColor(0);
  Legend->SetFillStyle(0);
  Legend->SetLineWidth(0);
  Legend->SetLineColor(0);
  Legend->AddEntry(TempLine.get(), Form("Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue), "l");
  Legend->AddEntry(Histogram, Form("LLH, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
  std::string Title = Histogram->GetName();
  Title += "_canv";
  TCanvas *TempCanvas = new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024);
  TempCanvas->SetGridx();
  TempCanvas->SetGridy();
  Histogram->Draw();
  TempHistogram->Draw("same");
  TempLine->Draw("same");
  Legend->Draw("same");
 
  TempCanvas->Write();
 
  delete TempHistogram;
  delete TempCanvas;
}

MakeFluctuatedHistogram() [1/2]

void SampleSummary::MakeFluctuatedHistogram ( TH1D *  FluctHist, TH1D *  PolyHist  )

inlineprivate

Make Poisson fluctuation of TH1D hist.

Definition at line 2183 of file SampleSummary.cpp.

                                                                          {
// ****************
  if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
  else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
}

MakeFluctuatedHistogram() [2/2]

void SampleSummary::MakeFluctuatedHistogram ( TH2Poly *  FluctHist, TH2Poly *  PolyHist  )

inlineprivate

Make Poisson fluctuation of TH2Poly hist.

Definition at line 2191 of file SampleSummary.cpp.

                                                                                {
// ****************
  if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
  else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
}

MakePredictive()

void SampleSummary::MakePredictive ( )

inlineprivate

Finalise the distributions from the thrown samples.

Definition at line 978 of file SampleSummary.cpp.

                                   {
// *******************
  // First make the projection on the z axis of the TH3D* for every pmu cosmu bin
  double llh_total_temp = 0.0;
 
  // Loop over the samples
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:llh_total_temp)
  #endif
  for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
  {
    // Skip disabled samples
    if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
 
    // Count the -2LLH for each histogram
    double negLogL_Mean = 0.0;
    double negLogL_Mode = 0.0;
 
    // Loop over each pmu cosmu bin
    for (int j = 1; j < maxBins[SampleNum]+1; ++j)
    {
      TH1D *Projection = PosteriorHist[SampleNum][j].get();
      TH1D *W2Projection = w2Hist[SampleNum][j].get();
 
      // Data content for the j,kth bin
      const double nData = DataHist[SampleNum]->GetBinContent(j);
      
      // Get the mean for this projection for all the samples
      // This is the mean prediction for this given j,k bin
      const double nMean = Projection->GetMean();
      const double nMeanError = Projection->GetRMS();
      const double nMode = Projection->GetBinCenter(Projection->GetMaximumBin());
      const double nModeError = GetModeError(Projection);
      
      const double nW2Mean = W2Projection->GetMean();
      const double nW2Mode = W2Projection->GetBinCenter(W2Projection->GetMaximumBin());
      
      double TempLLH_Mean = 0.0;
      double TempLLH_Mode = 0.0;
      
      //KS:Get LLH contribution getTestStatLLH can calculate Barlow Beeston/IceCube or Poisson
      TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
      TempLLH_Mode = SampleHandler->GetTestStatLLH(nData, nMode, nW2Mode);
 
      // Increment -2LLH
      //KS: do times 2 because banff reports chi2
      negLogL_Mean += 2*TempLLH_Mean;
      negLogL_Mode += 2*TempLLH_Mode;
      
      // Set the content and error to the mean in the bin
      MeanHist[SampleNum]->SetBinContent(j, MeanHist[SampleNum]->GetBinContent(j)+nMean);
      // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
      MeanHist[SampleNum]->SetBinError(j, nMeanError);
 
      if(DoBetaParam)
      {
        TH1D *BetaTemp = BetaHist[SampleNum][j].get();
        const double nBetaMean = BetaTemp->GetMean();
        const double nBetaMeanError = BetaTemp->GetRMS();
        //KS: Here we modify predictions by beta parameter from Barlow-Beeston
        MeanHistCorrected[SampleNum]->SetBinContent(j, MeanHistCorrected[SampleNum]->GetBinContent(j)+nMean*nBetaMean);
        //KS: Use error propagation to calcuate error
        const double ErrorTemp = std::sqrt( (nBetaMean*nMeanError) * (nBetaMean*nMeanError) + (nMean*nBetaMeanError) * (nMean*nBetaMeanError));
        // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
        MeanHistCorrected[SampleNum]->SetBinError(j, ErrorTemp);
      }
      // Set the content to the mode in the bin
      ModeHist[SampleNum]->SetBinContent(j, ModeHist[SampleNum]->GetBinContent(j)+nMode);
      // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
      ModeHist[SampleNum]->SetBinError(j, nModeError);
      // Set the content to the mean in the bin
      W2MeanHist[SampleNum]->SetBinContent(j, W2MeanHist[SampleNum]->GetBinContent(j)+nW2Mean);
      // Set the content to the mode in the bin
      W2ModeHist[SampleNum]->SetBinContent(j, W2ModeHist[SampleNum]->GetBinContent(j)+nW2Mode);
      
      // Set the mean and average LLH for this given bin
      // Can use these hists to see where the largest -2LLH hists come from
      lnLHist_Mean[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
      lnLHist_Mode[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mode);
      
      lnLHist_Mean1D[SampleNum]->Fill(2.0*TempLLH_Mean);
      lnLHist_Mode1D[SampleNum]->Fill(2.0*TempLLH_Mode);
    } // End loop over bins
    if(DoByModePlots)
    {
      for (int j = 0; j < Modes->GetNModes()+1; j++)
      {
        // Loop over each pmu cosmu bin
        for (int i = 1; i < maxBins[SampleNum]+1; ++i)
        {
          // Make the posterior/prior predictive projection on z
          // The z axis of Predictive is the bin content
          // Essentially zooming in on one bin and looking at the mean and mode of that bin
          TH1D *Projection = PosteriorHist_ByMode[SampleNum][j][i];
 
          // Get the mean for this projection for all the samples
          const double nMean = Projection->GetMean();
          const double nMeanError = Projection->GetRMS();
 
          // Set the content and error to the mean in the bin
          MeanHist_ByMode[SampleNum][j]->SetBinContent(i, MeanHist_ByMode[SampleNum][j]->GetBinContent(i)+nMean);
          // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
          MeanHist_ByMode[SampleNum][j]->SetBinError(i, nMeanError);
        } // End loop over bins
      }
    }
    llh_total_temp += negLogL_Mean;
  } // End loop over samples
 
  // This is not multithreaded as due to ProjectPoly it is not safe
  for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
  {
    // Skip disabled samples
    if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
 
    //KS:: Might consider caching it as we use it once agian much later
    TH1D *MeanProjectX = ProjectPoly(MeanHist[SampleNum], true, SampleNum, true);
    TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[SampleNum], true, SampleNum);
    // Loop over each pmu bin for 1D projection
    for (int j = 1; j <= lnLHist_Mean_ProjectX[SampleNum]->GetXaxis()->GetNbins(); ++j)
    {
      // Data content for the j,kth bin
      const double nData = DataHist_ProjectX[SampleNum]->GetBinContent(j);
      const double nMean = MeanProjectX->GetBinContent(j);
      const double nW2Mean = W2MeanProjectX->GetBinContent(j);
 
      double TempLLH_Mean = 0.0;
      TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
 
      //KS: do times 2 because banff reports chi2
      lnLHist_Mean_ProjectX[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
    }// End loop over  bins
 
    delete MeanProjectX;
    delete W2MeanProjectX;
  } // End loop over samples
 
  llh_total = llh_total_temp;
  // Now we have our posterior predictive histogram and it's LLH
  MACH3LOG_INFO("Prior/Posterior predictive LLH mean (sample only) = {:.2f}", llh_total);
  std::stringstream ss;
  ss << llh_total;
  lnLHist->SetTitle((std::string(lnLHist->GetTitle())+"_"+ss.str()).c_str());
 
  // Now make the fluctuated hists of the MeanHist and ModeHist
  MakeChi2Hists();
 
  // Get the 1D LLH dists
  MakeCutLLH();
} // End MakePredictive() function

PlotBetaParameters()

void SampleSummary::PlotBetaParameters ( )

inlineprivate

KS: In Barlow Beeston we have Beta Parameters which scale generated MC.

Definition at line 1612 of file SampleSummary.cpp.

                                       {
// ****************
  // Make a new directory
  TDirectory *BetaDir = Outputfile->mkdir("BetaParameters");
  BetaDir->cd();
 
  int originalErrorLevel = gErrorIgnoreLevel;
 
  //To avoid Warning in <Fit>: Fit data is empty
  gErrorIgnoreLevel = kFatal;
 
  MACH3LOG_INFO("Writing Beta parameters");
  std::vector<TDirectory *> DirBeta(nSamples);
  for (int i = 0; i < nSamples; ++i)
  {
    // Make a new directory
    DirBeta[i] = BetaDir->mkdir((SampleNames[i]).c_str());
    DirBeta[i]->cd();
 
    // Loop over each pmu cosmu bin
    for (int j = 1; j < maxBins[i]+1; ++j)
    {
      const double data = DataHist[i]->GetBinContent(j);
      const double mc = NominalHist[i]->GetBinContent(j);
      const double w2 = W2NomHist[i]->GetBinContent(j);
 
      const double BetaPrior = GetBetaParameter(data, mc, w2, likelihood);
 
      auto TempLine = std::unique_ptr<TLine>(new TLine(BetaPrior, BetaHist[i][j]->GetMinimum(), BetaPrior, BetaHist[i][j]->GetMaximum()));
      TempLine->SetLineColor(kRed);
      TempLine->SetLineWidth(2);
 
      // Also fit a Gaussian because why not?
      TF1 *Fitter = new TF1("Fit", "gaus", BetaHist[i][j]->GetBinLowEdge(1), BetaHist[i][j]->GetBinLowEdge(BetaHist[i][j]->GetNbinsX()+1));
      BetaHist[i][j]->Fit(Fitter, "RQ");
      Fitter->SetLineColor(kRed-5);
 
      auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
      Legend->SetFillColor(0);
      Legend->SetFillStyle(0);
      Legend->SetLineWidth(0);
      Legend->SetLineColor(0);
      Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.4f, N_{data}=%.0f", BetaPrior, data), "l");
      Legend->AddEntry(BetaHist[i][j].get(), Form("Post, #mu=%.4f#pm%.4f", BetaHist[i][j]->GetMean(), BetaHist[i][j]->GetRMS()), "l");
      Legend->AddEntry(Fitter, Form("Gauss, #mu=%.4f#pm%.4f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
      std::string TempTitle = std::string(BetaHist[i][j]->GetName());
 
      TempTitle += "_canv";
      TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
      TempCanvas->SetGridx();
      TempCanvas->SetGridy();
      TempCanvas->SetRightMargin(0.03);
      TempCanvas->SetBottomMargin(0.08);
      TempCanvas->SetLeftMargin(0.10);
      TempCanvas->SetTopMargin(0.06);
      TempCanvas->cd();
      BetaHist[i][j]->Draw();
      TempLine->Draw("same");
      Fitter->Draw("same");
      Legend->Draw("same");
      TempCanvas->Write();
      BetaHist[i][j]->Write();
 
      delete TempCanvas;
      delete Fitter;
    }
    DirBeta[i]->Write();
    delete DirBeta[i];
  }
  BetaDir->Write();
  delete BetaDir;
 
  gErrorIgnoreLevel = originalErrorLevel;
  Outputfile->cd();
}

PrepareOutput()

void SampleSummary::PrepareOutput ( )

inlineprivate

KS: Prepare output tree and necessary variables.

Definition at line 556 of file SampleSummary.cpp.

                                  {
// **********************
 
  // Make the output file (MakePosterioPredictive call writes to this)
  std::string TempString = OutputName;
  TempString.replace(TempString.find(".root"), 5, std::string("_procsW2.root"));
  Outputfile = new TFile(TempString.c_str(), "RECREATE");
 
  // The array of doubles we write to the TTree
  // Data vs Draw
  llh_data_draw.resize(nSamples);
  // Fluctuated Draw vs Draw
  llh_drawfluc_draw.resize(nSamples);
  // Fluctuated Predicitve vs Draw
  llh_predfluc_draw.resize(nSamples);
 
  // Data vs Draw using Rate
  llh_rate_data_draw.resize(nSamples);
  // Data vs Fluctuated Predictive using Rate
  llh_rate_predfluc_draw.resize(nSamples);
 
  // Data vs Fluctuated Draw
  llh_data_drawfluc.resize(nSamples);
  // Data vs Fluctuated Predictive
  llh_data_predfluc.resize(nSamples);
  // Draw vs Predictive
  llh_draw_pred.resize(nSamples);
  // Fluctuated Draw vs Predictive
  llh_drawfluc_pred.resize(nSamples);
  // Fluctuated Draw vs Fluctuated Predictive
  llh_drawfluc_predfluc.resize(nSamples);
 
  // Fluctuated Predictive vs Predictive
  llh_predfluc_pred.resize(nSamples);
  // Fluctuated Data vs Draw
  llh_datafluc_draw.resize(nSamples);
  
  // Data vs Draw for 1D projection
  llh_data_draw_ProjectX.resize(nSamples);
  llh_drawfluc_draw_ProjectX.resize(nSamples);
    
  // The output tree we're going to write to
  OutputTree = new TTree("LLH_draws", "LLH_draws");
  SampleNames.resize(nSamples);
  // Loop over the samples and set the addresses of the variables to write to file
  for (int i = 0; i < nSamples; ++i)
  {
    // Get the name
    std::string SampleName = SampleHandler->GetSampleName(i);
    // Strip out spaces
    while (SampleName.find(" ") != std::string::npos) {
      SampleName.replace(SampleName.find(" "), 1, std::string("_"));
    }
    SampleNames[i] = SampleName;
    //CW: Also strip out - signs because it messes up TBranches
    while (SampleName.find("-") != std::string::npos) {
      SampleName.replace(SampleName.find("-"), 1, std::string("_"));
    }
//  All LLH below are used for actual p-value calculations
    OutputTree->Branch((SampleName+"_data_draw").c_str(),     &llh_data_draw[i]);
    OutputTree->Branch((SampleName+"_drawfluc_draw").c_str(), &llh_drawfluc_draw[i]);
    OutputTree->Branch((SampleName+"_predfluc_draw").c_str(), &llh_predfluc_draw[i]);
 
//  All LLH below are used for actual p-value calculations however using rate only
    OutputTree->Branch((SampleName+"_rate_data_draw").c_str(), &llh_rate_data_draw[i]);
    OutputTree->Branch((SampleName+"_rate_predfluc_draw").c_str(), &llh_rate_predfluc_draw[i]);
 
//  All LLH below are for validation reason but not used for final P-Value
    OutputTree->Branch((SampleName+"_data_drawfluc").c_str(), &llh_data_drawfluc[i]);
    OutputTree->Branch((SampleName+"_data_predfluc").c_str(), &llh_data_predfluc[i]);
    OutputTree->Branch((SampleName+"_draw_pred").c_str(),     &llh_draw_pred[i]);
    OutputTree->Branch((SampleName+"_drawfluc_pred").c_str(), &llh_drawfluc_pred[i]);
    OutputTree->Branch((SampleName+"_drawfluc_predfluc").c_str(), &llh_drawfluc_predfluc[i]);
    OutputTree->Branch((SampleName+"_predfluc_pred").c_str(), &llh_predfluc_pred[i]);
    OutputTree->Branch((SampleName+"_datafluc_draw").c_str(), &llh_datafluc_draw[i]);
 
//  All LLH below are used for calcauting P-Value but using 1D projections
    OutputTree->Branch((SampleName+"_data_draw_ProjectX").c_str(), &llh_data_draw_ProjectX[i]);
    OutputTree->Branch((SampleName+"_drawfluc_draw_ProjectX").c_str(), &llh_drawfluc_draw_ProjectX[i]);
  }
//All LLH below are used for actual p-value calculations
  OutputTree->Branch("LLH_Penalty",         &llh_penalty);
  OutputTree->Branch("Total_LLH_Data_Draw", &total_llh_data_draw);
  OutputTree->Branch("Total_LLH_DrawFluc_Draw", &total_llh_drawfluc_draw);
  OutputTree->Branch("Total_LLH_PredFluc_Draw", &total_llh_predfluc_draw);
  
//  All LLH below are used for actual p-value calculations however using rate only
  OutputTree->Branch("Total_LLH_Rate_PredFluc_Draw", &total_llh_rate_predfluc_draw);
 
//All LLH below are for validation reason but not used for final P-Value
  OutputTree->Branch("Total_LLH_Data_DrawFluc", &total_llh_data_drawfluc);
  OutputTree->Branch("Total_LLH_Data_PredFluc", &total_llh_data_predfluc);
  OutputTree->Branch("Total_LLH_Draw_Pred",     &total_llh_draw_pred);
  OutputTree->Branch("Total_LLH_DrawFluc_Pred", &total_llh_drawfluc_pred);
  OutputTree->Branch("Total_LLH_DrawFluc_PredFluc", &total_llh_drawfluc_predfluc);
  OutputTree->Branch("Total_LLH_PredFluc_Pred", &total_llh_predfluc_pred);
  OutputTree->Branch("Total_LLH_DataFluc_Draw", &total_llh_datafluc_draw);
  
//All LLH below are used for calcauting P-Value but 1D projections
  OutputTree->Branch("total_llh_data_draw_ProjectX", &total_llh_data_draw_ProjectX);
  OutputTree->Branch("total_llh_drawfluc_draw_ProjectX", &total_llh_drawfluc_draw_ProjectX);
 
  Outputfile->cd();
  Dir.resize(nSamples);
  for (int i = 0; i < nSamples; ++i)
  {
    // Make a new directory
    Dir[i] = Outputfile->mkdir((SampleNames[i]).c_str());
  }
}

ProjectHist()

TH1D * SampleSummary::ProjectHist ( TH2D *  Histogram, const bool  ProjectX  )

inlineprivate

Helper to project TH2D onto axis.

Definition at line 2147 of file SampleSummary.cpp.

                                                                     {
// ****************
  TH1D* Projection = nullptr;
  std::string name;
  if (ProjectX) {
    name = std::string(Histogram->GetName()) + "_x";
    Projection = Histogram->ProjectionX(name.c_str(), 1, Histogram->GetYaxis()->GetNbins(), "e");
  } else {
    name = std::string(Histogram->GetName()) + "_y";
    Projection = Histogram->ProjectionY(name.c_str(), 1, Histogram->GetXaxis()->GetNbins(), "e");
  }
  return Projection;
}

ProjectPoly()

TH1D * SampleSummary::ProjectPoly ( TH2Poly *  Histogram, const bool  ProjectX, const int  selection, const bool  MakeErrorHist = false  )

inlineprivate

Helper to project TH2Poly onto axis.

Definition at line 2163 of file SampleSummary.cpp.

                                                                                                                       {
// ****************
  std::vector<double> xbins;
  std::vector<double> ybins;
 
  SampleHandler->SetupBinning(M3::int_t(selection), xbins, ybins);
  TH1D* Projection = nullptr;
  std::string name;
  if (ProjectX) {
    name = std::string(Histogram->GetName()) + "_x";
    Projection = PolyProjectionX(Histogram, name.c_str(), xbins, MakeErrorHist);
  }  else {
    name = std::string(Histogram->GetName()) + "_y";
    Projection = PolyProjectionY(Histogram, name.c_str(), ybins, MakeErrorHist);
  }
  return Projection;
}

SetLikelihood()

void SampleSummary::SetLikelihood ( const TestStatistic  TestStat )

inline

KS: Set likelihood type.

Definition at line 48 of file SampleSummary.h.

{ likelihood = TestStat;};

SetNModelParams()

void SampleSummary::SetNModelParams ( const int  nPars )

inline

Set number of model params used for BIC.

Definition at line 50 of file SampleSummary.h.

{ nModelParams = nPars;};

StudyBIC()

void SampleSummary::StudyBIC ( )

inlineprivate

Study Bayesian Information Criterion (BIC) [11].

Definition at line 2227 of file SampleSummary.cpp.

                            {
// ****************
  //make fancy event rate histogram
  double DataRate = 0.0;
  double BinsRate = 0.0;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:DataRate, BinsRate)
  #endif
  for (int i = 0; i < nSamples; ++i)
  {
    if (DataHist[i] == nullptr) continue;
    DataRate += NoOverflowIntegral(DataHist[i]);
    BinsRate += maxBins[i];
  }
 
  const double EventRateBIC = GetBIC(llh_total, DataRate, nModelParams);
  const double BinBasedBIC = GetBIC(llh_total, BinsRate, nModelParams);
  MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of events: {:.2f}", EventRateBIC);
  MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of bins: {:.2f}", BinBasedBIC);
  MACH3LOG_INFO("Additional info: nModelParams {} DataRate: {:.2f} BinsRate: {:.2f}", nModelParams, DataRate, BinsRate);
}

StudyDIC()

void SampleSummary::StudyDIC ( )

inlineprivate

KS: Get the Deviance Information Criterion (DIC) [29] [31].

Definition at line 2251 of file SampleSummary.cpp.

                             {
// ****************
  //The posterior mean of the deviance
  double Dbar = 0.;
 
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:Dbar)
  #endif
  for (unsigned int i = 0; i < nThrows; ++i)
  {
    double LLH_temp = 0.;
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      // Get -2*log-likelihood
      LLH_temp += GetLLH(DataHist[SampleNum], MCVector[i][SampleNum], W2MCVector[i][SampleNum]);
    }
    Dbar += LLH_temp;
  }
  Dbar = Dbar / nThrows;
 
  // A point estimate of the deviance
  const double Dhat = llh_total;
 
  //Effective number of parameters
  const double p_D = std::fabs(Dbar - Dhat);
 
  //Actual test stat
  const double DIC_stat = Dhat + 2 * p_D;
  MACH3LOG_INFO("Effective number of parameters following DIC formalism is equal to: {:.2f}", p_D);
  MACH3LOG_INFO("DIC test statistic = {:.2f}", DIC_stat);
}

StudyInformationCriterion()

void SampleSummary::StudyInformationCriterion ( M3::kInfCrit  Criterion )

inlineprivate

Information Criterion.

Definition at line 2199 of file SampleSummary.cpp.

                                                                  {
// ****************
  MACH3LOG_INFO("******************************");
  switch(Criterion) {
    case M3::kInfCrit::kBIC:
      // Study Bayesian Information Criterion
      StudyBIC();
      break;
    case M3::kInfCrit::kDIC:
      // Study Deviance Information Criterion
      StudyDIC();
      break;
    case M3::kInfCrit::kWAIC:
      // Study Watanabe-Akaike information criterion (WAIC)
      StudyWAIC();
      break;
    case M3::kInfCrit::kInfCrits:
      MACH3LOG_ERROR("kInfCrits is not a valid kInfCrit!");
      throw MaCh3Exception(__FILE__, __LINE__);
    default:
      MACH3LOG_ERROR("UNKNOWN Information Criterion SPECIFIED!");
      MACH3LOG_ERROR("You gave {}", static_cast<int>(Criterion));
      throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  MACH3LOG_INFO("******************************");
}

StudyKinematicCorrelations()

void SampleSummary::StudyKinematicCorrelations ( )

inlineprivate

KS: Study how correlated are sample or kinematic bins.

Definition at line 1690 of file SampleSummary.cpp.

                                               {
// ****************
  MACH3LOG_INFO("Calculating Correlations");
  TStopwatch timer;
  timer.Start();
 
    // Data vs Draw for 1D projection
  std::vector<double> NEvents_Sample(nSamples);
  double event_rate = 0.;
 
  // The output tree we're going to write to
  TTree* Event_Rate_Tree = new TTree("Event_Rate_draws", "Event_Rate_draws");
  Event_Rate_Tree->Branch("Event_Rate", &event_rate);
  // Loop over the samples and set the addresses of the variables to write to file
  for (int i = 0; i < nSamples; ++i)
  {
    // Get the name
    std::string SampleName = SampleNames[i];
    //CW: Also strip out - signs because it messes up TBranches
    while (SampleName.find("-") != std::string::npos) {
      SampleName.replace(SampleName.find("-"), 1, std::string("_"));
    }
    Event_Rate_Tree->Branch((SampleName+"_Event_Rate").c_str(), &NEvents_Sample[i]);
  }
 
  // Holds the total event rate
  auto EventHist = std::make_unique<TH1D>("EventHist", "Total Event Rate", 100, 1, -1);
  EventHist->SetDirectory(nullptr);
  EventHist->GetXaxis()->SetTitle("Total event rate");
  EventHist->GetYaxis()->SetTitle("Counts");
  EventHist->SetLineWidth(2);
 
  // Holds the event rate for the distribution
  std::vector<std::unique_ptr<TH1D>> SumHist(nSamples);
  for (int i = 0; i < nSamples; ++i)
  {
    std::string name = std::string(NominalHist[i]->GetName());
    name = name.substr(0, name.find("_nom"));
 
    SumHist[i] = std::make_unique<TH1D>((name+"_sum").c_str(),(name+"_sum").c_str(), 100, 1, -1);
    SumHist[i]->GetXaxis()->SetTitle("N_{events}");
    SumHist[i]->GetYaxis()->SetTitle("Counts");
    double Integral = NoOverflowIntegral(DataHist[i]);
    std::stringstream ss;
    ss << Integral;
    SumHist[i]->SetTitle((std::string(SumHist[i]->GetTitle())+"_"+ss.str()).c_str());
  }
 
  for (unsigned int it = 0; it < nThrows; ++it)
  {
    double event_rate_temp = 0.;
    // Loop over the samples
    #ifdef MULTITHREAD
    #pragma omp parallel for reduction(+:event_rate_temp)
    #endif
    for (int SampleNum = 0;  SampleNum < nSamples; ++SampleNum)
    {
      NEvents_Sample[SampleNum] = NoOverflowIntegral(MCVector[it][SampleNum]);
      // Fill the sum histogram with the integral of the sampled distribution
      SumHist[SampleNum]->Fill(NEvents_Sample[SampleNum], WeightVector[it]);
 
      event_rate_temp += NEvents_Sample[SampleNum];
    } // end samples loop
    event_rate = event_rate_temp;
    EventHist->Fill(event_rate);
    Event_Rate_Tree->Fill();
  } //end loops over throws
  Event_Rate_Tree->Write();
  delete Event_Rate_Tree;
 
  double DataRate = 0.0;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:DataRate)
  #endif
  for (int i = 0; i < nSamples; ++i)
  {
    DataRate += NoOverflowIntegral(DataHist[i]);
  }
  MakeCutEventRate(EventHist.get(), DataRate);
 
  for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
  {
    Dir[SampleNum]->cd();
    //Make fancy event rate histogram
    MakeCutEventRate(SumHist[SampleNum].get(), NoOverflowIntegral(DataHist[SampleNum]));
  }
 
  // Make a new directory
  TDirectory *CorrDir = Outputfile->mkdir("Correlations");
  CorrDir->cd();
 
  TMatrixDSym* SampleCorrelation = new TMatrixDSym(nSamples);
  std::vector<std::vector<std::unique_ptr<TH2D>>> SamCorr(nSamples);
  for (int i = 0; i < nSamples; ++i)
  {
    SamCorr[i].resize(nSamples);
 
    (*SampleCorrelation)(i,i) = 1.0;
    const double Min_i = SumHist[i]->GetXaxis()->GetBinLowEdge(1);
    const double Max_i = SumHist[i]->GetXaxis()->GetBinUpEdge(SumHist[i]->GetNbinsX()+1);
    for (int j = 0; j < nSamples; ++j)
    {
      const double Min_j = SumHist[j]->GetXaxis()->GetBinLowEdge(1);
      const double Max_j = SumHist[j]->GetXaxis()->GetBinUpEdge(SumHist[j]->GetNbinsX()+1);
 
      // TH2D to hold the Correlation
      SamCorr[i][j] = std::make_unique<TH2D>(Form("SamCorr_%i_%i", i, j), Form("SamCorr_%i_%i", i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
      SamCorr[i][j]->SetDirectory(nullptr);
      SamCorr[i][j]->SetMinimum(0);
      SamCorr[i][j]->GetXaxis()->SetTitle(SampleNames[i].c_str());
      SamCorr[i][j]->GetYaxis()->SetTitle(SampleNames[j].c_str());
      SamCorr[i][j]->GetZaxis()->SetTitle("Events");
    }
  }
 
  // Now we are sure we have the diagonal elements, let's make the off-diagonals
  #ifdef MULTITHREAD
  #pragma omp parallel for
  #endif
  for (int i = 0; i < nSamples; ++i)
  {
    for (int j = 0; j <= i; ++j)
    {
      // Skip the diagonal elements which we've already done above
      if (j == i) continue;
 
      for (unsigned int it = 0; it < nThrows; ++it)
      {
        SamCorr[i][j]->Fill(NoOverflowIntegral(MCVector[it][i]), NoOverflowIntegral(MCVector[it][j]));
      }
      SamCorr[i][j]->Smooth();
 
      // Get the Covariance for these two parameters
      (*SampleCorrelation)(i,j) = SamCorr[i][j]->GetCorrelationFactor();
      (*SampleCorrelation)(j,i) = (*SampleCorrelation)(i,j);
    }// End j loop
  }// End i loop
 
  auto hSamCorr = std::make_unique<TH2D>("Sample Correlation", "Sample Correlation", nSamples, 0, nSamples, nSamples, 0, nSamples);
  hSamCorr->SetDirectory(nullptr);
  hSamCorr->GetZaxis()->SetTitle("Correlation");
  hSamCorr->SetMinimum(-1);
  hSamCorr->SetMaximum(1);
  hSamCorr->GetXaxis()->SetLabelSize(0.015);
  hSamCorr->GetYaxis()->SetLabelSize(0.015);
 
  // Loop over the Covariance matrix entries
  for (int i = 0; i < nSamples; ++i)
  {
    hSamCorr->GetXaxis()->SetBinLabel(i+1, SampleNames[i].c_str());
 
    for (int j = 0; j < nSamples; ++j)
    {
      hSamCorr->GetYaxis()->SetBinLabel(j+1, SampleNames[j].c_str());
      // The value of the Covariance
      const double corr = (*SampleCorrelation)(i,j);
      hSamCorr->SetBinContent(i+1, j+1, corr);
    }
  }
  hSamCorr->Draw("colz");
  hSamCorr->Write("Sample_Corr");
 
  SampleCorrelation->Write("Sample_Correlation");
  delete SampleCorrelation;
 
  for (int i = 0; i < nSamples; ++i)
  {
    for (int j = 0; j <= i; ++j)
    {
      // Skip the diagonal elements which we've already done above
      if (j == i) continue;
      SamCorr[i][j]->Write();
    }// End j loop
  }// End i loop
 
  //KS: This can take ages so better turn it off by default
  bool DoPerKinemBin = false;
  if(DoPerKinemBin)
  {
    //KS: Now the same but for kinematic bin of each sample
    for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
    {
      TMatrixDSym* KinCorrelation = new TMatrixDSym(maxBins[SampleNum]);
      std::vector<std::vector<std::unique_ptr<TH2D>>> KinCorr(maxBins[SampleNum]);
      for (int i = 0; i < maxBins[SampleNum]; ++i)
      {
        KinCorr[i].resize(maxBins[SampleNum]);
        (*KinCorrelation)(i,i) = 1.0;
 
        const double Min_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinLowEdge(1);
        const double Max_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][i+1]->GetNbinsX()+1);
 
        //Get PolyBin
        TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
        // Just make a little fancy name
        std::stringstream ss2;
        ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
        ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
 
        for (int j = 0; j < maxBins[SampleNum]; ++j)
        {
          const double Min_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinLowEdge(1);
          const double Max_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][j+1]->GetNbinsX()+1);
 
          // TH2D to hold the Correlation
          KinCorr[i][j] = std::make_unique<TH2D>( Form("Kin_%i_%i_%i", SampleNum, i, j),
                    Form("Kin_%i_%i_%i", SampleNum, i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
          KinCorr[i][j]->SetDirectory(nullptr);
          KinCorr[i][j]->SetMinimum(0);
 
          KinCorr[i][j]->GetXaxis()->SetTitle(ss2.str().c_str());
 
          bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
          // Just make a little fancy name
          std::stringstream ss3;
          ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
          ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
          KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
          KinCorr[i][j]->GetZaxis()->SetTitle("Events");
        }
      }
      // Now we are sure we have the diagonal elements, let's make the off-diagonals
      #ifdef MULTITHREAD
      #pragma omp parallel for
      #endif
      for (int i = 0; i < maxBins[SampleNum]; ++i)
      {
        for (int j = 0; j <= i; ++j)
        {
          // Skip the diagonal elements which we've already done above
          if (j == i) continue;
 
          for (unsigned int it = 0; it < nThrows; ++it)
          {
            KinCorr[i][j]->Fill(MCVector[it][SampleNum]->GetBinContent(i+1), MCVector[it][SampleNum]->GetBinContent(j+1));
          }
          KinCorr[i][j]->Smooth();
 
          // Get the Covariance for these two parameters
          (*KinCorrelation)(i,j) = KinCorr[i][j]->GetCorrelationFactor();
          (*KinCorrelation)(j,i) = (*KinCorrelation)(i,j);
        }// End j loop
      }// End i loop
 
      auto hKinCorr = std::make_unique<TH2D>(SampleNames[SampleNum].c_str(), SampleNames[SampleNum].c_str(),
                                          maxBins[SampleNum], 0, maxBins[SampleNum], maxBins[SampleNum], 0, maxBins[SampleNum]);
      hKinCorr->SetDirectory(nullptr);
      hKinCorr->GetZaxis()->SetTitle("Correlation");
      hKinCorr->SetMinimum(-1);
      hKinCorr->SetMaximum(1);
      hKinCorr->GetXaxis()->SetLabelSize(0.015);
      hKinCorr->GetYaxis()->SetLabelSize(0.015);
 
      // Loop over the Covariance matrix entries
      for (int i = 0; i < maxBins[SampleNum]; ++i)
      {
        //Get PolyBin
        TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
        // Just make a little fancy name
        std::stringstream ss2;
        ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
        ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
        hKinCorr->GetXaxis()->SetBinLabel(i+1, ss2.str().c_str());
 
        for (int j = 0; j < maxBins[SampleNum]; ++j)
        {
          bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
          // Just make a little fancy name
          std::stringstream ss3;
          ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
          ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
          KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
 
          hKinCorr->GetYaxis()->SetBinLabel(j+1, ss3.str().c_str());
          // The value of the Covariance
          const double corr = (*KinCorrelation)(i,j);
          hKinCorr->SetBinContent(i+1, j+1, corr);
        }
      }
      hKinCorr->Draw("colz");
      hKinCorr->Write((SampleNames[SampleNum] + "_Corr").c_str());
 
      KinCorrelation->Write((SampleNames[SampleNum] + "_Correlation").c_str());
      delete KinCorrelation;
 
      /*
      for (int i = 0; i < maxBins[SampleNum]; ++i)
      {
        for (int j = 0; j <= i; ++j)
        {
          // Skip the diagonal elements which we've already done above
          if (j == i) continue;
          KinCorr[i][j]->Write();
        }// End j loop
      }// End i loop
      */
    }//end loop over samples
  }//end if DoPerKinemBin
  else
  {
    MACH3LOG_INFO("Not calculating correlations per each kinematic bin");
  }
 
  if(DoByModePlots)
  {
    // Holds the total event rate by mode
    std::vector<TH1D*> EventHist_ByMode(Modes->GetNModes()+1);
    for (int j = 0; j < Modes->GetNModes()+1; j++)
    {
      std::string ModeName = Modes->GetMaCh3ModeName(j);
      EventHist_ByMode[j] = new TH1D(Form("EventHist_%s", ModeName.c_str()), Form("Total Event Rate %s", ModeName.c_str()), 100, 1, -1);
      EventHist_ByMode[j]->GetXaxis()->SetTitle("Total event rate");
      EventHist_ByMode[j]->GetYaxis()->SetTitle("Counts");
      EventHist_ByMode[j]->SetLineWidth(2);
    }
 
    //KS: Here we calculate total event rates for each mode, maybe not most efficient but can be improved in the future
    for (unsigned int it = 0; it < nThrows; ++it)
    {
      for (int j = 0; j < Modes->GetNModes()+1; j++)
      {
        double event_rate_temp = 0.;
        #ifdef MULTITHREAD
        #pragma omp parallel for reduction(+:event_rate_temp)
        #endif
        for (int SampleNum = 0;  SampleNum < nSamples; SampleNum++)
        {
          event_rate_temp += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
        }
        EventHist_ByMode[j]->Fill(event_rate_temp);
      }
    }
 
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      MakeCutEventRate(EventHist_ByMode[i], DataRate);
    }
 
    TMatrixDSym* ModeCorrelation = new TMatrixDSym(Modes->GetNModes()+1);
 
    TH2D*** ModeCorr = new TH2D**[Modes->GetNModes()+1]();
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      ModeCorr[i] = new TH2D*[Modes->GetNModes()+1]();
 
      (*ModeCorrelation)(i,i) = 1.0;
 
      const double Min_i = EventHist_ByMode[i]->GetXaxis()->GetBinLowEdge(1);
      const double Max_i = EventHist_ByMode[i]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[i]->GetNbinsX()+1);
      for (int j = 0; j < Modes->GetNModes()+1; ++j)
      {
        const double Min_j = EventHist_ByMode[j]->GetXaxis()->GetBinLowEdge(1);
        const double Max_j = EventHist_ByMode[j]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[j]->GetNbinsX()+1);
 
        // TH2D to hold the Correlation
        ModeCorr[i][j] = new TH2D(Form("ModeCorr_%i_%i",i,j), Form("ModeCorr_%i_%i",i,j), 70, Min_i, Max_i, 70, Min_j, Max_j);
        ModeCorr[i][j]->SetDirectory(nullptr);
        ModeCorr[i][j]->SetMinimum(0);
        ModeCorr[i][j]->GetXaxis()->SetTitle(Modes->GetMaCh3ModeName(i).c_str());
        ModeCorr[i][j]->GetYaxis()->SetTitle(Modes->GetMaCh3ModeName(j).c_str());
        ModeCorr[i][j]->GetZaxis()->SetTitle("Events");
      }
    }
 
    // Now we are sure we have the diagonal elements, let's make the off-diagonals
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      for (int j = 0; j <= i; ++j)
      {
        // Skip the diagonal elements which we've already done above
        if (j == i) continue;
 
        for (unsigned int it = 0; it < nThrows; ++it)
        {
          double Integral_X = 0.;
          double Integral_Y = 0.;
          for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
          {
            Integral_X += NoOverflowIntegral(MCVectorByMode[it][SampleNum][i]);
            Integral_Y += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
          }
          ModeCorr[i][j]->Fill(Integral_X, Integral_Y);
        }
        ModeCorr[i][j]->Smooth();
 
        // Get the Covariance for these two parameters
        (*ModeCorrelation)(i,j) = ModeCorr[i][j]->GetCorrelationFactor();
        (*ModeCorrelation)(j,i) = (*ModeCorrelation)(i,j);
      }// End j loop
    }// End i loop
 
    TH2D* hModeCorr = new TH2D("Mode Correlation", "Mode Correlation", Modes->GetNModes()+1, 0, Modes->GetNModes()+1, Modes->GetNModes()+1, 0, Modes->GetNModes()+1);
    hModeCorr->SetDirectory(nullptr);
    hModeCorr->GetZaxis()->SetTitle("Correlation");
    hModeCorr->SetMinimum(-1);
    hModeCorr->SetMaximum(1);
    hModeCorr->GetXaxis()->SetLabelSize(0.015);
    hModeCorr->GetYaxis()->SetLabelSize(0.015);
 
    // Loop over the Covariance matrix entries
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      hModeCorr->GetXaxis()->SetBinLabel(i+1, Modes->GetMaCh3ModeName(i).c_str());
 
      for (int j = 0; j < Modes->GetNModes()+1; ++j)
      {
        hModeCorr->GetYaxis()->SetBinLabel(j+1, Modes->GetMaCh3ModeName(j).c_str());
        // The value of the Covariance
        const double corr = (*ModeCorrelation)(i,j);
        hModeCorr->SetBinContent(i+1, j+1, corr);
      }
    }
    hModeCorr->Draw("colz");
    hModeCorr->Write("Mode_Corr");
    delete hModeCorr;
 
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      for (int j = 0; j <= i; ++j)
      {
        // Skip the diagonal elements which we've already done above
        if (j == i) continue;
        ModeCorr[i][j]->Write();
      }// End j loop
    }// End i loop
 
    for (int i = 0; i < Modes->GetNModes()+1; ++i)
    {
      for (int j = 0; j < Modes->GetNModes()+1; ++j)
      {
        delete ModeCorr[i][j];
      }
      delete[] ModeCorr[i];
    }
    delete[] ModeCorr;
    ModeCorrelation->Write("Mode_Correlation");
    delete ModeCorrelation;
 
    for (int j = 0; j < Modes->GetNModes()+1; j++)
    {
      delete EventHist_ByMode[j];
    }
  }
 
  CorrDir->Close();
  delete CorrDir;
 
  timer.Stop();
  MACH3LOG_INFO("Calculating correlations took  {:.2f}s", timer.RealTime());
  Outputfile->cd();
}

StudyWAIC()

void SampleSummary::StudyWAIC ( )

inlineprivate

KS: Get the Watanabe-Akaike information criterion (WAIC) [11] [15].

Definition at line 2285 of file SampleSummary.cpp.

                              {
// ****************
  // log pointwise predictive density
  double lppd = 0.;
  // effective number of parameters
  double p_WAIC = 0.;
 
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:lppd, p_WAIC)
  #endif
  for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum) {
    int nBins = maxBins[SampleNum];
    for (int i = 1; i <= nBins; ++i) {
      double mean_llh = 0.;
      double sum_exp_llh = 0;
      double mean_llh_squared = 0.;
 
      for (unsigned int s = 0; s < nThrows; ++s) {
        const double data = DataHist[SampleNum]->GetBinContent(i);
        const double mc = MCVector[s][SampleNum]->GetBinContent(i);
        const double w2 = W2MCVector[s][SampleNum]->GetBinContent(i);
 
        // Get the -log-likelihood for this sample and bin
        double neg_LLH_temp = SampleHandler->GetTestStatLLH(data, mc, w2);
 
        // Negate the negative log-likelihood to get the actual log-likelihood
        double LLH_temp = -neg_LLH_temp;
 
        mean_llh += LLH_temp;
        mean_llh_squared += LLH_temp * LLH_temp;
        sum_exp_llh += std::exp(LLH_temp);
      }
 
      // Compute the mean log-likelihood and the squared mean
      mean_llh /= nThrows;
      mean_llh_squared /= nThrows;
      sum_exp_llh /= nThrows;
      sum_exp_llh = std::log(sum_exp_llh);
 
      // Log pointwise predictive density based on Eq. 4 in Gelman2014
      lppd += sum_exp_llh;
 
      // Compute the effective number of parameters for WAIC
      p_WAIC += mean_llh_squared - (mean_llh * mean_llh);
    }
  }
 
  // Compute WAIC, see Eq. 13 in Gelman2014
  double WAIC = -2 * (lppd - p_WAIC);
  MACH3LOG_INFO("Effective number of parameters following WAIC formalism is equal to: {:.2f}", p_WAIC);
  MACH3LOG_INFO("WAIC = {:.2f}", WAIC);
}

Write()

void SampleSummary::Write

(

)

KS: Write results into root file.

Definition at line 669 of file SampleSummary.cpp.

                          {
// *******************
  // Prepare the output tree
  PrepareOutput();
 
  MACH3LOG_INFO("Summarising {} throws...", nThrows);
  // After all the throws are added finalise the sample
  TStopwatch timer;
  timer.Start();
  MakePredictive();
  timer.Stop();
  MACH3LOG_INFO("Made Prior/Posterior Predictive, it took {:.2f}s, now writing...", timer.RealTime());
 
  // Studying information criterion
  StudyInformationCriterion(M3::kWAIC);
 
  OutputTree->Write();
 
  // Make the various distributions
  lnLHist->Write();
  lnLHist_drawfluc->Write();
  lnLHist_drawflucdraw->Write();
  lnLHist_drawdata->Write();
  lnLDrawHist->Write();
  lnLFlucHist->Write();
  lnLDrawHistRate->Write();
  //KS: Only available for Posterior Predictive
  if(!isPriorPredictive) RandomHist->Write();
 
  lnLFlucHist_ProjectX->Write();
  
  // Loop over each sample and write to file
  //KS: Multithreading is tempting here but we also write to ROOT file, separating all LLH and poly projections from write could work well
  for (int i = 0; i < nSamples; ++i)
  {
     // Skip the null histograms
    if (DataHist[i] == nullptr || NoOverflowIntegral(DataHist[i]) == 0) continue;
    Dir[i]->cd();
 
    // Make the data/MC ratio histogram
    TH2Poly *RatioHistMean = RatioPolys(DataHist[i], MeanHist[i]);
    RatioHistMean->GetZaxis()->SetTitle("Data/Mean");
    TH2Poly *RatioHistMode = RatioPolys(DataHist[i], ModeHist[i]);
    RatioHistMode->GetZaxis()->SetTitle("Data/Mode");
    TH2Poly *RatioHistNom = RatioPolys(DataHist[i], NominalHist[i]);
    RatioHistNom->GetZaxis()->SetTitle("Data/Nom");
 
    // And the normalised data histogram
    TH2Poly *DataNormHist = NormalisePoly(DataHist[i]);
    // Last true refers to if project along x or y
    TH2Poly *MeanNormHist = NormalisePoly(MeanHist[i]);
    TH2Poly *ModeNormHist = NormalisePoly(ModeHist[i]);
    TH1D *MeanProjectX = ProjectPoly(MeanHist[i], true, i, true);
    TH1D *MeanProjectY = ProjectPoly(MeanHist[i], false, i, true);
    TH1D *ModeProjectX = ProjectPoly(ModeHist[i], true, i, true);
    TH1D *ModeProjectY = ProjectPoly(ModeHist[i], false, i, true);
 
    TH1D *MeanHistCorrectedProjectX = nullptr;
    if(DoBetaParam) MeanHistCorrectedProjectX = ProjectPoly(MeanHistCorrected[i], true, i, true);
    TH1D *MeanHistCorrectedProjectY = nullptr;
    if(DoBetaParam) MeanHistCorrectedProjectY = ProjectPoly(MeanHistCorrected[i], false, i, true);
 
    TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[i], true, i);
    TH1D *W2MeanProjectY = ProjectPoly(W2MeanHist[i], false, i);
    TH1D *W2ModeProjectX = ProjectPoly(W2ModeHist[i], true, i);
    TH1D *W2ModeProjectY = ProjectPoly(W2ModeHist[i], false, i);
 
    TH2Poly *NomNormHist = NormalisePoly(NominalHist[i]);
    TH1D *NomProjectX = ProjectPoly(NominalHist[i], true, i);
    TH1D *NomProjectY = ProjectPoly(NominalHist[i], false, i);
 
    TH1D *W2NomProjectX = ProjectPoly(W2NomHist[i], true, i);
    TH1D *W2NomProjectY = ProjectPoly(W2NomHist[i], false, i);
 
    // Same for the TH2Ds
    CalcLLH(DataHist[i], NominalHist[i], W2NomHist[i]);
    CalcLLH(DataHist[i], MeanHist[i], W2MeanHist[i]);
    CalcLLH(DataHist[i], ModeHist[i], W2ModeHist[i]);
 
    // Calculate the log likelihood for the 1D dists
    // Sets the title of the second TH1D to the -2LLH
    CalcLLH(DataHist_ProjectX[i], NomProjectX, W2NomProjectX);
    CalcLLH(DataHist_ProjectX[i], MeanProjectX, W2MeanProjectX);
    CalcLLH(DataHist_ProjectX[i], ModeProjectX, W2ModeProjectX);
    CalcLLH(DataHist_ProjectY[i], NomProjectY, W2NomProjectY);
    CalcLLH(DataHist_ProjectY[i], MeanProjectY, W2MeanProjectY);
    CalcLLH(DataHist_ProjectY[i], ModeProjectY, W2ModeProjectY);
 
    std::string SampleName = SampleNames[i];
    // Also strip out - signs because it messes up TBranches
    while (SampleName.find("-") != std::string::npos) {
      SampleName.replace(SampleName.find("-"), 1, std::string("_"));
    } 
    OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_drawfluc_draw>>htemp").c_str());
    TH2D *TempHistogram = static_cast<TH2D*>(gDirectory->Get("htemp")->Clone());
    TempHistogram->GetXaxis()->SetTitle("-2LLH(Draw Fluc, Draw)");
    TempHistogram->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
    TempHistogram->SetNameTitle((SampleNames[i]+"_drawfluc_draw").c_str(), (SampleNames[i]+"_drawfluc_draw").c_str());
    Get2DBayesianpValue(TempHistogram);
    TempHistogram->Write();
    delete TempHistogram;
 
    // Also write the 2D histograms for the p-value
    OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_predfluc_draw>>htemp2").c_str());
    TH2D *TempHistogram2 = static_cast<TH2D*>(gDirectory->Get("htemp2")->Clone());
    TempHistogram2->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
    TempHistogram2->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
    TempHistogram2->SetNameTitle((SampleNames[i]+"_predfluc_draw").c_str(), (SampleNames[i]+"_predfluc_draw").c_str());
    Get2DBayesianpValue(TempHistogram2);
    TempHistogram2->Write();
    delete TempHistogram2;
   
    // finally p-value for 1D projection
    OutputTree->Draw((SampleName+"_rate_data_draw:"+SampleName+"_rate_predfluc_draw>>htemp3").c_str());
    TH2D *TempHistogram3 = static_cast<TH2D*>(gDirectory->Get("htemp3")->Clone());
    TempHistogram3->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
    TempHistogram3->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
    TempHistogram3->SetNameTitle((SampleNames[i]+"_rate_predfluc_draw").c_str(), (SampleNames[i]+"_rate_predfluc_draw").c_str());
    Get2DBayesianpValue(TempHistogram3);
    TempHistogram3->Write();
    delete TempHistogram3;
 
    // finally p-value for 1D projection
    OutputTree->Draw((SampleName+"_data_draw_ProjectX:"+SampleName+"_drawfluc_draw_ProjectX>>htemp4").c_str());
    TH2D *TempHistogram4 = static_cast<TH2D*>(gDirectory->Get("htemp4")->Clone());
    TempHistogram4->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
    TempHistogram4->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
    TempHistogram4->SetNameTitle((SampleNames[i]+"_drawfluc_draw_ProjectX").c_str(), (SampleNames[i]+"_drawfluc_draw_ProjectX").c_str());
    Get2DBayesianpValue(TempHistogram4);
    TempHistogram4->Write();
    delete TempHistogram4;
    
    // Write the Histograms to each folder
    DataHist[i]->Write();
    NominalHist[i]->Write();
    MeanHist[i]->Write();
    ModeHist[i]->Write();
    RatioHistMean->Write();
    RatioHistMode->Write();
    RatioHistNom->Write();
    if(DoBetaParam) MeanHistCorrected[i]->Write();
 
    W2NomHist[i]->Write();
    W2MeanHist[i]->Write();
    W2ModeHist[i]->Write();
 
    DataNormHist->Write();
    NomNormHist->Write();
    MeanNormHist->Write();
    ModeNormHist->Write();
 
    DataHist_ProjectX[i]->Write();
    NomProjectX->Write();
    MeanProjectX->Write();
    ModeProjectX->Write();
    if(DoBetaParam) MeanHistCorrectedProjectX->Write();
    ViolinHists_ProjectX[i]->Write();
    
    DataHist_ProjectY[i]->Write();
    NomProjectY->Write();
    MeanProjectY->Write();
    ModeProjectY->Write();
    if(DoBetaParam) MeanHistCorrectedProjectY->Write();
    ViolinHists_ProjectY[i]->Write();
 
    W2NomProjectX->Write();
    W2MeanProjectX->Write();
    W2ModeProjectX->Write();
 
    W2NomProjectY->Write();
    W2MeanProjectY->Write();
    W2ModeProjectY->Write();
    
    //KS: This will dump lots of hists, use it only for debugging
    if(Debug > 0)
    {
      TDirectory* DebugDir = Dir[i]->mkdir("Debug");
      DebugDir->cd();
      for (int b = 1; b <= maxBins[i]; ++b)
      {
        PosteriorHist[i][b]->Write();
        std::string Title = PosteriorHist[i][b]->GetName();
 
        auto TempLine = std::make_unique<TLine>(NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
                                                NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
        TempLine->SetLineColor(kRed);
        TempLine->SetLineWidth(2);
 
        auto TempLineData = std::make_unique<TLine>(DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
                                                    DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
        TempLineData->SetLineColor(kGreen);
        TempLineData->SetLineWidth(2);
 
        // Also fit a Gaussian because why not?
        TF1 *Fitter = new TF1("Fit", "gaus", PosteriorHist[i][b]->GetBinLowEdge(1), PosteriorHist[i][b]->GetBinLowEdge(PosteriorHist[i][b]->GetNbinsX()+1));
        PosteriorHist[i][b]->Fit(Fitter, "RQ");
        Fitter->SetLineColor(kRed-5);
 
        auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
        Legend->SetFillColor(0);
        Legend->SetFillStyle(0);
        Legend->SetLineWidth(0);
        Legend->SetLineColor(0);
        Legend->AddEntry(TempLineData.get(), Form("Data #mu=%.2f", DataHist[i]->GetBinContent(b)), "l");
        Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.2f", NominalHist[i]->GetBinContent(b)), "l");
        Legend->AddEntry(PosteriorHist[i][b].get(), Form("Post, #mu=%.2f#pm%.2f", PosteriorHist[i][b]->GetMean(), PosteriorHist[i][b]->GetRMS()), "l");
        Legend->AddEntry(Fitter, Form("Gauss, #mu=%.2f#pm%.2f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
        std::string TempTitle = std::string(PosteriorHist[i][b]->GetName());
 
        TempTitle += "_canv";
        TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
        TempCanvas->SetGridx();
        TempCanvas->SetGridy();
        TempCanvas->SetRightMargin(0.03);
        TempCanvas->SetBottomMargin(0.08);
        TempCanvas->SetLeftMargin(0.10);
        TempCanvas->SetTopMargin(0.06);
        TempCanvas->cd();
        PosteriorHist[i][b]->Draw();
        TempLine->Draw("same");
        TempLineData->Draw("same");
        Fitter->Draw("same");
        Legend->Draw("same");
        TempCanvas->Write();
 
        delete TempCanvas;
        delete Fitter;
        //This isn't useful check only in desperation
        if(Debug > 1) w2Hist[i][b]->Write();
      }
      DebugDir->Close();
      delete DebugDir;
      Dir[i]->cd();
    }
    lnLHist_Mean[i]->Write();
    lnLHist_Mode[i]->Write();
 
    lnLHist_Mean_ProjectX[i]->Write();
 
    lnLHist_Mean1D[i]->Write();
    lnLHist_Mode1D[i]->Write();
 
    MakeCutLLH1D(lnLHist_Sample_DrawData[i], GetLLH(DataHist[i], MeanHist[i], W2MeanHist[i]));
    lnLHist_Sample_DrawData[i]->Write();
    MakeCutLLH1D(lnLHist_Sample_DrawflucDraw[i], GetLLH(DataHist[i], MeanHist[i], W2MeanHist[i]));
    lnLHist_Sample_DrawflucDraw[i]->Write();
    MakeCutLLH1D(lnLHist_Sample_PredflucDraw[i], GetLLH(DataHist[i], MeanHist[i], W2MeanHist[i]));
    lnLHist_Sample_PredflucDraw[i]->Write();
    
    if(DoByModePlots)
    {
      for (int j = 0; j < Modes->GetNModes()+1; ++j)
      {
        MeanHist_ByMode[i][j]->Write();
        TH1D *MeanProjectX_ByMode = ProjectPoly(MeanHist_ByMode[i][j], true, i, true);
        TH1D *MeanProjectY_ByMode = ProjectPoly(MeanHist_ByMode[i][j], false, i, true);
        MeanProjectX_ByMode->Write();
        MeanProjectY_ByMode->Write();
        //KS: This will dump lots of hists, use it only for debugging
        if(Debug > 0)
        {
          for (int b = 1; b <= maxBins[i]; ++b)
          {
            PosteriorHist_ByMode[i][j][b]->Write();
          }
        }
        delete MeanProjectX_ByMode;
        delete MeanProjectY_ByMode;
      } // End loop over bins
    }
    // Delete temporary objects
    delete RatioHistMean;
    delete RatioHistMode;
    delete RatioHistNom;
 
    delete DataNormHist;
    delete MeanNormHist;
    delete ModeNormHist;
    delete NomNormHist;
 
    delete DataHist_ProjectX[i];
    delete MeanProjectX;
    delete ModeProjectX;
    if(DoBetaParam) delete MeanHistCorrectedProjectX;
    delete NomProjectX;
 
    delete DataHist_ProjectY[i];
    delete MeanProjectY;
    delete ModeProjectY;
    if(DoBetaParam) delete MeanHistCorrectedProjectY;
    delete NomProjectY;
 
    delete W2NomProjectX;
    delete W2MeanProjectX;
    delete W2ModeProjectX;
  
    delete W2NomProjectY;
    delete W2MeanProjectY;
    delete W2ModeProjectY;
    MACH3LOG_INFO("");
  } //end loop over samples
  if(DoBetaParam) PlotBetaParameters();
 
  StudyKinematicCorrelations();
  MACH3LOG_INFO("Wrote to {}", Outputfile->GetName());
}

Member Data Documentation

BetaHist

std::vector<std::vector<std::unique_ptr<TH1D> > > SampleSummary::BetaHist

private

Distribution of beta parameters in Barlow Beeston formalisms.

Definition at line 230 of file SampleSummary.h.

DataHist

std::vector<TH2Poly*> SampleSummary::DataHist

private

The data histogram for the selection.

Definition at line 167 of file SampleSummary.h.

DataHist_ProjectX

std::vector<TH1D*> SampleSummary::DataHist_ProjectX

private

The data histogram for the selection X projection.

Definition at line 169 of file SampleSummary.h.

DataHist_ProjectY

std::vector<TH1D*> SampleSummary::DataHist_ProjectY

private

The data histogram for the selection Y projection.

Definition at line 171 of file SampleSummary.h.

Debug

int SampleSummary::Debug

private

Tells Debug level to save additional histograms.

Definition at line 348 of file SampleSummary.h.

Dir

std::vector<TDirectory*> SampleSummary::Dir

private

Directory for each sample.

Definition at line 257 of file SampleSummary.h.

DoBetaParam

bool SampleSummary::DoBetaParam

private

Are we making Beta Histograms.

Definition at line 232 of file SampleSummary.h.

DoByModePlots

bool SampleSummary::DoByModePlots

private

By mode variables.

Definition at line 329 of file SampleSummary.h.

doShapeOnly

bool SampleSummary::doShapeOnly

private

bool whether to normalise each toy to have shape based p-value and pos pred distribution

Definition at line 241 of file SampleSummary.h.

first_pass

bool SampleSummary::first_pass

private

KS: Hacky flag to let us know if this is first toy.

Definition at line 133 of file SampleSummary.h.

isPriorPredictive

bool SampleSummary::isPriorPredictive

private

bool whether we have Prior or Posterior Predictive

Definition at line 238 of file SampleSummary.h.

likelihood

TestStatistic SampleSummary::likelihood

private

Type of likelihood for example Poisson, Barlow-Beeston or Ice Cube.

Definition at line 342 of file SampleSummary.h.

llh_data_draw

std::vector<double> SampleSummary::llh_data_draw

private

Data vs Draw.

Definition at line 262 of file SampleSummary.h.

llh_data_draw_ProjectX

std::vector<double> SampleSummary::llh_data_draw_ProjectX

private

Projection X (most likely muon momentum) of LLH.

Definition at line 290 of file SampleSummary.h.

llh_data_drawfluc

std::vector<double> SampleSummary::llh_data_drawfluc

private

Data vs Fluctuated Draw.

Definition at line 274 of file SampleSummary.h.

llh_data_predfluc

std::vector<double> SampleSummary::llh_data_predfluc

private

Data vs Fluctuated Predictive.

Definition at line 276 of file SampleSummary.h.

llh_datafluc_draw

std::vector<double> SampleSummary::llh_datafluc_draw

private

Fluctuated Data vs Draw.

Definition at line 287 of file SampleSummary.h.

llh_draw_pred

std::vector<double> SampleSummary::llh_draw_pred

private

Draw vs Predictive.

Definition at line 278 of file SampleSummary.h.

llh_drawfluc_draw

std::vector<double> SampleSummary::llh_drawfluc_draw

private

Fluctuated Draw vs Draw.

Definition at line 264 of file SampleSummary.h.

llh_drawfluc_draw_ProjectX

std::vector<double> SampleSummary::llh_drawfluc_draw_ProjectX

private

Definition at line 291 of file SampleSummary.h.

llh_drawfluc_pred

std::vector<double> SampleSummary::llh_drawfluc_pred

private

Fluctuated Draw vs Predictive.

Definition at line 280 of file SampleSummary.h.

llh_drawfluc_predfluc

std::vector<double> SampleSummary::llh_drawfluc_predfluc

private

Fluctuated Draw vs Fluctuated Predictive.

Definition at line 285 of file SampleSummary.h.

llh_penalty

double SampleSummary::llh_penalty

private

LLH penalty for each throw.

Definition at line 294 of file SampleSummary.h.

llh_predfluc_draw

std::vector<double> SampleSummary::llh_predfluc_draw

private

Fluctuated Predictive vs Draw.

Definition at line 266 of file SampleSummary.h.

llh_predfluc_pred

std::vector<double> SampleSummary::llh_predfluc_pred

private

Fluctuated Predictive vs Predictive.

Definition at line 283 of file SampleSummary.h.

llh_rate_data_draw

std::vector<double> SampleSummary::llh_rate_data_draw

private

Data vs Draw using rate only.

Definition at line 269 of file SampleSummary.h.

llh_rate_predfluc_draw

std::vector<double> SampleSummary::llh_rate_predfluc_draw

private

Fluctuated Predictive vs Draw using rate only.

Definition at line 271 of file SampleSummary.h.

llh_total

double SampleSummary::llh_total

private

Total LLH for the posterior predictive distribution.

Definition at line 250 of file SampleSummary.h.

LLHPenaltyVector

std::vector<double> SampleSummary::LLHPenaltyVector

private

Vector to hold the penalty term.

Definition at line 146 of file SampleSummary.h.

lnLDrawHist

std::unique_ptr<TH2D> SampleSummary::lnLDrawHist

private

The 2D lnLhist, showing (draw vs data) and (draw vs fluct), anything above y=x axis is the p-value.

Definition at line 190 of file SampleSummary.h.

lnLDrawHistRate

std::unique_ptr<TH2D> SampleSummary::lnLDrawHistRate

private

The 2D lnLhist, showing (draw vs data) and (draw vs fluct), using rate, anything above y=x axis is the p-value.

Definition at line 195 of file SampleSummary.h.

lnLFlucHist

std::unique_ptr<TH2D> SampleSummary::lnLFlucHist

private

The 2D lnLHist, showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-value.

Definition at line 192 of file SampleSummary.h.

lnLFlucHist_ProjectX

std::unique_ptr<TH2D> SampleSummary::lnLFlucHist_ProjectX

private

The 2D lnLHist but for ProjectionX histogram (pmu), showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-value.

Definition at line 197 of file SampleSummary.h.

lnLHist

std::unique_ptr<TH1D> SampleSummary::lnLHist

private

The histogram containing the lnL for each throw.

Definition at line 182 of file SampleSummary.h.

lnLHist_drawdata

std::unique_ptr<TH1D> SampleSummary::lnLHist_drawdata

private

The lnLhist for the draw vs data.

Definition at line 188 of file SampleSummary.h.

lnLHist_drawfluc

std::unique_ptr<TH1D> SampleSummary::lnLHist_drawfluc

private

The lnLhist for the draw vs MC fluctuated.

Definition at line 184 of file SampleSummary.h.

lnLHist_drawflucdraw

std::unique_ptr<TH1D> SampleSummary::lnLHist_drawflucdraw

private

The lnLhist for the draw vs draw fluctuated.

Definition at line 186 of file SampleSummary.h.

lnLHist_Mean

std::vector<TH2Poly*> SampleSummary::lnLHist_Mean

private

The LLH distribution in pmu cosmu for using the mean in each bin.

Definition at line 207 of file SampleSummary.h.

lnLHist_Mean1D

std::vector<TH1D*> SampleSummary::lnLHist_Mean1D

private

Holds the bin-by-bin LLH for the mean posterior predictive vs the data.

Definition at line 222 of file SampleSummary.h.

lnLHist_Mean_ProjectX

std::vector<TH1D*> SampleSummary::lnLHist_Mean_ProjectX

private

The LLH distribution in pmu using the mean in each bin.

Definition at line 212 of file SampleSummary.h.

lnLHist_Mode

std::vector<TH2Poly*> SampleSummary::lnLHist_Mode

private

The LLH distribution in pmu cosmu for using the mode in each bin.

Definition at line 209 of file SampleSummary.h.

lnLHist_Mode1D

std::vector<TH1D*> SampleSummary::lnLHist_Mode1D

private

Holds the bin-by-bin LLH for the mode posterior predictive vs the data.

Definition at line 224 of file SampleSummary.h.

lnLHist_Sample_DrawData

std::vector<TH1D*> SampleSummary::lnLHist_Sample_DrawData

private

The histogram containing the lnL (draw vs data) for each throw for each sample.

Definition at line 200 of file SampleSummary.h.

lnLHist_Sample_DrawflucDraw

std::vector<TH1D*> SampleSummary::lnLHist_Sample_DrawflucDraw

private

The histogram containing the lnL (draw vs draw fluct) for each throw for each sample.

Definition at line 202 of file SampleSummary.h.

lnLHist_Sample_PredflucDraw

std::vector<TH1D*> SampleSummary::lnLHist_Sample_PredflucDraw

private

The histogram containing the lnL (draw vs pred fluct) for each throw for each sample.

Definition at line 204 of file SampleSummary.h.

maxBins

std::vector<int> SampleSummary::maxBins

private

Max Number of Bins per each sample.

Definition at line 247 of file SampleSummary.h.

MCVector

std::vector<std::vector<TH2Poly*> > SampleSummary::MCVector

private

Vector of vectors which holds the loaded MC histograms.

Definition at line 139 of file SampleSummary.h.

MCVectorByMode

std::vector<std::vector<std::vector<TH2Poly*> > > SampleSummary::MCVectorByMode

private

Vector of vectors which holds the loaded MC histograms for each mode.

Definition at line 143 of file SampleSummary.h.

MeanHist

std::vector<TH2Poly*> SampleSummary::MeanHist

private

The posterior predictive distribution in pmu cosmu using the mean.

Definition at line 215 of file SampleSummary.h.

MeanHist_ByMode

std::vector<std::vector<TH2Poly*> > SampleSummary::MeanHist_ByMode

private

The posterior predictive distribution in pmu cosmu using the mean.

Definition at line 331 of file SampleSummary.h.

MeanHistCorrected

std::vector<TH2Poly*> SampleSummary::MeanHistCorrected

private

The posterior predictive distribution in pmu cosmu using the mean after applying Barlow-Beeston Correction.

Definition at line 217 of file SampleSummary.h.

ModeHist

std::vector<TH2Poly*> SampleSummary::ModeHist

private

The posterior predictive distribution in pmu cosmu using the mode.

Definition at line 219 of file SampleSummary.h.

Modes

MaCh3Modes* SampleSummary::Modes

private

MaCh3 Modes.

Definition at line 339 of file SampleSummary.h.

nChainSteps

unsigned int SampleSummary::nChainSteps

private

Number of throws by user.

Definition at line 235 of file SampleSummary.h.

nModelParams

int SampleSummary::nModelParams

private

Number of parameters.

Definition at line 345 of file SampleSummary.h.

NominalHist

std::vector<TH2Poly*> SampleSummary::NominalHist

private

The nominal histogram for the selection.

Definition at line 173 of file SampleSummary.h.

nSamples

int SampleSummary::nSamples

private

Number of samples.

Definition at line 151 of file SampleSummary.h.

nThrows

unsigned int SampleSummary::nThrows

private

Number of throws.

Definition at line 244 of file SampleSummary.h.

Outputfile

TFile* SampleSummary::Outputfile

private

Output filename.

Definition at line 255 of file SampleSummary.h.

OutputName

std::string SampleSummary::OutputName

private

Output filename.

Definition at line 253 of file SampleSummary.h.

OutputTree

TTree* SampleSummary::OutputTree

private

TTree which we save useful information to.

Definition at line 260 of file SampleSummary.h.

PosteriorHist

std::vector<std::vector<std::unique_ptr<TH1D> > > SampleSummary::PosteriorHist

private

The posterior predictive for the whole selection: this gets built after adding in the toys. Now an array of Th1ds, 1 for each poly bin, for each sample.

Definition at line 157 of file SampleSummary.h.

PosteriorHist_ByMode

TH1D**** SampleSummary::PosteriorHist_ByMode

private

Histogram which corresponds to each bin in the sample's th2poly.

Definition at line 333 of file SampleSummary.h.

RandomHist

std::unique_ptr<TH1D> SampleSummary::RandomHist

private

Holds the history of which entries have been drawn in the MCMC file.

Definition at line 227 of file SampleSummary.h.

rnd

std::unique_ptr<TRandom3> SampleSummary::rnd

private

Random number generator.

Definition at line 131 of file SampleSummary.h.

SampleHandler

SampleHandlerBase* SampleSummary::SampleHandler

private

Pointer to SampleHandler object, mostly used to get sample names, binning etc.

Definition at line 336 of file SampleSummary.h.

SampleNames

std::vector<std::string> SampleSummary::SampleNames

private

name for each sample

Definition at line 154 of file SampleSummary.h.

StandardFluctuation

bool SampleSummary::StandardFluctuation

private

KS: We have two methods for Poissonian fluctuation.

Definition at line 136 of file SampleSummary.h.

total_llh_data_draw

double SampleSummary::total_llh_data_draw

private

Data vs Draw.

Definition at line 297 of file SampleSummary.h.

total_llh_data_draw_ProjectX

double SampleSummary::total_llh_data_draw_ProjectX

private

Data vs Draw for projection X (most likely muon momentum)

Definition at line 324 of file SampleSummary.h.

total_llh_data_drawfluc

double SampleSummary::total_llh_data_drawfluc

private

Data vs Fluctuated Draw.

Definition at line 311 of file SampleSummary.h.

total_llh_data_predfluc

double SampleSummary::total_llh_data_predfluc

private

Data vs Fluctuated Predictive.

Definition at line 309 of file SampleSummary.h.

total_llh_datafluc_draw

double SampleSummary::total_llh_datafluc_draw

private

Fluctuated Data vs Draw.

Definition at line 319 of file SampleSummary.h.

total_llh_draw_pred

double SampleSummary::total_llh_draw_pred

private

Draw vs Predictive.

Definition at line 313 of file SampleSummary.h.

total_llh_drawfluc_draw

double SampleSummary::total_llh_drawfluc_draw

private

Fluctuated Draw vs Draw.

Definition at line 299 of file SampleSummary.h.

total_llh_drawfluc_draw_ProjectX

double SampleSummary::total_llh_drawfluc_draw_ProjectX

private

Fluctuated Draw vs Draw for projection X (most likely muon momentum)

Definition at line 326 of file SampleSummary.h.

total_llh_drawfluc_pred

double SampleSummary::total_llh_drawfluc_pred

private

Fluctuated Draw vs Predictive.

Definition at line 315 of file SampleSummary.h.

total_llh_drawfluc_predfluc

double SampleSummary::total_llh_drawfluc_predfluc

private

Fluctuated Draw vs Fluctuated Predictive.

Definition at line 317 of file SampleSummary.h.

total_llh_predfluc_draw

double SampleSummary::total_llh_predfluc_draw

private

Fluctuated Predictive vs Draw.

Definition at line 301 of file SampleSummary.h.

total_llh_predfluc_pred

double SampleSummary::total_llh_predfluc_pred

private

Fluctuated Predictive vs Predictive.

Definition at line 321 of file SampleSummary.h.

total_llh_rate_data_draw

double SampleSummary::total_llh_rate_data_draw

private

Rate Data vs Draw.

Definition at line 304 of file SampleSummary.h.

total_llh_rate_predfluc_draw

double SampleSummary::total_llh_rate_predfluc_draw

private

Fluctuated Predictive vs Draw using Rate.

Definition at line 306 of file SampleSummary.h.

ViolinHists_ProjectX

std::vector<TH2D*> SampleSummary::ViolinHists_ProjectX

private

Posterior predictive but for X projection but as a violin plot.

Definition at line 162 of file SampleSummary.h.

ViolinHists_ProjectY

std::vector<TH2D*> SampleSummary::ViolinHists_ProjectY

private

Posterior predictive but for Y projection but as a violin plot.

Definition at line 164 of file SampleSummary.h.

w2Hist

std::vector<std::vector<std::unique_ptr<TH1D> > > SampleSummary::w2Hist

private

The posterior predictive for the whole selection: this gets built after adding in the toys. Now an array of Th1ds, 1 for each poly bin, for each sample for W2.

Definition at line 159 of file SampleSummary.h.

W2MCVector

std::vector<std::vector<TH2Poly*> > SampleSummary::W2MCVector

private

Vector of vectors which holds the loaded W2 histograms.

Definition at line 141 of file SampleSummary.h.

W2MeanHist

std::vector<TH2Poly*> SampleSummary::W2MeanHist

private

Pointer to the w2 histograms (for mean values).

Definition at line 177 of file SampleSummary.h.

W2ModeHist

std::vector<TH2Poly*> SampleSummary::W2ModeHist

private

Pointer to the w2 histograms (for mode values).

Definition at line 179 of file SampleSummary.h.

W2NomHist

std::vector<TH2Poly*> SampleSummary::W2NomHist

private

Pointer to the w2 histograms (for nominal values).

Definition at line 175 of file SampleSummary.h.

WeightVector

std::vector<double> SampleSummary::WeightVector

private

Vector holding weight.

Definition at line 148 of file SampleSummary.h.

---

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

• Fitters/SampleSummary.h
• Fitters/SampleSummary.cpp