MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
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:
[legend]

Public Member Functions

 SampleSummary (const int n_Samples, const std::string &Filename, SampleHandlerBase *const sample, const int nSteps)
 Constructor.
 
 ~SampleSummary ()
 Destructor.
 
void AddData (std::vector< TH2Poly * > &DataHist)
 KS: Add data histograms.
 
void AddNominal (std::vector< TH2Poly * > &NominalHist, std::vector< TH2Poly * > &W2Nom)
 KS: Add prior histograms.
 
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.
 
void AddThrowByMode (std::vector< std::vector< TH2Poly * > > &SampleVector_ByMode)
 KS: Add histograms for each mode.
 
void Write ()
 KS: Write results into root file.
 
void SetLikelihood (const TestStatistic TestStat)
 KS: Set likelihood type.
 
void SetNModelParams (const int nPars)
 Set number of model params used for BIC.
 

Private Member Functions

void MakePredictive ()
 Finalise the distributions from the thrown samples.
 
void PrepareOutput ()
 KS: Prepare output tree and necessary variables.
 
void CalcLLH (TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
 Helper functions to calculate likelihoods using TH2Poly, will modify MC hist tittle to include LLH.
 
void CalcLLH (TH1D *const &Data, TH1D *const &MC, TH1D *const &W2)
 Helper functions to calculate likelihoods using TH1D, will modify MC hist tittle to include LLH.
 
double GetLLH (TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
 Helper functions to calculate likelihoods using TH2Poly.
 
double GetLLH (TH1D *const &Data, TH1D *const &MC, TH1D *const &W2)
 Helper functions to calculate likelihoods using TH1D.
 
void PlotBetaParameters ()
 KS: In Barlow Beeston we have Beta Parameters which scale generated MC.
 
void StudyKinematicCorrelations ()
 KS: Study how correlated are sample or kinematic bins.
 
void MakeCutLLH ()
 Make the cut LLH histogram.
 
void MakeCutLLH1D (TH1D *Histogram, double llh_ref=-999)
 
void MakeCutEventRate (TH1D *Histogram, const double DataRate)
 Make the 1D Event Rate Hist.
 
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.
 
bool CheckSamples (const int Length)
 Check the length of samples agrees.
 
TH1D * ProjectHist (TH2D *Histogram, const bool ProjectX)
 Helper to project TH2D onto axis.
 
TH1D * ProjectPoly (TH2Poly *Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist=false)
 Helper to project TH2Poly onto axis.
 
void MakeFluctuatedHistogram (TH1D *FluctHist, TH1D *PolyHist)
 Make Poisson fluctuation of TH1D hist.
 
void MakeFluctuatedHistogram (TH2Poly *FluctHist, TH2Poly *PolyHist)
 Make Poisson fluctuation of TH2Poly hist.
 
void StudyInformationCriterion (M3::kInfCrit Criterion)
 Information Criterion.
 
void StudyBIC ()
 Study Bayesian Information Criterion (BIC) [11].
 
void StudyDIC ()
 KS: Get the Deviance Information Criterion (DIC) [26] [28].
 
void StudyWAIC ()
 KS: Get the Watanabe-Akaike information criterion (WAIC) [11] [14].
 

Private Attributes

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

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_Samplestotal number of samples
Filenamename of output file
samplepointer to sample PDF object
nChainStepsnumber of steps in a chain, 0 indicate prior predictive was used

Definition at line 9 of file SampleSummary.cpp.

9 {
10// *******************
11 MACH3LOG_DEBUG("Making sample summary class...");
12 #ifdef MULTITHREAD
13 MACH3LOG_DEBUG("With OpenMP and {} threads", omp_get_max_threads());
14 #endif
15
17
18 if(StandardFluctuation) MACH3LOG_INFO("Using standard method of statistical fluctuation");
19 else MACH3LOG_INFO("Using alternative method of statistical fluctuation, which is much slower");
20
21 //KS: If true it will print posterior predictive for every beta parameter it is quite useful but make root big number of plots
22 DoBetaParam = true;
23 if(DoBetaParam) MACH3LOG_INFO("I will calculate #beta parameters from Barlow-Beeston formalism");
24
25 //If true code will normalise each histogram, this way you can calculate shape only error. etc. pvalue will be completely wrong unfortunately
26 doShapeOnly = false;
27
28 nChainSteps = nSteps;
29 //KS: nChainSteps == 0 means we run PriorPredcitive
30 if(nChainSteps == 0) isPriorPredictive = true;
31 else isPriorPredictive = false;
32
33 OutputName = Filename;
34 nSamples = n_Samples;
35 SampleHandler = sample;
36
37 //Get mach3 modes from manager
39
40 nThrows = 0;
41 first_pass = true;
42 Outputfile = nullptr;
43 OutputTree = nullptr;
44 rnd = std::make_unique<TRandom3>();
45
46 DataHist.resize(nSamples);
49 NominalHist.resize(nSamples);
50 PosteriorHist.resize(nSamples);
51 W2NomHist.resize(nSamples);
52 w2Hist.resize(nSamples);
53
56
57 if(DoBetaParam) BetaHist.resize(nSamples);
58
59 maxBins.resize(nSamples);
60
61 lnLHist_Mean.resize(nSamples);
62 lnLHist_Mode.resize(nSamples);
64 MeanHist.resize(nSamples);;
66 ModeHist.resize(nSamples);
67 W2MeanHist.resize(nSamples);
68 W2ModeHist.resize(nSamples);
74
75 //KS: When a histogram is created with an axis lower limit greater or equal to its upper limit ROOT will automatically adjust histogram range
76 // https://root.cern.ch/doc/master/classTH1.html#auto-bin
77 lnLHist = std::make_unique<TH1D>("lnLHist_predpredfluc", "lnLHist_predpredfluc", 100, 1, -1);
78 lnLHist->SetDirectory(nullptr);
79 lnLHist->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Pred)");
80 lnLHist->GetYaxis()->SetTitle("Counts");
81
82 lnLHist_drawdata = std::make_unique<TH1D>("lnLHist_drawdata", "lnLHist_drawdata", 100, 1, -1);
83 lnLHist_drawdata->SetDirectory(nullptr);
84 lnLHist_drawdata->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
85 lnLHist_drawdata->GetYaxis()->SetTitle("Counts");
86
87 lnLHist_drawfluc = std::make_unique<TH1D>("lnLHist_drawpredfluc", "lnLHist_drawpredfluc", 100, 1, -1);
88 lnLHist_drawfluc->SetDirectory(nullptr);
89 lnLHist_drawfluc->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
90 lnLHist_drawfluc->GetYaxis()->SetTitle("Counts");
91
92 lnLHist_drawflucdraw = std::make_unique<TH1D>("lnLHist_drawflucdraw", "lnLHist_drawflucdraw", 100, 1, -1);
93 lnLHist_drawflucdraw->SetDirectory(nullptr);
94 lnLHist_drawflucdraw->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
95 lnLHist_drawflucdraw->GetYaxis()->SetTitle("Counts");
96
97 lnLDrawHist = std::make_unique<TH2D>("lnLDrawHist", "lnLDrawHist", 50, 1, -1, 50, 1, -1);
98 lnLDrawHist->SetDirectory(nullptr);
99 lnLDrawHist->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
100 lnLDrawHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
101
102 lnLFlucHist = std::make_unique<TH2D>("lnLFlucHist", "lnLFlucHist", 50, 1, -1, 50, 1, -1);
103 lnLFlucHist->SetDirectory(nullptr);
104 lnLFlucHist->GetXaxis()->SetTitle("-2LLH_{Draw Fluc, Draw}");
105 lnLFlucHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
106
107 lnLDrawHistRate = std::make_unique<TH2D>("lnLDrawHistRate", "lnLDrawHistRate", 50, 1, -1, 50, 1, -1);
108 lnLDrawHistRate->SetDirectory(nullptr);
109 lnLDrawHistRate->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
110 lnLDrawHistRate->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
111
112 //KS: This is silly as it assumes all samples uses same kinematics
113 lnLFlucHist_ProjectX = std::make_unique<TH2D>("lnLFlucHist_ProjectX", "lnLFlucHist_ProjectX", 50, 1, -1, 50, 1, -1);
114 lnLFlucHist_ProjectX->SetDirectory(nullptr);
115 lnLFlucHist_ProjectX->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
116 lnLFlucHist_ProjectX->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
117
118 // Holds the hist of random number draws, only works for posterior predictive
120 {
121 RandomHist = std::make_unique<TH1D>("RandomHist", "RandomHist", 100, 0, nChainSteps);
122 RandomHist->SetDirectory(nullptr);
123 RandomHist->GetXaxis()->SetTitle("Step");
124 const double binwidth = nChainSteps/RandomHist->GetNbinsX();
125 std::stringstream ss;
126 ss << "Draws/" << binwidth;
127 RandomHist->GetYaxis()->SetTitle(ss.str().c_str());
128 RandomHist->SetLineWidth(2);
129 }
130 else RandomHist = nullptr;
131
132 for (int i = 0; i < nSamples; ++i)
133 {
134 DataHist[i] = nullptr;
135 DataHist_ProjectX[i] = nullptr;
136 DataHist_ProjectY[i] = nullptr;
137 NominalHist[i] = nullptr;
138
139 MeanHist[i] = nullptr;
140 if(DoBetaParam) MeanHistCorrected[i] = nullptr;
141 W2MeanHist[i] = nullptr;
142 W2ModeHist[i] = nullptr;
143 lnLHist_Mean[i] = nullptr;
144 lnLHist_Mode[i] = nullptr;
145 lnLHist_Mean_ProjectX[i] = nullptr;
146 lnLHist_Mean1D[i] = nullptr;
147 lnLHist_Mode1D[i] = nullptr;
148 lnLHist_Sample_DrawData[i] = nullptr;
149 lnLHist_Sample_DrawflucDraw[i] = nullptr;
150 lnLHist_Sample_PredflucDraw[i] = nullptr;
151 }//end loop over samples
152
153 DoByModePlots = false;
154 PosteriorHist_ByMode = nullptr;
155
156 nModelParams = 0;
157
158 Debug = 0;
159}
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
virtual std::string GetKinVarLabel(const int sample, const int Dimension)
std::vector< TH2Poly * > lnLHist_Mode
The LLH distribution in pmu cosmu for using the mode in each bin.
SampleHandlerBase * SampleHandler
Pointer to SampleHandler object, mostly used to get sample names, binning etc.
std::vector< TH2Poly * > lnLHist_Mean
The LLH distribution in pmu cosmu for using the mean in each bin.
std::vector< TH2Poly * > W2MeanHist
Pointer to the w2 histograms (for mean values).
std::unique_ptr< TH1D > RandomHist
Holds the history of which entries have been drawn in the MCMC file.
std::unique_ptr< TRandom3 > rnd
Random number generator.
std::vector< TH1D * > lnLHist_Mode1D
Holds the bin-by-bin LLH for the mode posterior predictive vs the data.
std::string OutputName
Output filename.
std::vector< std::vector< std::unique_ptr< TH1D > > > BetaHist
Distribution of beta parameters in Barlow Beeston formalisms.
int nModelParams
Number of parameters.
std::vector< TH2D * > ViolinHists_ProjectX
Posterior predictive but for X projection but as a violin plot.
bool first_pass
KS: Hacky flag to let us know if this is first toy.
std::vector< TH2Poly * > NominalHist
The nominal histogram for the selection.
std::vector< std::vector< std::unique_ptr< TH1D > > > PosteriorHist
The posterior predictive for the whole selection: this gets built after adding in the toys....
std::vector< TH2Poly * > MeanHistCorrected
The posterior predictive distribution in pmu cosmu using the mean after applying Barlow-Beeston Corre...
MaCh3Modes * Modes
MaCh3 Modes.
std::unique_ptr< TH1D > lnLHist
The histogram containing the lnL for each throw.
std::vector< TH1D * > lnLHist_Mean1D
Holds the bin-by-bin LLH for the mean posterior predictive vs the data.
TTree * OutputTree
TTree which we save useful information to.
bool DoByModePlots
By mode variables.
std::vector< TH1D * > lnLHist_Mean_ProjectX
The LLH distribution in pmu using the mean in each bin.
std::vector< TH2Poly * > DataHist
The data histogram for the selection.
int Debug
Tells Debug level to save additional histograms.
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.
std::vector< TH2D * > ViolinHists_ProjectY
Posterior predictive but for Y projection but as a violin plot.
std::vector< TH2Poly * > ModeHist
The posterior predictive distribution in pmu cosmu using the mode.
std::vector< TH2Poly * > W2NomHist
Pointer to the w2 histograms (for nominal values).
unsigned int nThrows
Number of throws.
std::vector< TH1D * > lnLHist_Sample_DrawData
The histogram containing the lnL (draw vs data) for each throw for each sample.
std::vector< std::vector< std::unique_ptr< TH1D > > > w2Hist
The posterior predictive for the whole selection: this gets built after adding in the toys....
TH1D **** PosteriorHist_ByMode
Histogram which corresponds to each bin in the sample's th2poly.
std::unique_ptr< TH1D > lnLHist_drawflucdraw
The lnLhist for the draw vs draw fluctuated.
TFile * Outputfile
Output filename.
std::unique_ptr< TH2D > lnLDrawHistRate
The 2D lnLhist, showing (draw vs data) and (draw vs fluct), using rate, anything above y=x axis is th...
std::unique_ptr< TH1D > lnLHist_drawdata
The lnLhist for the draw vs data.
int nSamples
Number of samples.
std::vector< TH1D * > DataHist_ProjectX
The data histogram for the selection X projection.
std::vector< TH1D * > DataHist_ProjectY
The data histogram for the selection Y projection.
bool isPriorPredictive
bool whether we have Prior or Posterior Predictive
bool DoBetaParam
Are we making Beta Histograms.
std::vector< TH2Poly * > W2ModeHist
Pointer to the w2 histograms (for mode values).
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-val...
bool StandardFluctuation
KS: We have two methods for Poissonian fluctuation.
std::unique_ptr< TH2D > lnLFlucHist_ProjectX
The 2D lnLHist but for ProjectionX histogram (pmu), showing (draw vs data) and (draw vs draw fluct),...
std::unique_ptr< TH1D > lnLHist_drawfluc
The lnLhist for the draw vs MC fluctuated.
bool doShapeOnly
bool whether to normalise each toy to have shape based p-value and pos pred distribution
std::vector< TH1D * > lnLHist_Sample_PredflucDraw
The histogram containing the lnL (draw vs pred fluct) for each throw for each sample.
std::vector< TH2Poly * > MeanHist
The posterior predictive distribution in pmu cosmu using the mean.
unsigned int nChainSteps
Number of throws by user.
std::vector< int > maxBins
Max Number of Bins per each sample.
std::vector< TH1D * > lnLHist_Sample_DrawflucDraw
The histogram containing the lnL (draw vs draw fluct) for each throw for each sample.
MaCh3Modes * GetMaCh3Modes() const
Return pointer to MaCh3 modes.

◆ ~SampleSummary()

SampleSummary::~SampleSummary ( )

Destructor.

Definition at line 163 of file SampleSummary.cpp.

163 {
164// *******************
165 Outputfile->cd();
166
167 //ROOT is weird and once you write TFile claim ownership of histograms. Best is to first delete histograms and then close file
168 Outputfile->Close();
169 delete Outputfile;
170
171 if(DoByModePlots)
172 {
173 for (int i = 0; i < nSamples; ++i)
174 {
175 if(DataHist[i] == nullptr) continue;
176 for (int j = 0; j < Modes->GetNModes()+1; j++)
177 {
178 for (int k = 1; k <= maxBins[i]; ++k)
179 {
180 if(PosteriorHist_ByMode[i][j][k] != nullptr) delete PosteriorHist_ByMode[i][j][k];
181 }
182 delete[] PosteriorHist_ByMode[i][j];
183 if(MeanHist_ByMode[i][j] != nullptr) delete MeanHist_ByMode[i][j];
184 }
185 delete[] PosteriorHist_ByMode[i];
186 }
187 delete[] PosteriorHist_ByMode;
188 }
189
190 for (int i = 0; i < nSamples; ++i)
191 {
192 if(DataHist[i] == nullptr) continue;
193 if(DataHist[i] != nullptr) delete DataHist[i];
194 if(NominalHist[i] != nullptr) delete NominalHist[i];
195 if(MeanHist[i] != nullptr) delete MeanHist[i];
196 if(ModeHist[i] != nullptr) delete ModeHist[i];
197 if(DoBetaParam && MeanHistCorrected[i] != nullptr) delete MeanHistCorrected[i];
198 if(W2MeanHist[i] != nullptr) delete W2MeanHist[i];
199 if(W2ModeHist[i] != nullptr) delete W2ModeHist[i];
200
201 if(ViolinHists_ProjectX[i] != nullptr) delete ViolinHists_ProjectX[i];
202 if(ViolinHists_ProjectY[i] != nullptr) delete ViolinHists_ProjectY[i];
203
204 if(lnLHist_Mean[i] != nullptr) delete lnLHist_Mean[i];
205 if(lnLHist_Mode[i] != nullptr) delete lnLHist_Mode[i];
206 if(lnLHist_Mean_ProjectX[i] != nullptr) delete lnLHist_Mean_ProjectX[i];
207 if(lnLHist_Mean1D[i] != nullptr) delete lnLHist_Mean1D[i];
208 if(lnLHist_Mode1D[i] != nullptr) delete lnLHist_Mode1D[i];
209 if(lnLHist_Sample_DrawData[i] != nullptr) delete lnLHist_Sample_DrawData[i];
212 }
213}
int GetNModes() const
KS: Get number of modes, keep in mind actual number is +1 greater due to unknown category.
Definition: MaCh3Modes.h:52
std::vector< std::vector< TH2Poly * > > MeanHist_ByMode
The posterior predictive distribution in pmu cosmu using the mean.

Member Function Documentation

◆ AddData()

void SampleSummary::AddData ( std::vector< TH2Poly * > &  DataHist)

KS: Add data histograms.

Parameters
DataHistHistogram with data even rates for each sample

Definition at line 233 of file SampleSummary.cpp.

233 {
234// *******************
235 const int Length = int(Data.size());
236 // Check length of samples are OK
237 if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
238 for (int i = 0; i < Length; ++i) {
239 if (Data[i] == nullptr) {
240 DataHist[i] = nullptr;
241 DataHist_ProjectX[i] = nullptr;
242 DataHist_ProjectY[i] = nullptr;
243 maxBins[i] = 0;
244 } else {
245 std::string classname = std::string(DataHist[i]->Class_Name());
246 if(classname == "TH2Poly")
247 {
248 DataHist[i] = static_cast<TH2Poly*>(Data[i]->Clone());
250 DataHist_ProjectX[i] = ProjectPoly(DataHist[i], true, i);
251 DataHist_ProjectY[i] = ProjectPoly(DataHist[i], false, i);
252 maxBins[i] = DataHist[i]->GetNumberOfBins();
253 } else {
254 MACH3LOG_ERROR("Somehow sample {} doesn't use TH2Poly", SampleHandler->GetSampleName(i));
255 MACH3LOG_ERROR("Right now I only support TH2Poly but I am ambitious piece of code and surely will have more support in the future");
256 throw MaCh3Exception(__FILE__ , __LINE__ );
257 }
258 }
259 }
260}
void NormaliseTH2Poly(TH2Poly *Histogram)
Helper to Normalise histograms.
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
Custom exception class for MaCh3 errors.
TH1D * ProjectPoly(TH2Poly *Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist=false)
Helper to project TH2Poly onto axis.
bool CheckSamples(const int Length)
Check the length of samples agrees.
virtual std::string GetSampleName(int Sample) const =0

◆ AddNominal()

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

KS: Add prior histograms.

Definition at line 264 of file SampleSummary.cpp.

264 {
265// *******************
266 const int Length = int(Nominal.size());
267 if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
268
269 //KS: ROOT is super annoying and you cannot use clone with openMP, hence we have another loop below
270 for (int i = 0; i < Length; ++i)
271 {
272 if (Nominal[i] == nullptr) {
273 NominalHist[i] = nullptr;
274 W2NomHist[i] = nullptr;
275 lnLHist_Mean[i] = nullptr;
276 lnLHist_Mode[i] = nullptr;
277 lnLHist_Mean_ProjectX[i] = nullptr;
278 MeanHist[i] = nullptr;
279 if(DoBetaParam) MeanHistCorrected[i] = nullptr;
280 ModeHist[i] = nullptr;
281 W2MeanHist[i] = nullptr;
282 W2ModeHist[i] = nullptr;
283 lnLHist_Sample_DrawData[i] = nullptr;
284 lnLHist_Sample_DrawflucDraw[i] = nullptr;
285 lnLHist_Sample_PredflucDraw[i] = nullptr;
286 // If not nullptr it indicates the selection was turned on, so initialise the privates
287 } else {
288 NominalHist[i] = static_cast<TH2Poly*>(Nominal[i]->Clone());
290 W2NomHist[i] = static_cast<TH2Poly*>(NomW2[i]->Clone());
291
292 lnLHist_Mean[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
293 lnLHist_Mean[i]->SetDirectory(nullptr);
294 lnLHist_Mode[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
295 lnLHist_Mode[i]->SetDirectory(nullptr);
296 lnLHist_Mean_ProjectX[i] = static_cast<TH1D*>(DataHist_ProjectX[i]->Clone());
297 MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
298 if(DoBetaParam) MeanHistCorrected[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
299 ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
300 W2MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
301 W2ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
302 }
303 }
304
305 // Loop over the length of nominal and set the initial distributions up
306 //KS: Don't multithread, mostly due to fact that we initialise histograms
307 for (int i = 0; i < Length; ++i) {
308 // If NULL it indicates the selection was turned off, so initialise all the hists to NULL
309 if (Nominal[i] != nullptr)
310 {
311 std::string name = std::string(NominalHist[i]->GetName());
312 name = name.substr(0, name.find("_nom"));
313
314 PosteriorHist[i].resize(maxBins[i]+1);
315 w2Hist[i].resize(maxBins[i]+1);
316
317 if(DoBetaParam) BetaHist[i].resize(maxBins[i]+1);
318
319 for (int j = 0; j <= maxBins[i]; ++j)
320 {
321 PosteriorHist[i][j] = nullptr;
322 }
323 lnLHist_Mean[i]->SetNameTitle((name+"_MeanlnL").c_str(), (name+"_MeanlnL").c_str());
324 lnLHist_Mean[i]->Reset("");
325 lnLHist_Mean[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
326
327 lnLHist_Mode[i]->SetNameTitle((name+"_ModelnL").c_str(), (name+"_ModelnL").c_str());
328 lnLHist_Mode[i]->Reset("");
329 lnLHist_Mode[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
330
331 lnLHist_Mean_ProjectX[i]->SetNameTitle((name+"_MeanlnL_ProjectX").c_str(), (name+"_MeanlnL_ProjectX").c_str());
332 lnLHist_Mean_ProjectX[i]->Reset("");
333 lnLHist_Mean_ProjectX[i]->GetYaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
334
335 MeanHist[i]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
336 MeanHist[i]->Reset("");
337 MeanHist[i]->GetZaxis()->SetTitle("Mean");
338
339 if(DoBetaParam)
340 {
341 MeanHistCorrected[i]->SetNameTitle((name+"_mean_corrected").c_str(), (name+"_mean_corrected").c_str());
342 MeanHistCorrected[i]->Reset("");
343 MeanHistCorrected[i]->GetZaxis()->SetTitle("Mean");
344 }
345 std::vector<double> xbins;
346 std::vector<double> ybins;
347
348 SampleHandler->SetupBinning(M3::int_t(i), xbins, ybins);
349
350 //KS: Y axis is number of events to get estimate of maximal number we use integral
351 const int MaxBinning = doShapeOnly ? 1 : int(NoOverflowIntegral(NominalHist[i])/4);
352 ViolinHists_ProjectX[i] = new TH2D((name+"_Violin_ProjectX").c_str(), (name+"_Violin_ProjectX").c_str(), int(xbins.size()-1), &xbins[0] , 400, 0, MaxBinning);
353 ViolinHists_ProjectX[i]->GetYaxis()->SetTitle("Events");
354 ViolinHists_ProjectX[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetXaxis()->GetTitle()).c_str() );
355 ViolinHists_ProjectX[i]->SetDirectory(nullptr);
356
357 ViolinHists_ProjectY[i] = new TH2D((name+"_Violin_ProjectY").c_str(), (name+"_Violin_ProjectY").c_str(), int(ybins.size()-1), &ybins[0] , 400, 0, MaxBinning);
358 ViolinHists_ProjectY[i]->GetYaxis()->SetTitle("Events");
359 ViolinHists_ProjectY[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetYaxis()->GetTitle()).c_str());
360 ViolinHists_ProjectY[i]->SetDirectory(nullptr);
361
362 ModeHist[i]->SetNameTitle((name+"_mode").c_str(), (name+"_mode").c_str());
363 ModeHist[i]->Reset("");
364 ModeHist[i]->GetZaxis()->SetTitle("Mode");
365
366 W2MeanHist[i]->SetNameTitle((name+"_w2mean").c_str(), (name+"_w2mean").c_str());
367 W2MeanHist[i]->Reset("");
368 W2MeanHist[i]->GetZaxis()->SetTitle("W2Mean");
369
370 W2ModeHist[i]->SetNameTitle((name+"_w2mode").c_str(), (name+"_w2mode").c_str());
371 W2ModeHist[i]->Reset("");
372 W2ModeHist[i]->GetZaxis()->SetTitle("W2Mode");
373
374 // Declare the lnL histograms
375 lnLHist_Mean1D[i] = new TH1D((name+"_MeanlnL1D").c_str(),(name+"_MeanlnL1D").c_str(), 50, 1, -1);
376 lnLHist_Mean1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
377 lnLHist_Mean1D[i]->GetYaxis()->SetTitle("Counts");
378
379 lnLHist_Mode1D[i] = new TH1D((name+"_ModelnL1D").c_str(),(name+"_ModelnL1D").c_str(), 50, 1, -1);
380 lnLHist_Mode1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
381 lnLHist_Mode1D[i]->GetYaxis()->SetTitle("Counts");
382
383 lnLHist_Sample_DrawData[i] = new TH1D((name+"_lnLdrawdata").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
384 lnLHist_Sample_DrawData[i]->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
385 lnLHist_Sample_DrawData[i]->GetYaxis()->SetTitle("Counts");
386
387 lnLHist_Sample_DrawflucDraw[i] = new TH1D((name+"_lnLdrawfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
388 lnLHist_Sample_DrawflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
389 lnLHist_Sample_DrawflucDraw[i]->GetYaxis()->SetTitle("Counts");
390
391 lnLHist_Sample_PredflucDraw[i] = new TH1D((name+"_lnLpredfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
392 lnLHist_Sample_PredflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
393 lnLHist_Sample_PredflucDraw[i]->GetYaxis()->SetTitle("Counts");
394 }
395 }
396 //KS: Separate loop for thread safe reasons
397 for (int i = 0; i < Length; ++i)
398 {
399 //KS: We copy histograms so delete original
400 delete Nominal[i];
401 delete NomW2[i];
402 }
403}
double NoOverflowIntegral(TH2Poly *poly)
WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow.
virtual void SetupBinning(const M3::int_t Selection, std::vector< double > &BinningX, std::vector< double > &BinningY)
int int_t
Definition: Core.h:29

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

408 {
409// *******************
410 nThrows++;
411 //KS: Only make sense for PosteriorPredictive
412 if( !isPriorPredictive )RandomHist->Fill(DrawNumber);
413
414 const int size = int(SampleVector.size());
415 if (!CheckSamples(size)) throw MaCh3Exception(__FILE__ , __LINE__ );
416
417 // Push back the throw
418 MCVector.push_back(SampleVector);
419 LLHPenaltyVector.push_back(LLHPenalty);
420 WeightVector.push_back(Weight);
421 W2MCVector.push_back(W2Vec);
422
423 // Initialise the posterior hist
424 if (first_pass)
425 {
426 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
427 {
428 const int nXBins = 500;
429 //Initialise TH1D which corresponds to each bin in the sample's th2poly
430 std::string name = std::string(SampleVector[SampleNum]->GetName());
431 for (int i = 1; i <= maxBins[SampleNum]; ++i)
432 {
433 //Get PolyBin
434 TH2PolyBin* bin = static_cast<TH2PolyBin*>(SampleVector[SampleNum]->GetBins()->At(i-1));
435
436 // Just make a little fancy name
437 std::stringstream ss2;
438 ss2 << name << "_";
439 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
440 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
441
442 PosteriorHist[SampleNum][i] = std::make_unique<TH1D>(ss2.str().c_str(), ss2.str().c_str(),nXBins, 1, -1);
443 PosteriorHist[SampleNum][i]->SetDirectory(nullptr);
444 w2Hist[SampleNum][i] = std::make_unique<TH1D>(("w2_"+ss2.str()).c_str(), ("w2_"+ss2.str()).c_str(),nXBins, 1, -1);
445 w2Hist[SampleNum][i]->SetDirectory(nullptr);
446 if(DoBetaParam)
447 {
448 std::string betaName = "#beta_param_";
449 BetaHist[SampleNum][i] = std::make_unique<TH1D>((betaName + ss2.str()).c_str(), (betaName + ss2.str()).c_str(), 70, 1, -1);
450 BetaHist[SampleNum][i]->SetDirectory(nullptr);
451 BetaHist[SampleNum][i]->GetXaxis()->SetTitle("#beta parameter value");
452 BetaHist[SampleNum][i]->GetYaxis()->SetTitle("Counts");
453 }
454 } //end loop over bins
455 }//end loop over samples
456 }
457 first_pass = false;
458
459 // Loop over the samples
460 #ifdef MULTITHREAD
461 #pragma omp parallel for
462 #endif
463 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
464 {
465 if (SampleVector[SampleNum] == nullptr) continue;
466 if(doShapeOnly) NormaliseTH2Poly(SampleVector[SampleNum]);
467 // Loop over the distribution and fill the prior/posterior predictive
468 for (int i = 1; i <= maxBins[SampleNum]; ++i) {
469 const double Content = SampleVector[SampleNum]->GetBinContent(i);
470 PosteriorHist[SampleNum][i]->Fill(Content, Weight);
471 const double w2 = W2Vec[SampleNum]->GetBinContent(i);
472 w2Hist[SampleNum][i]->Fill(w2, Weight);
473 if(DoBetaParam)
474 {
475 const double data = DataHist[SampleNum]->GetBinContent(i);
476 const double BetaParam = GetBetaParameter(data, Content, w2, likelihood);
477 BetaHist[SampleNum][i]->Fill(BetaParam, Weight);
478 }
479 } // end bin loop
480 } // end samples loop
481} // end AddThrow
int size
double GetBetaParameter(const double data, const double mc, const double w2, TestStatistic TestStat)
KS: Calculate Beta parameter which will be different based on specified test statistic.
std::vector< double > LLHPenaltyVector
Vector to hold the penalty term.
std::vector< std::vector< TH2Poly * > > MCVector
Vector of vectors which holds the loaded MC histograms.
std::vector< std::vector< TH2Poly * > > W2MCVector
Vector of vectors which holds the loaded W2 histograms.
TestStatistic likelihood
Type of likelihood for example Poisson, Barlow-Beeston or Ice Cube.
std::vector< double > WeightVector
Vector holding weight.

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

486 {
487// *******************
488 MCVectorByMode.push_back(SampleVector_ByMode);
489
490 //KS: This means this is first time
491 if(!DoByModePlots)
492 {
493 MACH3LOG_INFO("Turning reaction breadkwon mode, brum brum");
494 PosteriorHist_ByMode = new TH1D***[nSamples];
496 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
497 {
498 if (DataHist[SampleNum] == nullptr) continue;
499
500 PosteriorHist_ByMode[SampleNum] = new TH1D**[Modes->GetNModes()+1];
501 MeanHist_ByMode[SampleNum].resize(Modes->GetNModes()+1);
502 for (int j = 0; j < Modes->GetNModes()+1; j++)
503 {
504 PosteriorHist_ByMode[SampleNum][j] = new TH1D*[maxBins[SampleNum]+1];
505 constexpr int nXBins = 500;
506
507 std::string name = std::string(NominalHist[SampleNum]->GetName());
508 name = name.substr(0, name.find("_nom"));
509 name = name + "_"+Modes->GetMaCh3ModeName(j);
510
511 for (int i = 1; i <= maxBins[SampleNum]; i++)
512 {
513 //Get PolyBin
514 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i-1));
515
516 // Just make a little fancy name
517 std::stringstream ss2;
518 ss2 << name << "_";
519 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
520 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
521
522 //Initialise TH1D which corresponds to each bin in the sample's th2poly
523 PosteriorHist_ByMode[SampleNum][j][i] = new TH1D((name+ss2.str()).c_str(),(name+ss2.str()).c_str(),nXBins, 1, -1);
524 }
525 MeanHist_ByMode[SampleNum][j] = static_cast<TH2Poly*>(NominalHist[SampleNum]->Clone());
526 MeanHist_ByMode[SampleNum][j]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
527 MeanHist_ByMode[SampleNum][j]->Reset("");
528 MeanHist_ByMode[SampleNum][j]->GetZaxis()->SetTitle("Mean");
529 }
530 }
531 }
532 DoByModePlots = true;
533 // Loop over the sameples
534 #ifdef MULTITHREAD
535 #pragma omp parallel for
536 #endif
537 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
538 {
539 if (DataHist[SampleNum] == nullptr) continue;
540
541 for (int j = 0; j < Modes->GetNModes()+1; j++)
542 {
543 if(doShapeOnly) NormaliseTH2Poly(SampleVector_ByMode[SampleNum][j]);
544 // Loop over the distribution and fill the prior/posterior predictive
545 for (int i = 1; i <= maxBins[SampleNum]; ++i)
546 {
547 const double Content = SampleVector_ByMode[SampleNum][j]->GetBinContent(i);
548 const int Entries = int(PosteriorHist_ByMode[SampleNum][j][i]->GetEntries());
549 PosteriorHist_ByMode[SampleNum][j][i]->Fill(Content, WeightVector[Entries]);
550 }
551 }
552 }
553} // end AddThrowByMode
std::string GetMaCh3ModeName(const int Index)
KS: Get normal name of mode, if mode not known you will get UNKNOWN_BAD.
Definition: MaCh3Modes.cpp:138
std::vector< std::vector< std::vector< TH2Poly * > > > MCVectorByMode
Vector of vectors which holds the loaded MC histograms for each mode.

◆ 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 tittle to include LLH.

Parameters
Datahistogram with data distribution for a single sample
MChistogram with MC distribution for a single sample
W2histogram with W2 distribution for a single sample

Definition at line 1561 of file SampleSummary.cpp.

1561 {
1562// ****************
1563 const double llh = GetLLH(DatHist, MCHist, W2Hist);
1564 std::stringstream ss;
1565 ss << "_2LLH=" << llh;
1566 MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
1567 MACH3LOG_INFO("{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), DatHist->Integral(), MCHist->Integral(), llh);
1568}
double GetLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly.

◆ 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 tittle to include LLH.

Parameters
Datahistogram with data distribution for a single sample
MChistogram with MC distribution for a single sample
W2histogram with W2 distribution for a single sample

Definition at line 1572 of file SampleSummary.cpp.

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

◆ CheckSamples()

bool SampleSummary::CheckSamples ( const int  Length)
inlineprivate

Check the length of samples agrees.

Definition at line 217 of file SampleSummary.cpp.

217 {
218// *******************
219 bool ok = (nSamples == Length);
220 if (!ok) {
221 MACH3LOG_ERROR("Size of SampleVector input != number of defined samples");
222 MACH3LOG_ERROR("Size of SampleVector: {}", Length);
223 MACH3LOG_ERROR("Size of defined samples: {}", nSamples);
224 MACH3LOG_ERROR("Something has gone wrong with making the Samples");
225 throw MaCh3Exception(__FILE__ , __LINE__ );
226 }
227 return ok;
228}

◆ GetLLH() [1/2]

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

Helper functions to calculate likelihoods using TH1D.

Parameters
Datahistogram with data distribution for a single sample
MChistogram with MC distribution for a single sample
W2histogram with W2 distribution for a single sample

Definition at line 1597 of file SampleSummary.cpp.

1597 {
1598// ****************
1599 double llh = 0.0;
1600 for (int i = 1; i <= DatHist->GetXaxis()->GetNbins(); ++i)
1601 {
1602 const double data = DatHist->GetBinContent(i);
1603 const double mc = MCHist->GetBinContent(i);
1604 const double w2 = W2Hist->GetBinContent(i);
1605 llh += SampleHandler->GetTestStatLLH(data, mc, w2);
1606 }
1607 //KS: do times 2 because banff reports chi2
1608 return 2*llh;
1609}
double GetTestStatLLH(double data, double mc) const
Calculate test statistic for a single bin using Poisson.

◆ GetLLH() [2/2]

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

Helper functions to calculate likelihoods using TH2Poly.

Parameters
Datahistogram with data distribution for a single sample
MChistogram with MC distribution for a single sample
W2histogram with W2 distribution for a single sample

Definition at line 1582 of file SampleSummary.cpp.

1582 {
1583// ****************
1584 double llh = 0.0;
1585 for (int i = 1; i < DatHist->GetNumberOfBins()+1; ++i)
1586 {
1587 const double data = DatHist->GetBinContent(i);
1588 const double mc = MCHist->GetBinContent(i);
1589 const double w2 = W2Hist->GetBinContent(i);
1590 llh += SampleHandler->GetTestStatLLH(data, mc, w2);
1591 }
1592 //KS: do times 2 because banff reports chi2
1593 return 2*llh;
1594}

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

1134 {
1135// *******************
1136 MACH3LOG_INFO("Making the chi2 histograms...");
1137 // Have this to signify if we're doing the first pass
1138 first_pass = true;
1139
1140 double AveragePenalty = 0;
1141
1142 // Vectors to hold exact LLH
1143 std::vector<double> LLH_PredFluc_V(nThrows);
1144 std::vector<double> LLH_DataDraw_V(nThrows);
1145 std::vector<double> LLH_DrawFlucDraw_V(nThrows);
1146
1147 // Loop over the draws
1148 // Should look into multi-threading this. Would require temporary THxx structures from arrays
1149 //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
1150 for (unsigned int i = 0; i < nThrows; ++i)
1151 {
1152 if (i % (nThrows/10) == 0) {
1154 }
1155
1156 // Set the total LLH to zero to initialise
1157 double total_llh_data_draw_temp = 0.0;
1158 double total_llh_drawfluc_draw_temp = 0.0;
1159 double total_llh_predfluc_draw_temp = 0.0;
1160
1161 double total_llh_rate_data_draw_temp = 0.0;
1162 double total_llh_rate_predfluc_draw_temp = 0.0;
1163
1164 double total_llh_data_drawfluc_temp = 0.0;
1165 double total_llh_data_predfluc_temp = 0.0;
1166 double total_llh_draw_pred_temp = 0.0;
1167 double total_llh_drawfluc_pred_temp = 0.0;
1168 double total_llh_drawfluc_predfluc_temp = 0.0;
1169 double total_llh_predfluc_pred_temp = 0.0;
1170 double total_llh_datafluc_draw_temp = 0.0;
1171
1172 double total_llh_data_draw_ProjectX_temp = 0.0;
1173 double total_llh_drawfluc_draw_ProjectX_temp = 0.0;
1174
1175 // Save the double that gets written to file
1177 AveragePenalty += llh_penalty;
1178
1179 // Make the Poisson fluctuated hist
1180 std::vector<TH2Poly*> FluctHist(nSamples);
1181 // Also Poisson fluctuate the drawn MCMC hist
1182 std::vector<TH2Poly*> FluctDrawHist(nSamples);
1183 // Finally Poisson fluctuate the data histogram
1184 std::vector<TH2Poly*> DataFlucHist(nSamples);
1185
1186 // Finally Poisson fluctuate the data histogram
1187 std::vector<TH1D*> FluctDrawHistProjectX(nSamples);
1188 std::vector<TH1D*> DrawHistProjectX(nSamples);
1189 std::vector<TH1D*> DrawHistProjectY(nSamples);
1190 std::vector<TH1D*> DrawW2HistProjectX(nSamples);
1191
1192 //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
1193 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1194 {
1195 FluctHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1196 FluctDrawHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1197 DataFlucHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1198
1199 FluctDrawHistProjectX[SampleNum] = static_cast<TH1D*>(DataHist_ProjectX[SampleNum]->Clone());
1200
1201 // Get the ith draw for the jth sample
1202 TH2Poly *DrawHist = MCVector[i][SampleNum];
1203 TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
1204
1205 //ProjectPoly calls new TH1D under the hood, never define new ROOT object under MP...
1206 DrawHistProjectX[SampleNum] = ProjectPoly(DrawHist, true, SampleNum);
1207 DrawW2HistProjectX[SampleNum] = ProjectPoly(DrawW2Hist, true, SampleNum);
1208 DrawHistProjectY[SampleNum] = ProjectPoly(DrawHist, false, SampleNum);
1209 }
1210 #ifdef MULTITHREAD
1211 //KS: might be most obscure OMP reduction I have ever made...
1212 #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)
1213 #endif
1214 // Loop over the samples
1215 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1216 {
1217 // Get the ith draw for the jth sample
1218 TH2Poly *DrawHist = MCVector[i][SampleNum];
1219 TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
1220 // Skip empty samples
1221 if (DrawHist == nullptr) continue;
1222
1223 // Add LLH penalties from the systematics to the LLH that use the drawn histogram
1224 // Data vs Draw
1225 llh_data_draw[SampleNum] = llh_penalty;
1226 // Fluctuated Draw vs Draw
1227 llh_drawfluc_draw[SampleNum] = llh_penalty;
1228 // Fluctuated Predicitve vs Draw
1229 llh_predfluc_draw[SampleNum] = llh_penalty;
1230
1231 // Data vs Draw using rate
1232 llh_rate_data_draw[SampleNum] = llh_penalty;
1233 // Fluctuated Predicitve vs Draw using rate
1235
1236 // Data vs Fluctuated Draw
1237 llh_data_drawfluc[SampleNum] = llh_penalty;
1238 // Draw vs Predictive
1239 llh_draw_pred[SampleNum] = llh_penalty;
1240 // Fluctuated Draw vs Predictive
1241 llh_drawfluc_pred[SampleNum] = llh_penalty;
1242 // Fluctuated Draw vs Fluctuated Predictive
1244 // Fluctuated Data vs Draw
1245 llh_datafluc_draw[SampleNum] = llh_penalty;
1246
1247 //Some LLH for 1D projections
1250
1251 //Other get 0 penalty term
1252 // Fluctuated Predictive vs Predictive
1253 llh_predfluc_pred[SampleNum] = 0.0;
1254 // Data vs Fluctuated Predictive
1255 llh_data_predfluc[SampleNum] = 0.0;
1256
1257 // Make the Poisson fluctuated hist
1258 MakeFluctuatedHistogram(FluctHist[SampleNum], MeanHist[SampleNum]);
1259 // Also Poisson fluctuate the drawn MCMC hist
1260 MakeFluctuatedHistogram(FluctDrawHist[SampleNum], DrawHist);
1261 // Finally Poisson fluctuate the data histogram
1262 MakeFluctuatedHistogram(DataFlucHist[SampleNum], DataHist[SampleNum]);
1263
1264 // Likelihood between the drawn histogram and the data
1265 const double DataDrawLLH = GetLLH(DataHist[SampleNum], DrawHist, DrawW2Hist);
1266 llh_data_draw[SampleNum] += DataDrawLLH;
1267 total_llh_data_draw_temp += DataDrawLLH;
1268
1269 // Likelihood between drawn histogram and fluctuated drawn histogram
1270 const double DrawFlucDrawLLH = GetLLH(FluctDrawHist[SampleNum], DrawHist, DrawW2Hist);
1271 llh_drawfluc_draw[SampleNum] += DrawFlucDrawLLH;
1272 total_llh_drawfluc_draw_temp += DrawFlucDrawLLH;
1273
1274 // Likelihood between drawn histogram and fluctuated posterior predictive distribution
1275 const double PredFlucDrawLLH = GetLLH(FluctHist[SampleNum], DrawHist, DrawW2Hist);
1276 llh_predfluc_draw[SampleNum] += PredFlucDrawLLH;
1277 total_llh_predfluc_draw_temp += PredFlucDrawLLH;
1278
1279//Rate Based p-value
1280 // Likelihood between the drawn histogram and the data
1281 const double RateDataDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(DataHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
1282 llh_rate_data_draw[SampleNum] += RateDataDrawLLH;
1283 total_llh_rate_data_draw_temp += RateDataDrawLLH;
1284
1285 // Likelihood between drawn histogram and fluctuated posterior predictive distribution using rate
1286 const double RatePredFlucDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(FluctHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
1287 llh_rate_predfluc_draw[SampleNum] += RatePredFlucDrawLLH;
1288 total_llh_rate_predfluc_draw_temp += RatePredFlucDrawLLH;
1289
1290// All LLH below are for validation reason but not used for final P-Value
1291 // Likelihood between the fluctuated drawn histogram and the data
1292 const double DataDrawFlucLLH = GetLLH(DataHist[SampleNum], FluctDrawHist[SampleNum], DrawW2Hist);
1293 llh_data_drawfluc[SampleNum] += DataDrawFlucLLH;
1294 total_llh_data_drawfluc_temp += DataDrawFlucLLH;
1295
1296 // Likelihood between the drawn histogram and the data
1297 const double DataPredFlucLLH = GetLLH(DataHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
1298 llh_data_predfluc[SampleNum] += DataPredFlucLLH;
1299 total_llh_data_predfluc_temp += DataPredFlucLLH;
1300
1301 // Likelihood between the drawn hist and the Posterior Predictive
1302 const double DrawPredLLH = GetLLH(DrawHist, MeanHist[SampleNum], W2MeanHist[SampleNum]);
1303 llh_draw_pred[SampleNum] += DrawPredLLH;
1304 total_llh_draw_pred_temp += DrawPredLLH;
1305
1306 // Likelihood between fluctuated drawn and predictive
1307 const double DrawFlucPredLLH = GetLLH(FluctDrawHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
1308 llh_drawfluc_pred[SampleNum] += DrawFlucPredLLH;
1309 total_llh_drawfluc_pred_temp += DrawFlucPredLLH;
1310
1311 // Likelihood between drawn histogram and fluctuated drawn histogram
1312 const double DrawFlucPredFlucLLH = GetLLH(FluctDrawHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
1313 llh_drawfluc_predfluc[SampleNum] += DrawFlucPredFlucLLH;
1314 total_llh_drawfluc_predfluc_temp += DrawFlucPredFlucLLH;
1315
1316 // Likelihood between the fluctuated drawn histogram and the posterior predictive
1317 const double PredFlucPredLLH = GetLLH(FluctHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
1318 llh_predfluc_pred[SampleNum] += PredFlucPredLLH;
1319 total_llh_predfluc_pred_temp += PredFlucPredLLH;
1320
1321 // Likelihood between fluctuated data histogram and drawn histogram
1322 const double DataFlucDrawLLH = GetLLH(DataFlucHist[SampleNum], DrawHist, DrawW2Hist);
1323 llh_datafluc_draw[SampleNum] += DataFlucDrawLLH;
1324 total_llh_datafluc_draw_temp += DataFlucDrawLLH;
1325
1326 lnLHist_Sample_DrawData[SampleNum]->Fill(DataDrawLLH);
1327 lnLHist_Sample_DrawflucDraw[SampleNum]->Fill(DrawFlucDrawLLH);
1328 lnLHist_Sample_PredflucDraw[SampleNum]->Fill(PredFlucDrawLLH);
1329
1330// At the end we leave LLH for 1D projections
1331 MakeFluctuatedHistogram(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum]);
1332
1333 // Likelihood between the drawn histogram and the data for muon momentum
1334 const double DataDrawLLH_ProjectX = GetLLH(DataHist_ProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1335 llh_data_draw_ProjectX[SampleNum] += DataDrawLLH_ProjectX;
1336 total_llh_data_draw_ProjectX_temp += DataDrawLLH_ProjectX;
1337
1338 const double DrawFlucDrawLLH_ProjectX = GetLLH(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1339 llh_drawfluc_draw_ProjectX[SampleNum] += DrawFlucDrawLLH_ProjectX;
1340 total_llh_drawfluc_draw_ProjectX_temp += DrawFlucDrawLLH_ProjectX;
1341
1342 //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
1343 FastViolinFill(ViolinHists_ProjectX[SampleNum], DrawHistProjectX[SampleNum]);
1344 FastViolinFill(ViolinHists_ProjectY[SampleNum], DrawHistProjectY[SampleNum]);
1345 } // End loop over samples (still looping throws)
1346
1347 // Delete the temporary histograms
1348 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1349 {
1350 delete FluctHist[SampleNum];
1351 delete FluctDrawHist[SampleNum];
1352 delete DataFlucHist[SampleNum];
1353 delete FluctDrawHistProjectX[SampleNum];
1354 delete DrawHistProjectX[SampleNum];
1355 delete DrawHistProjectY[SampleNum];
1356 delete DrawW2HistProjectX[SampleNum];
1357 }
1358
1359 total_llh_data_draw = total_llh_data_draw_temp;
1360 total_llh_drawfluc_draw = total_llh_drawfluc_draw_temp;
1361 total_llh_predfluc_draw = total_llh_predfluc_draw_temp;
1362
1363 total_llh_rate_data_draw = total_llh_rate_data_draw_temp;
1364 total_llh_rate_predfluc_draw = total_llh_rate_predfluc_draw_temp;
1365
1366 total_llh_data_drawfluc = total_llh_data_drawfluc_temp;
1367 total_llh_data_predfluc = total_llh_data_predfluc_temp;
1368 total_llh_draw_pred = total_llh_draw_pred_temp;
1369 total_llh_drawfluc_pred = total_llh_drawfluc_pred_temp;
1370 total_llh_drawfluc_predfluc = total_llh_drawfluc_predfluc_temp;
1371 total_llh_predfluc_pred = total_llh_predfluc_pred_temp;
1372 total_llh_datafluc_draw = total_llh_datafluc_draw_temp;
1373
1374 total_llh_data_draw_ProjectX = total_llh_data_draw_ProjectX_temp;
1375 total_llh_drawfluc_draw_ProjectX = total_llh_drawfluc_draw_ProjectX_temp;
1376
1377 // Add LLH penalties from the systematics to the LLH that use the drawn histogram
1381 //Rate based
1384
1389
1392
1397
1400
1402
1404
1405 // Also save to arrays to make sure we have the utmost super accuracy
1406 LLH_PredFluc_V[i] = total_llh_predfluc_draw;
1407 LLH_DataDraw_V[i] = total_llh_data_draw;
1408 LLH_DrawFlucDraw_V[i] = total_llh_drawfluc_draw;
1409
1410 // Write to the output tree
1411 OutputTree->Fill();
1412 } // End loop over throws
1413
1414 AveragePenalty = AveragePenalty/double(nThrows);
1415 MACH3LOG_INFO("Average LLH penalty over toys is {:.2f}", AveragePenalty);
1416 // Calculate exact p-value instead of binned
1417 unsigned int Accept_PredFluc = 0;
1418 unsigned int Accept_DrawFluc = 0;
1419 for (unsigned int i = 0; i < nThrows; ++i)
1420 {
1421 if (LLH_DataDraw_V[i] > LLH_DrawFlucDraw_V[i]) Accept_DrawFluc++;
1422 if (LLH_DataDraw_V[i] > LLH_PredFluc_V[i]) Accept_PredFluc++;
1423 }
1424 const double pvalue_DrawFluc = double(Accept_DrawFluc)/double(nThrows);
1425 const double pvalue_PredFluc = double(Accept_PredFluc)/double(nThrows);
1426
1427 MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Draw: {:.2f}", pvalue_DrawFluc);
1428 MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Prediction: {:.2f}", pvalue_PredFluc);
1429}
void FastViolinFill(TH2D *violin, TH1D *hist_1d)
KS: Fill Violin histogram with entry from a toy.
std::vector< double > llh_data_draw_ProjectX
Projection X (most likely muon momentum) of LLH.
std::vector< double > llh_predfluc_draw
Fluctuated Predictive vs Draw.
double total_llh_predfluc_pred
Fluctuated Predictive vs Predictive.
double total_llh_drawfluc_draw_ProjectX
Fluctuated Draw vs Draw for projection X (most likely muon momentum)
double total_llh_data_draw
Data vs Draw.
double total_llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using Rate.
std::vector< double > llh_data_predfluc
Data vs Fluctuated Predictive.
std::vector< double > llh_predfluc_pred
Fluctuated Predictive vs Predictive.
double total_llh_draw_pred
Draw vs Predictive.
std::vector< double > llh_drawfluc_draw
Fluctuated Draw vs Draw.
double total_llh_data_draw_ProjectX
Data vs Draw for projection X (most likely muon momentum)
std::vector< double > llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
std::vector< double > llh_data_draw
Data vs Draw.
double total_llh_data_drawfluc
Data vs Fluctuated Draw.
double total_llh_drawfluc_draw
Fluctuated Draw vs Draw.
std::vector< double > llh_draw_pred
Draw vs Predictive.
double total_llh_predfluc_draw
Fluctuated Predictive vs Draw.
std::vector< double > llh_drawfluc_draw_ProjectX
std::vector< double > llh_data_drawfluc
Data vs Fluctuated Draw.
double total_llh_drawfluc_pred
Fluctuated Draw vs Predictive.
std::vector< double > llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using rate only.
double total_llh_rate_data_draw
Rate Data vs Draw.
std::vector< double > llh_drawfluc_pred
Fluctuated Draw vs Predictive.
std::vector< double > llh_datafluc_draw
Fluctuated Data vs Draw.
void MakeFluctuatedHistogram(TH1D *FluctHist, TH1D *PolyHist)
Make Poisson fluctuation of TH1D hist.
std::vector< double > llh_rate_data_draw
Data vs Draw using rate only.
double total_llh_data_predfluc
Data vs Fluctuated Predictive.
double llh_penalty
LLH penalty for each throw.
double total_llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
double total_llh_datafluc_draw
Fluctuated Data vs Draw.
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:212

◆ MakeCutEventRate()

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

Make the 1D Event Rate Hist.

Definition at line 1511 of file SampleSummary.cpp.

1511 {
1512// ****************
1513 // For the event rate histogram add a TLine to the data rate
1514 auto TempLine = std::make_unique<TLine>(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum());
1515 TempLine->SetLineColor(kRed);
1516 TempLine->SetLineWidth(2);
1517 // Also fit a Gaussian because why not?
1518 TF1 *Fitter = new TF1("Fit", "gaus", Histogram->GetBinLowEdge(1), Histogram->GetBinLowEdge(Histogram->GetNbinsX()+1));
1519 Histogram->Fit(Fitter, "RQ");
1520 Fitter->SetLineColor(kRed-5);
1521 // Calculate a p-value
1522 double Above = 0.0;
1523 for (int z = 0; z < Histogram->GetNbinsX(); ++z) {
1524 const double xvalue = Histogram->GetBinCenter(z+1);
1525 if (xvalue >= DataRate) {
1526 Above += Histogram->GetBinContent(z+1);
1527 }
1528 }
1529 const double pvalue = Above/Histogram->Integral();
1530 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1531 Legend->SetFillColor(0);
1532 Legend->SetFillStyle(0);
1533 Legend->SetLineWidth(0);
1534 Legend->SetLineColor(0);
1535 Legend->AddEntry(TempLine.get(), Form("Data, %.0f, p-value=%.2f", DataRate, pvalue), "l");
1536 Legend->AddEntry(Histogram, Form("MC, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
1537 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.1f#pm%.1f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
1538 std::string TempTitle = std::string(Histogram->GetName());
1539 TempTitle += "_canv";
1540 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1541 TempCanvas->SetGridx();
1542 TempCanvas->SetGridy();
1543 TempCanvas->SetRightMargin(0.03);
1544 TempCanvas->SetBottomMargin(0.08);
1545 TempCanvas->SetLeftMargin(0.10);
1546 TempCanvas->SetTopMargin(0.06);
1547 TempCanvas->cd();
1548 Histogram->Draw();
1549 TempLine->Draw("same");
1550 Fitter->Draw("same");
1551 Legend->Draw("same");
1552 TempCanvas->Write();
1553 Histogram->Write();
1554
1555 delete TempCanvas;
1556 delete Fitter;
1557}

◆ MakeCutLLH()

void SampleSummary::MakeCutLLH ( )
inlineprivate

Make the cut LLH histogram.

Definition at line 1433 of file SampleSummary.cpp.

1433 {
1434// *******************
1435 Outputfile->cd();
1436 MakeCutLLH1D(lnLHist.get());
1440
1445}
void Get2DBayesianpValue(TH2D *Histogram)
Calculates the 2D Bayesian p-value and generates a visualization.
void MakeCutLLH1D(TH1D *Histogram, double llh_ref=-999)

◆ MakeCutLLH1D()

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

Definition at line 1449 of file SampleSummary.cpp.

1449 {
1450// ****************
1451 const double TotalIntegral = Histogram->Integral();
1452 double Above = 0.0;
1453 // Get the LLH reference from total llh or some reference histogram
1454 double llh_reference = 0.0;
1455 if (llh_ref >= 0) {
1456 llh_reference = llh_ref;
1457 } else {
1458 llh_reference = llh_total;
1459 }
1460 for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) {
1461 const double xvalue = Histogram->GetBinCenter(i+1);
1462 if (xvalue >= llh_reference) {
1463 Above += Histogram->GetBinContent(i+1);
1464 }
1465 }
1466 const double pvalue = Above/TotalIntegral;
1467 std::stringstream ss;
1468 ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue;
1469 Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str());
1470
1471 // Write a TCanvas and make a line and a filled histogram
1472 auto TempLine = std::make_unique<TLine>(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum());
1473 TempLine->SetLineColor(kBlack);
1474 TempLine->SetLineWidth(2);
1475
1476 // Make the fill histogram
1477 TH1D *TempHistogram = static_cast<TH1D*>(Histogram->Clone());
1478 TempHistogram->SetFillStyle(1001);
1479 TempHistogram->SetFillColor(kRed);
1480 for (int i = 0; i < TempHistogram->GetNbinsX(); ++i) {
1481 if (TempHistogram->GetBinCenter(i+1) < llh_reference) {
1482 TempHistogram->SetBinContent(i+1, 0.0);
1483 }
1484 }
1485
1486 auto Legend = std::make_unique<TLegend>(0.6, 0.6, 0.9, 0.9);
1487 Legend->SetFillColor(0);
1488 Legend->SetFillStyle(0);
1489 Legend->SetLineWidth(0);
1490 Legend->SetLineColor(0);
1491 Legend->AddEntry(TempLine.get(), Form("Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue), "l");
1492 Legend->AddEntry(Histogram, Form("LLH, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
1493 std::string Title = Histogram->GetName();
1494 Title += "_canv";
1495 TCanvas *TempCanvas = new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024);
1496 TempCanvas->SetGridx();
1497 TempCanvas->SetGridy();
1498 Histogram->Draw();
1499 TempHistogram->Draw("same");
1500 TempLine->Draw("same");
1501 Legend->Draw("same");
1502
1503 TempCanvas->Write();
1504
1505 delete TempHistogram;
1506 delete TempCanvas;
1507}
double llh_total
Total LLH for the posterior predictive distribution.

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

2183 {
2184// ****************
2185 if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
2186 else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
2187}
void MakeFluctuatedHistogramAlternative(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check.
void MakeFluctuatedHistogramStandard(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using default fast method.

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

2191 {
2192// ****************
2193 if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
2194 else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
2195}

◆ MakePredictive()

void SampleSummary::MakePredictive ( )
inlineprivate

Finalise the distributions from the thrown samples.

Definition at line 978 of file SampleSummary.cpp.

978 {
979// *******************
980 // First make the projection on the z axis of the TH3D* for every pmu cosmu bin
981 double llh_total_temp = 0.0;
982
983 // Loop over the samples
984 #ifdef MULTITHREAD
985 #pragma omp parallel for reduction(+:llh_total_temp)
986 #endif
987 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
988 {
989 // Skip disabled samples
990 if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
991
992 // Count the -2LLH for each histogram
993 double negLogL_Mean = 0.0;
994 double negLogL_Mode = 0.0;
995
996 // Loop over each pmu cosmu bin
997 for (int j = 1; j < maxBins[SampleNum]+1; ++j)
998 {
999 TH1D *Projection = PosteriorHist[SampleNum][j].get();
1000 TH1D *W2Projection = w2Hist[SampleNum][j].get();
1001
1002 // Data content for the j,kth bin
1003 const double nData = DataHist[SampleNum]->GetBinContent(j);
1004
1005 // Get the mean for this projection for all the samples
1006 // This is the mean prediction for this given j,k bin
1007 const double nMean = Projection->GetMean();
1008 const double nMeanError = Projection->GetRMS();
1009 const double nMode = Projection->GetBinCenter(Projection->GetMaximumBin());
1010 const double nModeError = GetModeError(Projection);
1011
1012 const double nW2Mean = W2Projection->GetMean();
1013 const double nW2Mode = W2Projection->GetBinCenter(W2Projection->GetMaximumBin());
1014
1015 double TempLLH_Mean = 0.0;
1016 double TempLLH_Mode = 0.0;
1017
1018 //KS:Get LLH contribution getTestStatLLH can calculate Barlow Beeston/IceCube or Poisson
1019 TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
1020 TempLLH_Mode = SampleHandler->GetTestStatLLH(nData, nMode, nW2Mode);
1021
1022 // Increment -2LLH
1023 //KS: do times 2 because banff reports chi2
1024 negLogL_Mean += 2*TempLLH_Mean;
1025 negLogL_Mode += 2*TempLLH_Mode;
1026
1027 // Set the content and error to the mean in the bin
1028 MeanHist[SampleNum]->SetBinContent(j, MeanHist[SampleNum]->GetBinContent(j)+nMean);
1029 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1030 MeanHist[SampleNum]->SetBinError(j, nMeanError);
1031
1032 if(DoBetaParam)
1033 {
1034 TH1D *BetaTemp = BetaHist[SampleNum][j].get();
1035 const double nBetaMean = BetaTemp->GetMean();
1036 const double nBetaMeanError = BetaTemp->GetRMS();
1037 //KS: Here we modify predictions by beta parameter from Barlow-Beeston
1038 MeanHistCorrected[SampleNum]->SetBinContent(j, MeanHistCorrected[SampleNum]->GetBinContent(j)+nMean*nBetaMean);
1039 //KS: Use error propagation to calcuate error
1040 const double ErrorTemp = std::sqrt( (nBetaMean*nMeanError) * (nBetaMean*nMeanError) + (nMean*nBetaMeanError) * (nMean*nBetaMeanError));
1041 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1042 MeanHistCorrected[SampleNum]->SetBinError(j, ErrorTemp);
1043 }
1044 // Set the content to the mode in the bin
1045 ModeHist[SampleNum]->SetBinContent(j, ModeHist[SampleNum]->GetBinContent(j)+nMode);
1046 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1047 ModeHist[SampleNum]->SetBinError(j, nModeError);
1048 // Set the content to the mean in the bin
1049 W2MeanHist[SampleNum]->SetBinContent(j, W2MeanHist[SampleNum]->GetBinContent(j)+nW2Mean);
1050 // Set the content to the mode in the bin
1051 W2ModeHist[SampleNum]->SetBinContent(j, W2ModeHist[SampleNum]->GetBinContent(j)+nW2Mode);
1052
1053 // Set the mean and average LLH for this given bin
1054 // Can use these hists to see where the largest -2LLH hists come from
1055 lnLHist_Mean[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
1056 lnLHist_Mode[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mode);
1057
1058 lnLHist_Mean1D[SampleNum]->Fill(2.0*TempLLH_Mean);
1059 lnLHist_Mode1D[SampleNum]->Fill(2.0*TempLLH_Mode);
1060 } // End loop over bins
1061 if(DoByModePlots)
1062 {
1063 for (int j = 0; j < Modes->GetNModes()+1; j++)
1064 {
1065 // Loop over each pmu cosmu bin
1066 for (int i = 1; i < maxBins[SampleNum]+1; ++i)
1067 {
1068 // Make the posterior/prior predictive projection on z
1069 // The z axis of Predictive is the bin content
1070 // Essentially zooming in on one bin and looking at the mean and mode of that bin
1071 TH1D *Projection = PosteriorHist_ByMode[SampleNum][j][i];
1072
1073 // Get the mean for this projection for all the samples
1074 const double nMean = Projection->GetMean();
1075 const double nMeanError = Projection->GetRMS();
1076
1077 // Set the content and error to the mean in the bin
1078 MeanHist_ByMode[SampleNum][j]->SetBinContent(i, MeanHist_ByMode[SampleNum][j]->GetBinContent(i)+nMean);
1079 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1080 MeanHist_ByMode[SampleNum][j]->SetBinError(i, nMeanError);
1081 } // End loop over bins
1082 }
1083 }
1084 llh_total_temp += negLogL_Mean;
1085 } // End loop over samples
1086
1087 // This is not multithreaded as due to ProjectPoly it is not safe
1088 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1089 {
1090 // Skip disabled samples
1091 if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
1092
1093 //KS:: Might consider caching it as we use it once agian much later
1094 TH1D *MeanProjectX = ProjectPoly(MeanHist[SampleNum], true, SampleNum, true);
1095 TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[SampleNum], true, SampleNum);
1096 // Loop over each pmu bin for 1D projection
1097 for (int j = 1; j <= lnLHist_Mean_ProjectX[SampleNum]->GetXaxis()->GetNbins(); ++j)
1098 {
1099 // Data content for the j,kth bin
1100 const double nData = DataHist_ProjectX[SampleNum]->GetBinContent(j);
1101 const double nMean = MeanProjectX->GetBinContent(j);
1102 const double nW2Mean = W2MeanProjectX->GetBinContent(j);
1103
1104 double TempLLH_Mean = 0.0;
1105 TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
1106
1107 //KS: do times 2 because banff reports chi2
1108 lnLHist_Mean_ProjectX[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
1109 }// End loop over bins
1110
1111 delete MeanProjectX;
1112 delete W2MeanProjectX;
1113 } // End loop over samples
1114
1115 llh_total = llh_total_temp;
1116 // Now we have our posterior predictive histogram and it's LLH
1117 MACH3LOG_INFO("Prior/Posterior predictive LLH mean (sample only) = {:.2f}", llh_total);
1118 std::stringstream ss;
1119 ss << llh_total;
1120 lnLHist->SetTitle((std::string(lnLHist->GetTitle())+"_"+ss.str()).c_str());
1121
1122 // Now make the fluctuated hists of the MeanHist and ModeHist
1123 MakeChi2Hists();
1124
1125 // Get the 1D LLH dists
1126 MakeCutLLH();
1127} // End MakePredictive() function
double GetModeError(TH1D *hpost)
Get the mode error from a TH1D.
void MakeChi2Hists()
Make the fluctuated histograms (2D and 1D) for the chi2s Essentially taking the MCMC draws and calcul...
void MakeCutLLH()
Make the cut LLH histogram.

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

1612 {
1613// ****************
1614 // Make a new directory
1615 TDirectory *BetaDir = Outputfile->mkdir("BetaParameters");
1616 BetaDir->cd();
1617
1618 int originalErrorLevel = gErrorIgnoreLevel;
1619
1620 //To avoid Warning in <Fit>: Fit data is empty
1621 gErrorIgnoreLevel = kFatal;
1622
1623 MACH3LOG_INFO("Writing Beta parameters");
1624 std::vector<TDirectory *> DirBeta(nSamples);
1625 for (int i = 0; i < nSamples; ++i)
1626 {
1627 // Make a new directory
1628 DirBeta[i] = BetaDir->mkdir((SampleNames[i]).c_str());
1629 DirBeta[i]->cd();
1630
1631 // Loop over each pmu cosmu bin
1632 for (int j = 1; j < maxBins[i]+1; ++j)
1633 {
1634 const double data = DataHist[i]->GetBinContent(j);
1635 const double mc = NominalHist[i]->GetBinContent(j);
1636 const double w2 = W2NomHist[i]->GetBinContent(j);
1637
1638 const double BetaPrior = GetBetaParameter(data, mc, w2, likelihood);
1639
1640 auto TempLine = std::unique_ptr<TLine>(new TLine(BetaPrior, BetaHist[i][j]->GetMinimum(), BetaPrior, BetaHist[i][j]->GetMaximum()));
1641 TempLine->SetLineColor(kRed);
1642 TempLine->SetLineWidth(2);
1643
1644 // Also fit a Gaussian because why not?
1645 TF1 *Fitter = new TF1("Fit", "gaus", BetaHist[i][j]->GetBinLowEdge(1), BetaHist[i][j]->GetBinLowEdge(BetaHist[i][j]->GetNbinsX()+1));
1646 BetaHist[i][j]->Fit(Fitter, "RQ");
1647 Fitter->SetLineColor(kRed-5);
1648
1649 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1650 Legend->SetFillColor(0);
1651 Legend->SetFillStyle(0);
1652 Legend->SetLineWidth(0);
1653 Legend->SetLineColor(0);
1654 Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.4f, N_{data}=%.0f", BetaPrior, data), "l");
1655 Legend->AddEntry(BetaHist[i][j].get(), Form("Post, #mu=%.4f#pm%.4f", BetaHist[i][j]->GetMean(), BetaHist[i][j]->GetRMS()), "l");
1656 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.4f#pm%.4f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
1657 std::string TempTitle = std::string(BetaHist[i][j]->GetName());
1658
1659 TempTitle += "_canv";
1660 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1661 TempCanvas->SetGridx();
1662 TempCanvas->SetGridy();
1663 TempCanvas->SetRightMargin(0.03);
1664 TempCanvas->SetBottomMargin(0.08);
1665 TempCanvas->SetLeftMargin(0.10);
1666 TempCanvas->SetTopMargin(0.06);
1667 TempCanvas->cd();
1668 BetaHist[i][j]->Draw();
1669 TempLine->Draw("same");
1670 Fitter->Draw("same");
1671 Legend->Draw("same");
1672 TempCanvas->Write();
1673 BetaHist[i][j]->Write();
1674
1675 delete TempCanvas;
1676 delete Fitter;
1677 }
1678 DirBeta[i]->Write();
1679 delete DirBeta[i];
1680 }
1681 BetaDir->Write();
1682 delete BetaDir;
1683
1684 gErrorIgnoreLevel = originalErrorLevel;
1685 Outputfile->cd();
1686}
std::vector< std::string > SampleNames
name for each sample

◆ PrepareOutput()

void SampleSummary::PrepareOutput ( )
inlineprivate

KS: Prepare output tree and necessary variables.

Definition at line 556 of file SampleSummary.cpp.

556 {
557// **********************
558
559 // Make the output file (MakePosterioPredictive call writes to this)
560 std::string TempString = OutputName;
561 TempString.replace(TempString.find(".root"), 5, std::string("_procsW2.root"));
562 Outputfile = new TFile(TempString.c_str(), "RECREATE");
563
564 // The array of doubles we write to the TTree
565 // Data vs Draw
566 llh_data_draw.resize(nSamples);
567 // Fluctuated Draw vs Draw
569 // Fluctuated Predicitve vs Draw
571
572 // Data vs Draw using Rate
574 // Data vs Fluctuated Predictive using Rate
576
577 // Data vs Fluctuated Draw
579 // Data vs Fluctuated Predictive
581 // Draw vs Predictive
582 llh_draw_pred.resize(nSamples);
583 // Fluctuated Draw vs Predictive
585 // Fluctuated Draw vs Fluctuated Predictive
587
588 // Fluctuated Predictive vs Predictive
590 // Fluctuated Data vs Draw
592
593 // Data vs Draw for 1D projection
596
597 // The output tree we're going to write to
598 OutputTree = new TTree("LLH_draws", "LLH_draws");
599 SampleNames.resize(nSamples);
600 // Loop over the samples and set the addresses of the variables to write to file
601 for (int i = 0; i < nSamples; ++i)
602 {
603 // Get the name
604 std::string SampleName = SampleHandler->GetSampleName(i);
605 // Strip out spaces
606 while (SampleName.find(" ") != std::string::npos) {
607 SampleName.replace(SampleName.find(" "), 1, std::string("_"));
608 }
609 SampleNames[i] = SampleName;
610 //CW: Also strip out - signs because it messes up TBranches
611 while (SampleName.find("-") != std::string::npos) {
612 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
613 }
614// All LLH below are used for actual p-value calculations
615 OutputTree->Branch((SampleName+"_data_draw").c_str(), &llh_data_draw[i]);
616 OutputTree->Branch((SampleName+"_drawfluc_draw").c_str(), &llh_drawfluc_draw[i]);
617 OutputTree->Branch((SampleName+"_predfluc_draw").c_str(), &llh_predfluc_draw[i]);
618
619// All LLH below are used for actual p-value calculations however using rate only
620 OutputTree->Branch((SampleName+"_rate_data_draw").c_str(), &llh_rate_data_draw[i]);
621 OutputTree->Branch((SampleName+"_rate_predfluc_draw").c_str(), &llh_rate_predfluc_draw[i]);
622
623// All LLH below are for validation reason but not used for final P-Value
624 OutputTree->Branch((SampleName+"_data_drawfluc").c_str(), &llh_data_drawfluc[i]);
625 OutputTree->Branch((SampleName+"_data_predfluc").c_str(), &llh_data_predfluc[i]);
626 OutputTree->Branch((SampleName+"_draw_pred").c_str(), &llh_draw_pred[i]);
627 OutputTree->Branch((SampleName+"_drawfluc_pred").c_str(), &llh_drawfluc_pred[i]);
628 OutputTree->Branch((SampleName+"_drawfluc_predfluc").c_str(), &llh_drawfluc_predfluc[i]);
629 OutputTree->Branch((SampleName+"_predfluc_pred").c_str(), &llh_predfluc_pred[i]);
630 OutputTree->Branch((SampleName+"_datafluc_draw").c_str(), &llh_datafluc_draw[i]);
631
632// All LLH below are used for calcauting P-Value but using 1D projections
633 OutputTree->Branch((SampleName+"_data_draw_ProjectX").c_str(), &llh_data_draw_ProjectX[i]);
634 OutputTree->Branch((SampleName+"_drawfluc_draw_ProjectX").c_str(), &llh_drawfluc_draw_ProjectX[i]);
635 }
636//All LLH below are used for actual p-value calculations
637 OutputTree->Branch("LLH_Penalty", &llh_penalty);
638 OutputTree->Branch("Total_LLH_Data_Draw", &total_llh_data_draw);
639 OutputTree->Branch("Total_LLH_DrawFluc_Draw", &total_llh_drawfluc_draw);
640 OutputTree->Branch("Total_LLH_PredFluc_Draw", &total_llh_predfluc_draw);
641
642// All LLH below are used for actual p-value calculations however using rate only
643 OutputTree->Branch("Total_LLH_Rate_PredFluc_Draw", &total_llh_rate_predfluc_draw);
644
645//All LLH below are for validation reason but not used for final P-Value
646 OutputTree->Branch("Total_LLH_Data_DrawFluc", &total_llh_data_drawfluc);
647 OutputTree->Branch("Total_LLH_Data_PredFluc", &total_llh_data_predfluc);
648 OutputTree->Branch("Total_LLH_Draw_Pred", &total_llh_draw_pred);
649 OutputTree->Branch("Total_LLH_DrawFluc_Pred", &total_llh_drawfluc_pred);
650 OutputTree->Branch("Total_LLH_DrawFluc_PredFluc", &total_llh_drawfluc_predfluc);
651 OutputTree->Branch("Total_LLH_PredFluc_Pred", &total_llh_predfluc_pred);
652 OutputTree->Branch("Total_LLH_DataFluc_Draw", &total_llh_datafluc_draw);
653
654//All LLH below are used for calcauting P-Value but 1D projections
655 OutputTree->Branch("total_llh_data_draw_ProjectX", &total_llh_data_draw_ProjectX);
656 OutputTree->Branch("total_llh_drawfluc_draw_ProjectX", &total_llh_drawfluc_draw_ProjectX);
657
658 Outputfile->cd();
659 Dir.resize(nSamples);
660 for (int i = 0; i < nSamples; ++i)
661 {
662 // Make a new directory
663 Dir[i] = Outputfile->mkdir((SampleNames[i]).c_str());
664 }
665}
std::vector< TDirectory * > Dir
Directory for each sample.

◆ ProjectHist()

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

Helper to project TH2D onto axis.

Definition at line 2147 of file SampleSummary.cpp.

2147 {
2148// ****************
2149 TH1D* Projection = nullptr;
2150 std::string name;
2151 if (ProjectX) {
2152 name = std::string(Histogram->GetName()) + "_x";
2153 Projection = Histogram->ProjectionX(name.c_str(), 1, Histogram->GetYaxis()->GetNbins(), "e");
2154 } else {
2155 name = std::string(Histogram->GetName()) + "_y";
2156 Projection = Histogram->ProjectionY(name.c_str(), 1, Histogram->GetXaxis()->GetNbins(), "e");
2157 }
2158 return Projection;
2159}

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

2163 {
2164// ****************
2165 std::vector<double> xbins;
2166 std::vector<double> ybins;
2167
2168 SampleHandler->SetupBinning(M3::int_t(selection), xbins, ybins);
2169 TH1D* Projection = nullptr;
2170 std::string name;
2171 if (ProjectX) {
2172 name = std::string(Histogram->GetName()) + "_x";
2173 Projection = PolyProjectionX(Histogram, name.c_str(), xbins, MakeErrorHist);
2174 } else {
2175 name = std::string(Histogram->GetName()) + "_y";
2176 Projection = PolyProjectionY(Histogram, name.c_str(), ybins, MakeErrorHist);
2177 }
2178 return Projection;
2179}
TH1D * PolyProjectionX(TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
WP: Poly Projectors.
TH1D * PolyProjectionY(TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
WP: Poly Projectors.

◆ SetLikelihood()

void SampleSummary::SetLikelihood ( const TestStatistic  TestStat)
inline

KS: Set likelihood type.

Definition at line 48 of file SampleSummary.h.

48{ 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.

50{ nModelParams = nPars;};

◆ StudyBIC()

void SampleSummary::StudyBIC ( )
inlineprivate

Study Bayesian Information Criterion (BIC) [11].

Definition at line 2227 of file SampleSummary.cpp.

2227 {
2228// ****************
2229 //make fancy event rate histogram
2230 double DataRate = 0.0;
2231 double BinsRate = 0.0;
2232 #ifdef MULTITHREAD
2233 #pragma omp parallel for reduction(+:DataRate, BinsRate)
2234 #endif
2235 for (int i = 0; i < nSamples; ++i)
2236 {
2237 if (DataHist[i] == nullptr) continue;
2238 DataRate += NoOverflowIntegral(DataHist[i]);
2239 BinsRate += maxBins[i];
2240 }
2241
2242 const double EventRateBIC = GetBIC(llh_total, DataRate, nModelParams);
2243 const double BinBasedBIC = GetBIC(llh_total, BinsRate, nModelParams);
2244 MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of events: {:.2f}", EventRateBIC);
2245 MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of bins: {:.2f}", BinBasedBIC);
2246 MACH3LOG_INFO("Additional info: nModelParams {} DataRate: {:.2f} BinsRate: {:.2f}", nModelParams, DataRate, BinsRate);
2247}
double GetBIC(const double llh, const int data, const int nPars)
Get the Bayesian Information Criterion (BIC) or Schwarz information criterion (also SIC,...

◆ StudyDIC()

void SampleSummary::StudyDIC ( )
inlineprivate

KS: Get the Deviance Information Criterion (DIC) [26] [28].

Definition at line 2251 of file SampleSummary.cpp.

2251 {
2252// ****************
2253 //The posterior mean of the deviance
2254 double Dbar = 0.;
2255
2256 #ifdef MULTITHREAD
2257 #pragma omp parallel for reduction(+:Dbar)
2258 #endif
2259 for (unsigned int i = 0; i < nThrows; ++i)
2260 {
2261 double LLH_temp = 0.;
2262 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
2263 {
2264 // Get -2*log-likelihood
2265 LLH_temp += GetLLH(DataHist[SampleNum], MCVector[i][SampleNum], W2MCVector[i][SampleNum]);
2266 }
2267 Dbar += LLH_temp;
2268 }
2269 Dbar = Dbar / nThrows;
2270
2271 // A point estimate of the deviance
2272 const double Dhat = llh_total;
2273
2274 //Effective number of parameters
2275 const double p_D = std::fabs(Dbar - Dhat);
2276
2277 //Actual test stat
2278 const double DIC_stat = Dhat + 2 * p_D;
2279 MACH3LOG_INFO("Effective number of parameters following DIC formalism is equal to: {:.2f}", p_D);
2280 MACH3LOG_INFO("DIC test statistic = {:.2f}", DIC_stat);
2281}

◆ StudyInformationCriterion()

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

Information Criterion.

Definition at line 2199 of file SampleSummary.cpp.

2199 {
2200// ****************
2201 MACH3LOG_INFO("******************************");
2202 switch(Criterion) {
2203 case M3::kInfCrit::kBIC:
2204 // Study Bayesian Information Criterion
2205 StudyBIC();
2206 break;
2207 case M3::kInfCrit::kDIC:
2208 // Study Deviance Information Criterion
2209 StudyDIC();
2210 break;
2212 // Study Watanabe-Akaike information criterion (WAIC)
2213 StudyWAIC();
2214 break;
2216 MACH3LOG_ERROR("kInfCrits is not a valid kInfCrit!");
2217 throw MaCh3Exception(__FILE__, __LINE__);
2218 default:
2219 MACH3LOG_ERROR("UNKNOWN Information Criterion SPECIFIED!");
2220 MACH3LOG_ERROR("You gave {}", static_cast<int>(Criterion));
2221 throw MaCh3Exception(__FILE__ , __LINE__ );
2222 }
2223 MACH3LOG_INFO("******************************");
2224}
void StudyBIC()
Study Bayesian Information Criterion (BIC) .
void StudyDIC()
KS: Get the Deviance Information Criterion (DIC) .
void StudyWAIC()
KS: Get the Watanabe-Akaike information criterion (WAIC) .
@ kWAIC
Watanabe-Akaike information criterion.
Definition: SampleSummary.h:13
@ kInfCrits
This only enumerates.
Definition: SampleSummary.h:14
@ kBIC
Bayesian Information Criterion.
Definition: SampleSummary.h:11
@ kDIC
Deviance Information Criterion.
Definition: SampleSummary.h:12

◆ StudyKinematicCorrelations()

void SampleSummary::StudyKinematicCorrelations ( )
inlineprivate

KS: Study how correlated are sample or kinematic bins.

Definition at line 1690 of file SampleSummary.cpp.

1690 {
1691// ****************
1692 MACH3LOG_INFO("Calculating Correlations");
1693 TStopwatch timer;
1694 timer.Start();
1695
1696 // Data vs Draw for 1D projection
1697 std::vector<double> NEvents_Sample(nSamples);
1698 double event_rate = 0.;
1699
1700 // The output tree we're going to write to
1701 TTree* Event_Rate_Tree = new TTree("Event_Rate_draws", "Event_Rate_draws");
1702 Event_Rate_Tree->Branch("Event_Rate", &event_rate);
1703 // Loop over the samples and set the addresses of the variables to write to file
1704 for (int i = 0; i < nSamples; ++i)
1705 {
1706 // Get the name
1707 std::string SampleName = SampleNames[i];
1708 //CW: Also strip out - signs because it messes up TBranches
1709 while (SampleName.find("-") != std::string::npos) {
1710 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
1711 }
1712 Event_Rate_Tree->Branch((SampleName+"_Event_Rate").c_str(), &NEvents_Sample[i]);
1713 }
1714
1715 // Holds the total event rate
1716 auto EventHist = std::make_unique<TH1D>("EventHist", "Total Event Rate", 100, 1, -1);
1717 EventHist->SetDirectory(nullptr);
1718 EventHist->GetXaxis()->SetTitle("Total event rate");
1719 EventHist->GetYaxis()->SetTitle("Counts");
1720 EventHist->SetLineWidth(2);
1721
1722 // Holds the event rate for the distribution
1723 std::vector<std::unique_ptr<TH1D>> SumHist(nSamples);
1724 for (int i = 0; i < nSamples; ++i)
1725 {
1726 std::string name = std::string(NominalHist[i]->GetName());
1727 name = name.substr(0, name.find("_nom"));
1728
1729 SumHist[i] = std::make_unique<TH1D>((name+"_sum").c_str(),(name+"_sum").c_str(), 100, 1, -1);
1730 SumHist[i]->GetXaxis()->SetTitle("N_{events}");
1731 SumHist[i]->GetYaxis()->SetTitle("Counts");
1732 double Integral = NoOverflowIntegral(DataHist[i]);
1733 std::stringstream ss;
1734 ss << Integral;
1735 SumHist[i]->SetTitle((std::string(SumHist[i]->GetTitle())+"_"+ss.str()).c_str());
1736 }
1737
1738 for (unsigned int it = 0; it < nThrows; ++it)
1739 {
1740 double event_rate_temp = 0.;
1741 // Loop over the samples
1742 #ifdef MULTITHREAD
1743 #pragma omp parallel for reduction(+:event_rate_temp)
1744 #endif
1745 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1746 {
1747 NEvents_Sample[SampleNum] = NoOverflowIntegral(MCVector[it][SampleNum]);
1748 // Fill the sum histogram with the integral of the sampled distribution
1749 SumHist[SampleNum]->Fill(NEvents_Sample[SampleNum], WeightVector[it]);
1750
1751 event_rate_temp += NEvents_Sample[SampleNum];
1752 } // end samples loop
1753 event_rate = event_rate_temp;
1754 EventHist->Fill(event_rate);
1755 Event_Rate_Tree->Fill();
1756 } //end loops over throws
1757 Event_Rate_Tree->Write();
1758 delete Event_Rate_Tree;
1759
1760 double DataRate = 0.0;
1761 #ifdef MULTITHREAD
1762 #pragma omp parallel for reduction(+:DataRate)
1763 #endif
1764 for (int i = 0; i < nSamples; ++i)
1765 {
1766 DataRate += NoOverflowIntegral(DataHist[i]);
1767 }
1768 MakeCutEventRate(EventHist.get(), DataRate);
1769
1770 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1771 {
1772 Dir[SampleNum]->cd();
1773 //Make fancy event rate histogram
1774 MakeCutEventRate(SumHist[SampleNum].get(), NoOverflowIntegral(DataHist[SampleNum]));
1775 }
1776
1777 // Make a new directory
1778 TDirectory *CorrDir = Outputfile->mkdir("Correlations");
1779 CorrDir->cd();
1780
1781 TMatrixDSym* SampleCorrelation = new TMatrixDSym(nSamples);
1782 std::vector<std::vector<std::unique_ptr<TH2D>>> SamCorr(nSamples);
1783 for (int i = 0; i < nSamples; ++i)
1784 {
1785 SamCorr[i].resize(nSamples);
1786
1787 (*SampleCorrelation)(i,i) = 1.0;
1788 const double Min_i = SumHist[i]->GetXaxis()->GetBinLowEdge(1);
1789 const double Max_i = SumHist[i]->GetXaxis()->GetBinUpEdge(SumHist[i]->GetNbinsX()+1);
1790 for (int j = 0; j < nSamples; ++j)
1791 {
1792 const double Min_j = SumHist[j]->GetXaxis()->GetBinLowEdge(1);
1793 const double Max_j = SumHist[j]->GetXaxis()->GetBinUpEdge(SumHist[j]->GetNbinsX()+1);
1794
1795 // TH2D to hold the Correlation
1796 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);
1797 SamCorr[i][j]->SetDirectory(nullptr);
1798 SamCorr[i][j]->SetMinimum(0);
1799 SamCorr[i][j]->GetXaxis()->SetTitle(SampleNames[i].c_str());
1800 SamCorr[i][j]->GetYaxis()->SetTitle(SampleNames[j].c_str());
1801 SamCorr[i][j]->GetZaxis()->SetTitle("Events");
1802 }
1803 }
1804
1805 // Now we are sure we have the diagonal elements, let's make the off-diagonals
1806 #ifdef MULTITHREAD
1807 #pragma omp parallel for
1808 #endif
1809 for (int i = 0; i < nSamples; ++i)
1810 {
1811 for (int j = 0; j <= i; ++j)
1812 {
1813 // Skip the diagonal elements which we've already done above
1814 if (j == i) continue;
1815
1816 for (unsigned int it = 0; it < nThrows; ++it)
1817 {
1818 SamCorr[i][j]->Fill(NoOverflowIntegral(MCVector[it][i]), NoOverflowIntegral(MCVector[it][j]));
1819 }
1820 SamCorr[i][j]->Smooth();
1821
1822 // Get the Covariance for these two parameters
1823 (*SampleCorrelation)(i,j) = SamCorr[i][j]->GetCorrelationFactor();
1824 (*SampleCorrelation)(j,i) = (*SampleCorrelation)(i,j);
1825 }// End j loop
1826 }// End i loop
1827
1828 auto hSamCorr = std::make_unique<TH2D>("Sample Correlation", "Sample Correlation", nSamples, 0, nSamples, nSamples, 0, nSamples);
1829 hSamCorr->SetDirectory(nullptr);
1830 hSamCorr->GetZaxis()->SetTitle("Correlation");
1831 hSamCorr->SetMinimum(-1);
1832 hSamCorr->SetMaximum(1);
1833 hSamCorr->GetXaxis()->SetLabelSize(0.015);
1834 hSamCorr->GetYaxis()->SetLabelSize(0.015);
1835
1836 // Loop over the Covariance matrix entries
1837 for (int i = 0; i < nSamples; ++i)
1838 {
1839 hSamCorr->GetXaxis()->SetBinLabel(i+1, SampleNames[i].c_str());
1840
1841 for (int j = 0; j < nSamples; ++j)
1842 {
1843 hSamCorr->GetYaxis()->SetBinLabel(j+1, SampleNames[j].c_str());
1844 // The value of the Covariance
1845 const double corr = (*SampleCorrelation)(i,j);
1846 hSamCorr->SetBinContent(i+1, j+1, corr);
1847 }
1848 }
1849 hSamCorr->Draw("colz");
1850 hSamCorr->Write("Sample_Corr");
1851
1852 SampleCorrelation->Write("Sample_Correlation");
1853 delete SampleCorrelation;
1854
1855 for (int i = 0; i < nSamples; ++i)
1856 {
1857 for (int j = 0; j <= i; ++j)
1858 {
1859 // Skip the diagonal elements which we've already done above
1860 if (j == i) continue;
1861 SamCorr[i][j]->Write();
1862 }// End j loop
1863 }// End i loop
1864
1865 //KS: This can take ages so better turn it off by default
1866 bool DoPerKinemBin = false;
1867 if(DoPerKinemBin)
1868 {
1869 //KS: Now the same but for kinematic bin of each sample
1870 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1871 {
1872 TMatrixDSym* KinCorrelation = new TMatrixDSym(maxBins[SampleNum]);
1873 std::vector<std::vector<std::unique_ptr<TH2D>>> KinCorr(maxBins[SampleNum]);
1874 for (int i = 0; i < maxBins[SampleNum]; ++i)
1875 {
1876 KinCorr[i].resize(maxBins[SampleNum]);
1877 (*KinCorrelation)(i,i) = 1.0;
1878
1879 const double Min_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinLowEdge(1);
1880 const double Max_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][i+1]->GetNbinsX()+1);
1881
1882 //Get PolyBin
1883 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
1884 // Just make a little fancy name
1885 std::stringstream ss2;
1886 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1887 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1888
1889 for (int j = 0; j < maxBins[SampleNum]; ++j)
1890 {
1891 const double Min_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinLowEdge(1);
1892 const double Max_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][j+1]->GetNbinsX()+1);
1893
1894 // TH2D to hold the Correlation
1895 KinCorr[i][j] = std::make_unique<TH2D>( Form("Kin_%i_%i_%i", SampleNum, i, j),
1896 Form("Kin_%i_%i_%i", SampleNum, i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
1897 KinCorr[i][j]->SetDirectory(nullptr);
1898 KinCorr[i][j]->SetMinimum(0);
1899
1900 KinCorr[i][j]->GetXaxis()->SetTitle(ss2.str().c_str());
1901
1902 bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
1903 // Just make a little fancy name
1904 std::stringstream ss3;
1905 ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1906 ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1907 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1908 KinCorr[i][j]->GetZaxis()->SetTitle("Events");
1909 }
1910 }
1911 // Now we are sure we have the diagonal elements, let's make the off-diagonals
1912 #ifdef MULTITHREAD
1913 #pragma omp parallel for
1914 #endif
1915 for (int i = 0; i < maxBins[SampleNum]; ++i)
1916 {
1917 for (int j = 0; j <= i; ++j)
1918 {
1919 // Skip the diagonal elements which we've already done above
1920 if (j == i) continue;
1921
1922 for (unsigned int it = 0; it < nThrows; ++it)
1923 {
1924 KinCorr[i][j]->Fill(MCVector[it][SampleNum]->GetBinContent(i+1), MCVector[it][SampleNum]->GetBinContent(j+1));
1925 }
1926 KinCorr[i][j]->Smooth();
1927
1928 // Get the Covariance for these two parameters
1929 (*KinCorrelation)(i,j) = KinCorr[i][j]->GetCorrelationFactor();
1930 (*KinCorrelation)(j,i) = (*KinCorrelation)(i,j);
1931 }// End j loop
1932 }// End i loop
1933
1934 auto hKinCorr = std::make_unique<TH2D>(SampleNames[SampleNum].c_str(), SampleNames[SampleNum].c_str(),
1935 maxBins[SampleNum], 0, maxBins[SampleNum], maxBins[SampleNum], 0, maxBins[SampleNum]);
1936 hKinCorr->SetDirectory(nullptr);
1937 hKinCorr->GetZaxis()->SetTitle("Correlation");
1938 hKinCorr->SetMinimum(-1);
1939 hKinCorr->SetMaximum(1);
1940 hKinCorr->GetXaxis()->SetLabelSize(0.015);
1941 hKinCorr->GetYaxis()->SetLabelSize(0.015);
1942
1943 // Loop over the Covariance matrix entries
1944 for (int i = 0; i < maxBins[SampleNum]; ++i)
1945 {
1946 //Get PolyBin
1947 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
1948 // Just make a little fancy name
1949 std::stringstream ss2;
1950 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1951 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1952 hKinCorr->GetXaxis()->SetBinLabel(i+1, ss2.str().c_str());
1953
1954 for (int j = 0; j < maxBins[SampleNum]; ++j)
1955 {
1956 bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
1957 // Just make a little fancy name
1958 std::stringstream ss3;
1959 ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1960 ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1961 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1962
1963 hKinCorr->GetYaxis()->SetBinLabel(j+1, ss3.str().c_str());
1964 // The value of the Covariance
1965 const double corr = (*KinCorrelation)(i,j);
1966 hKinCorr->SetBinContent(i+1, j+1, corr);
1967 }
1968 }
1969 hKinCorr->Draw("colz");
1970 hKinCorr->Write((SampleNames[SampleNum] + "_Corr").c_str());
1971
1972 KinCorrelation->Write((SampleNames[SampleNum] + "_Correlation").c_str());
1973 delete KinCorrelation;
1974
1975 /*
1976 for (int i = 0; i < maxBins[SampleNum]; ++i)
1977 {
1978 for (int j = 0; j <= i; ++j)
1979 {
1980 // Skip the diagonal elements which we've already done above
1981 if (j == i) continue;
1982 KinCorr[i][j]->Write();
1983 }// End j loop
1984 }// End i loop
1985 */
1986 }//end loop over samples
1987 }//end if DoPerKinemBin
1988 else
1989 {
1990 MACH3LOG_INFO("Not calculating correlations per each kinematic bin");
1991 }
1992
1993 if(DoByModePlots)
1994 {
1995 // Holds the total event rate by mode
1996 std::vector<TH1D*> EventHist_ByMode(Modes->GetNModes()+1);
1997 for (int j = 0; j < Modes->GetNModes()+1; j++)
1998 {
1999 std::string ModeName = Modes->GetMaCh3ModeName(j);
2000 EventHist_ByMode[j] = new TH1D(Form("EventHist_%s", ModeName.c_str()), Form("Total Event Rate %s", ModeName.c_str()), 100, 1, -1);
2001 EventHist_ByMode[j]->GetXaxis()->SetTitle("Total event rate");
2002 EventHist_ByMode[j]->GetYaxis()->SetTitle("Counts");
2003 EventHist_ByMode[j]->SetLineWidth(2);
2004 }
2005
2006 //KS: Here we calculate total event rates for each mode, maybe not most efficient but can be improved in the future
2007 for (unsigned int it = 0; it < nThrows; ++it)
2008 {
2009 for (int j = 0; j < Modes->GetNModes()+1; j++)
2010 {
2011 double event_rate_temp = 0.;
2012 #ifdef MULTITHREAD
2013 #pragma omp parallel for reduction(+:event_rate_temp)
2014 #endif
2015 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
2016 {
2017 event_rate_temp += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
2018 }
2019 EventHist_ByMode[j]->Fill(event_rate_temp);
2020 }
2021 }
2022
2023 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2024 {
2025 MakeCutEventRate(EventHist_ByMode[i], DataRate);
2026 }
2027
2028 TMatrixDSym* ModeCorrelation = new TMatrixDSym(Modes->GetNModes()+1);
2029
2030 TH2D*** ModeCorr = new TH2D**[Modes->GetNModes()+1]();
2031 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2032 {
2033 ModeCorr[i] = new TH2D*[Modes->GetNModes()+1]();
2034
2035 (*ModeCorrelation)(i,i) = 1.0;
2036
2037 const double Min_i = EventHist_ByMode[i]->GetXaxis()->GetBinLowEdge(1);
2038 const double Max_i = EventHist_ByMode[i]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[i]->GetNbinsX()+1);
2039 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2040 {
2041 const double Min_j = EventHist_ByMode[j]->GetXaxis()->GetBinLowEdge(1);
2042 const double Max_j = EventHist_ByMode[j]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[j]->GetNbinsX()+1);
2043
2044 // TH2D to hold the Correlation
2045 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);
2046 ModeCorr[i][j]->SetDirectory(nullptr);
2047 ModeCorr[i][j]->SetMinimum(0);
2048 ModeCorr[i][j]->GetXaxis()->SetTitle(Modes->GetMaCh3ModeName(i).c_str());
2049 ModeCorr[i][j]->GetYaxis()->SetTitle(Modes->GetMaCh3ModeName(j).c_str());
2050 ModeCorr[i][j]->GetZaxis()->SetTitle("Events");
2051 }
2052 }
2053
2054 // Now we are sure we have the diagonal elements, let's make the off-diagonals
2055 #ifdef MULTITHREAD
2056 #pragma omp parallel for
2057 #endif
2058 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2059 {
2060 for (int j = 0; j <= i; ++j)
2061 {
2062 // Skip the diagonal elements which we've already done above
2063 if (j == i) continue;
2064
2065 for (unsigned int it = 0; it < nThrows; ++it)
2066 {
2067 double Integral_X = 0.;
2068 double Integral_Y = 0.;
2069 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
2070 {
2071 Integral_X += NoOverflowIntegral(MCVectorByMode[it][SampleNum][i]);
2072 Integral_Y += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
2073 }
2074 ModeCorr[i][j]->Fill(Integral_X, Integral_Y);
2075 }
2076 ModeCorr[i][j]->Smooth();
2077
2078 // Get the Covariance for these two parameters
2079 (*ModeCorrelation)(i,j) = ModeCorr[i][j]->GetCorrelationFactor();
2080 (*ModeCorrelation)(j,i) = (*ModeCorrelation)(i,j);
2081 }// End j loop
2082 }// End i loop
2083
2084 TH2D* hModeCorr = new TH2D("Mode Correlation", "Mode Correlation", Modes->GetNModes()+1, 0, Modes->GetNModes()+1, Modes->GetNModes()+1, 0, Modes->GetNModes()+1);
2085 hModeCorr->SetDirectory(nullptr);
2086 hModeCorr->GetZaxis()->SetTitle("Correlation");
2087 hModeCorr->SetMinimum(-1);
2088 hModeCorr->SetMaximum(1);
2089 hModeCorr->GetXaxis()->SetLabelSize(0.015);
2090 hModeCorr->GetYaxis()->SetLabelSize(0.015);
2091
2092 // Loop over the Covariance matrix entries
2093 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2094 {
2095 hModeCorr->GetXaxis()->SetBinLabel(i+1, Modes->GetMaCh3ModeName(i).c_str());
2096
2097 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2098 {
2099 hModeCorr->GetYaxis()->SetBinLabel(j+1, Modes->GetMaCh3ModeName(j).c_str());
2100 // The value of the Covariance
2101 const double corr = (*ModeCorrelation)(i,j);
2102 hModeCorr->SetBinContent(i+1, j+1, corr);
2103 }
2104 }
2105 hModeCorr->Draw("colz");
2106 hModeCorr->Write("Mode_Corr");
2107 delete hModeCorr;
2108
2109 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2110 {
2111 for (int j = 0; j <= i; ++j)
2112 {
2113 // Skip the diagonal elements which we've already done above
2114 if (j == i) continue;
2115 ModeCorr[i][j]->Write();
2116 }// End j loop
2117 }// End i loop
2118
2119 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2120 {
2121 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2122 {
2123 delete ModeCorr[i][j];
2124 }
2125 delete[] ModeCorr[i];
2126 }
2127 delete[] ModeCorr;
2128 ModeCorrelation->Write("Mode_Correlation");
2129 delete ModeCorrelation;
2130
2131 for (int j = 0; j < Modes->GetNModes()+1; j++)
2132 {
2133 delete EventHist_ByMode[j];
2134 }
2135 }
2136
2137 CorrDir->Close();
2138 delete CorrDir;
2139
2140 timer.Stop();
2141 MACH3LOG_INFO("Calculating correlations took {:.2f}s", timer.RealTime());
2142 Outputfile->cd();
2143}
void MakeCutEventRate(TH1D *Histogram, const double DataRate)
Make the 1D Event Rate Hist.

◆ StudyWAIC()

void SampleSummary::StudyWAIC ( )
inlineprivate

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

Definition at line 2285 of file SampleSummary.cpp.

2285 {
2286// ****************
2287 // log pointwise predictive density
2288 double lppd = 0.;
2289 // effective number of parameters
2290 double p_WAIC = 0.;
2291
2292 #ifdef MULTITHREAD
2293 #pragma omp parallel for reduction(+:lppd, p_WAIC)
2294 #endif
2295 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum) {
2296 int nBins = maxBins[SampleNum];
2297 for (int i = 1; i <= nBins; ++i) {
2298 double mean_llh = 0.;
2299 double sum_exp_llh = 0;
2300 double mean_llh_squared = 0.;
2301
2302 for (unsigned int s = 0; s < nThrows; ++s) {
2303 const double data = DataHist[SampleNum]->GetBinContent(i);
2304 const double mc = MCVector[s][SampleNum]->GetBinContent(i);
2305 const double w2 = W2MCVector[s][SampleNum]->GetBinContent(i);
2306
2307 // Get the -log-likelihood for this sample and bin
2308 double neg_LLH_temp = SampleHandler->GetTestStatLLH(data, mc, w2);
2309
2310 // Negate the negative log-likelihood to get the actual log-likelihood
2311 double LLH_temp = -neg_LLH_temp;
2312
2313 mean_llh += LLH_temp;
2314 mean_llh_squared += LLH_temp * LLH_temp;
2315 sum_exp_llh += std::exp(LLH_temp);
2316 }
2317
2318 // Compute the mean log-likelihood and the squared mean
2319 mean_llh /= nThrows;
2320 mean_llh_squared /= nThrows;
2321 sum_exp_llh /= nThrows;
2322 sum_exp_llh = std::log(sum_exp_llh);
2323
2324 // Log pointwise predictive density based on Eq. 4 in Gelman2014
2325 lppd += sum_exp_llh;
2326
2327 // Compute the effective number of parameters for WAIC
2328 p_WAIC += mean_llh_squared - (mean_llh * mean_llh);
2329 }
2330 }
2331
2332 // Compute WAIC, see Eq. 13 in Gelman2014
2333 double WAIC = -2 * (lppd - p_WAIC);
2334 MACH3LOG_INFO("Effective number of parameters following WAIC formalism is equal to: {:.2f}", p_WAIC);
2335 MACH3LOG_INFO("WAIC = {:.2f}", WAIC);
2336}

◆ Write()

void SampleSummary::Write ( )

KS: Write results into root file.

Definition at line 669 of file SampleSummary.cpp.

669 {
670// *******************
671 // Prepare the output tree
673
674 MACH3LOG_INFO("Summarising {} throws...", nThrows);
675 // After all the throws are added finalise the sample
676 TStopwatch timer;
677 timer.Start();
679 timer.Stop();
680 MACH3LOG_INFO("Made Prior/Posterior Predictive, it took {:.2f}s, now writing...", timer.RealTime());
681
682 // Studying information criterion
684
685 OutputTree->Write();
686
687 // Make the various distributions
688 lnLHist->Write();
689 lnLHist_drawfluc->Write();
690 lnLHist_drawflucdraw->Write();
691 lnLHist_drawdata->Write();
692 lnLDrawHist->Write();
693 lnLFlucHist->Write();
694 lnLDrawHistRate->Write();
695 //KS: Only available for Posterior Predictive
696 if(!isPriorPredictive) RandomHist->Write();
697
698 lnLFlucHist_ProjectX->Write();
699
700 // Loop over each sample and write to file
701 //KS: Multithreading is tempting here but we also write to ROOT file, separating all LLH and poly projections from write could work well
702 for (int i = 0; i < nSamples; ++i)
703 {
704 // Skip the null histograms
705 if (DataHist[i] == nullptr || NoOverflowIntegral(DataHist[i]) == 0) continue;
706 Dir[i]->cd();
707
708 // Make the data/MC ratio histogram
709 TH2Poly *RatioHistMean = RatioPolys(DataHist[i], MeanHist[i]);
710 RatioHistMean->GetZaxis()->SetTitle("Data/Mean");
711 TH2Poly *RatioHistMode = RatioPolys(DataHist[i], ModeHist[i]);
712 RatioHistMode->GetZaxis()->SetTitle("Data/Mode");
713 TH2Poly *RatioHistNom = RatioPolys(DataHist[i], NominalHist[i]);
714 RatioHistNom->GetZaxis()->SetTitle("Data/Nom");
715
716 // And the normalised data histogram
717 TH2Poly *DataNormHist = NormalisePoly(DataHist[i]);
718 // Last true refers to if project along x or y
719 TH2Poly *MeanNormHist = NormalisePoly(MeanHist[i]);
720 TH2Poly *ModeNormHist = NormalisePoly(ModeHist[i]);
721 TH1D *MeanProjectX = ProjectPoly(MeanHist[i], true, i, true);
722 TH1D *MeanProjectY = ProjectPoly(MeanHist[i], false, i, true);
723 TH1D *ModeProjectX = ProjectPoly(ModeHist[i], true, i, true);
724 TH1D *ModeProjectY = ProjectPoly(ModeHist[i], false, i, true);
725
726 TH1D *MeanHistCorrectedProjectX = nullptr;
727 if(DoBetaParam) MeanHistCorrectedProjectX = ProjectPoly(MeanHistCorrected[i], true, i, true);
728 TH1D *MeanHistCorrectedProjectY = nullptr;
729 if(DoBetaParam) MeanHistCorrectedProjectY = ProjectPoly(MeanHistCorrected[i], false, i, true);
730
731 TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[i], true, i);
732 TH1D *W2MeanProjectY = ProjectPoly(W2MeanHist[i], false, i);
733 TH1D *W2ModeProjectX = ProjectPoly(W2ModeHist[i], true, i);
734 TH1D *W2ModeProjectY = ProjectPoly(W2ModeHist[i], false, i);
735
736 TH2Poly *NomNormHist = NormalisePoly(NominalHist[i]);
737 TH1D *NomProjectX = ProjectPoly(NominalHist[i], true, i);
738 TH1D *NomProjectY = ProjectPoly(NominalHist[i], false, i);
739
740 TH1D *W2NomProjectX = ProjectPoly(W2NomHist[i], true, i);
741 TH1D *W2NomProjectY = ProjectPoly(W2NomHist[i], false, i);
742
743 // Same for the TH2Ds
747
748 // Calculate the log likelihood for the 1D dists
749 // Sets the title of the second TH1D to the -2LLH
750 CalcLLH(DataHist_ProjectX[i], NomProjectX, W2NomProjectX);
751 CalcLLH(DataHist_ProjectX[i], MeanProjectX, W2MeanProjectX);
752 CalcLLH(DataHist_ProjectX[i], ModeProjectX, W2ModeProjectX);
753 CalcLLH(DataHist_ProjectY[i], NomProjectY, W2NomProjectY);
754 CalcLLH(DataHist_ProjectY[i], MeanProjectY, W2MeanProjectY);
755 CalcLLH(DataHist_ProjectY[i], ModeProjectY, W2ModeProjectY);
756
757 std::string SampleName = SampleNames[i];
758 // Also strip out - signs because it messes up TBranches
759 while (SampleName.find("-") != std::string::npos) {
760 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
761 }
762 OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_drawfluc_draw>>htemp").c_str());
763 TH2D *TempHistogram = static_cast<TH2D*>(gDirectory->Get("htemp")->Clone());
764 TempHistogram->GetXaxis()->SetTitle("-2LLH(Draw Fluc, Draw)");
765 TempHistogram->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
766 TempHistogram->SetNameTitle((SampleNames[i]+"_drawfluc_draw").c_str(), (SampleNames[i]+"_drawfluc_draw").c_str());
767 Get2DBayesianpValue(TempHistogram);
768 TempHistogram->Write();
769 delete TempHistogram;
770
771 // Also write the 2D histograms for the p-value
772 OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_predfluc_draw>>htemp2").c_str());
773 TH2D *TempHistogram2 = static_cast<TH2D*>(gDirectory->Get("htemp2")->Clone());
774 TempHistogram2->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
775 TempHistogram2->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
776 TempHistogram2->SetNameTitle((SampleNames[i]+"_predfluc_draw").c_str(), (SampleNames[i]+"_predfluc_draw").c_str());
777 Get2DBayesianpValue(TempHistogram2);
778 TempHistogram2->Write();
779 delete TempHistogram2;
780
781 // finally p-value for 1D projection
782 OutputTree->Draw((SampleName+"_rate_data_draw:"+SampleName+"_rate_predfluc_draw>>htemp3").c_str());
783 TH2D *TempHistogram3 = static_cast<TH2D*>(gDirectory->Get("htemp3")->Clone());
784 TempHistogram3->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
785 TempHistogram3->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
786 TempHistogram3->SetNameTitle((SampleNames[i]+"_rate_predfluc_draw").c_str(), (SampleNames[i]+"_rate_predfluc_draw").c_str());
787 Get2DBayesianpValue(TempHistogram3);
788 TempHistogram3->Write();
789 delete TempHistogram3;
790
791 // finally p-value for 1D projection
792 OutputTree->Draw((SampleName+"_data_draw_ProjectX:"+SampleName+"_drawfluc_draw_ProjectX>>htemp4").c_str());
793 TH2D *TempHistogram4 = static_cast<TH2D*>(gDirectory->Get("htemp4")->Clone());
794 TempHistogram4->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
795 TempHistogram4->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
796 TempHistogram4->SetNameTitle((SampleNames[i]+"_drawfluc_draw_ProjectX").c_str(), (SampleNames[i]+"_drawfluc_draw_ProjectX").c_str());
797 Get2DBayesianpValue(TempHistogram4);
798 TempHistogram4->Write();
799 delete TempHistogram4;
800
801 // Write the Histograms to each folder
802 DataHist[i]->Write();
803 NominalHist[i]->Write();
804 MeanHist[i]->Write();
805 ModeHist[i]->Write();
806 RatioHistMean->Write();
807 RatioHistMode->Write();
808 RatioHistNom->Write();
809 if(DoBetaParam) MeanHistCorrected[i]->Write();
810
811 W2NomHist[i]->Write();
812 W2MeanHist[i]->Write();
813 W2ModeHist[i]->Write();
814
815 DataNormHist->Write();
816 NomNormHist->Write();
817 MeanNormHist->Write();
818 ModeNormHist->Write();
819
820 DataHist_ProjectX[i]->Write();
821 NomProjectX->Write();
822 MeanProjectX->Write();
823 ModeProjectX->Write();
824 if(DoBetaParam) MeanHistCorrectedProjectX->Write();
825 ViolinHists_ProjectX[i]->Write();
826
827 DataHist_ProjectY[i]->Write();
828 NomProjectY->Write();
829 MeanProjectY->Write();
830 ModeProjectY->Write();
831 if(DoBetaParam) MeanHistCorrectedProjectY->Write();
832 ViolinHists_ProjectY[i]->Write();
833
834 W2NomProjectX->Write();
835 W2MeanProjectX->Write();
836 W2ModeProjectX->Write();
837
838 W2NomProjectY->Write();
839 W2MeanProjectY->Write();
840 W2ModeProjectY->Write();
841
842 //KS: This will dump lots of hists, use it only for debugging
843 if(Debug > 0)
844 {
845 TDirectory* DebugDir = Dir[i]->mkdir("Debug");
846 DebugDir->cd();
847 for (int b = 1; b <= maxBins[i]; ++b)
848 {
849 PosteriorHist[i][b]->Write();
850 std::string Title = PosteriorHist[i][b]->GetName();
851
852 auto TempLine = std::make_unique<TLine>(NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
853 NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
854 TempLine->SetLineColor(kRed);
855 TempLine->SetLineWidth(2);
856
857 auto TempLineData = std::make_unique<TLine>(DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
858 DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
859 TempLineData->SetLineColor(kGreen);
860 TempLineData->SetLineWidth(2);
861
862 // Also fit a Gaussian because why not?
863 TF1 *Fitter = new TF1("Fit", "gaus", PosteriorHist[i][b]->GetBinLowEdge(1), PosteriorHist[i][b]->GetBinLowEdge(PosteriorHist[i][b]->GetNbinsX()+1));
864 PosteriorHist[i][b]->Fit(Fitter, "RQ");
865 Fitter->SetLineColor(kRed-5);
866
867 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
868 Legend->SetFillColor(0);
869 Legend->SetFillStyle(0);
870 Legend->SetLineWidth(0);
871 Legend->SetLineColor(0);
872 Legend->AddEntry(TempLineData.get(), Form("Data #mu=%.2f", DataHist[i]->GetBinContent(b)), "l");
873 Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.2f", NominalHist[i]->GetBinContent(b)), "l");
874 Legend->AddEntry(PosteriorHist[i][b].get(), Form("Post, #mu=%.2f#pm%.2f", PosteriorHist[i][b]->GetMean(), PosteriorHist[i][b]->GetRMS()), "l");
875 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.2f#pm%.2f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
876 std::string TempTitle = std::string(PosteriorHist[i][b]->GetName());
877
878 TempTitle += "_canv";
879 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
880 TempCanvas->SetGridx();
881 TempCanvas->SetGridy();
882 TempCanvas->SetRightMargin(0.03);
883 TempCanvas->SetBottomMargin(0.08);
884 TempCanvas->SetLeftMargin(0.10);
885 TempCanvas->SetTopMargin(0.06);
886 TempCanvas->cd();
887 PosteriorHist[i][b]->Draw();
888 TempLine->Draw("same");
889 TempLineData->Draw("same");
890 Fitter->Draw("same");
891 Legend->Draw("same");
892 TempCanvas->Write();
893
894 delete TempCanvas;
895 delete Fitter;
896 //This isn't useful check only in desperation
897 if(Debug > 1) w2Hist[i][b]->Write();
898 }
899 DebugDir->Close();
900 delete DebugDir;
901 Dir[i]->cd();
902 }
903 lnLHist_Mean[i]->Write();
904 lnLHist_Mode[i]->Write();
905
906 lnLHist_Mean_ProjectX[i]->Write();
907
908 lnLHist_Mean1D[i]->Write();
909 lnLHist_Mode1D[i]->Write();
910
912 lnLHist_Sample_DrawData[i]->Write();
914 lnLHist_Sample_DrawflucDraw[i]->Write();
916 lnLHist_Sample_PredflucDraw[i]->Write();
917
918 if(DoByModePlots)
919 {
920 for (int j = 0; j < Modes->GetNModes()+1; ++j)
921 {
922 MeanHist_ByMode[i][j]->Write();
923 TH1D *MeanProjectX_ByMode = ProjectPoly(MeanHist_ByMode[i][j], true, i, true);
924 TH1D *MeanProjectY_ByMode = ProjectPoly(MeanHist_ByMode[i][j], false, i, true);
925 MeanProjectX_ByMode->Write();
926 MeanProjectY_ByMode->Write();
927 //KS: This will dump lots of hists, use it only for debugging
928 if(Debug > 0)
929 {
930 for (int b = 1; b <= maxBins[i]; ++b)
931 {
932 PosteriorHist_ByMode[i][j][b]->Write();
933 }
934 }
935 delete MeanProjectX_ByMode;
936 delete MeanProjectY_ByMode;
937 } // End loop over bins
938 }
939 // Delete temporary objects
940 delete RatioHistMean;
941 delete RatioHistMode;
942 delete RatioHistNom;
943
944 delete DataNormHist;
945 delete MeanNormHist;
946 delete ModeNormHist;
947 delete NomNormHist;
948
949 delete DataHist_ProjectX[i];
950 delete MeanProjectX;
951 delete ModeProjectX;
952 if(DoBetaParam) delete MeanHistCorrectedProjectX;
953 delete NomProjectX;
954
955 delete DataHist_ProjectY[i];
956 delete MeanProjectY;
957 delete ModeProjectY;
958 if(DoBetaParam) delete MeanHistCorrectedProjectY;
959 delete NomProjectY;
960
961 delete W2NomProjectX;
962 delete W2MeanProjectX;
963 delete W2ModeProjectX;
964
965 delete W2NomProjectY;
966 delete W2MeanProjectY;
967 delete W2ModeProjectY;
968 MACH3LOG_INFO("");
969 } //end loop over samples
971
973 MACH3LOG_INFO("Wrote to {}", Outputfile->GetName());
974}
TH2Poly * RatioPolys(TH2Poly *NumHist, TH2Poly *DenomHist)
Helper to make ratio of TH2Polys.
TH2Poly * NormalisePoly(TH2Poly *Histogram)
WP: Helper to Normalise histograms.
void StudyInformationCriterion(M3::kInfCrit Criterion)
Information Criterion.
void CalcLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly, will modify MC hist tittle to include LLH.
void MakePredictive()
Finalise the distributions from the thrown samples.
void PrepareOutput()
KS: Prepare output tree and necessary variables.
void PlotBetaParameters()
KS: In Barlow Beeston we have Beta Parameters which scale generated MC.
void StudyKinematicCorrelations()
KS: Study how correlated are sample or kinematic bins.

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: