Class responsible for handling implementation of samples used in analysis, reweighting and returning LLH. More...

#include <Samples/SampleHandlerInterface.h>

Inheritance diagram for SampleHandlerInterface:

Collaboration diagram for SampleHandlerInterface:

Public Member Functions
	SampleHandlerInterface ()
	The main constructor. More...

virtual	~SampleHandlerInterface ()
	destructor More...

virtual M3::int_t	GetNSamples ()
	returns total number of samples More...

virtual std::string	GetSampleTitle (const int iSample) const =0
	Get fancy title for specified samples. More...

virtual std::string	GetName () const =0
	Get name for Sample Handler. More...

virtual double	GetSampleLikelihood (const int iSample) const =0
	Get likelihood (-logL) for a single sample. More...

virtual void	CleanMemoryBeforeFit ()=0
	Allow to clean not used memory before fit starts. More...

virtual void	SaveAdditionalInfo (TDirectory *Dir)
	Store additional info in a chain. More...

MaCh3Modes *	GetMaCh3Modes () const
	Return pointer to MaCh3 modes. More...

virtual void	Reweight ()=0
	main routine modifying MC prediction based on proposed parameter values More...

virtual double	GetLikelihood () const =0
	Return likelihood (-logL) for all samples. More...

virtual void	PrintRates (const bool DataOnly=false)=0
	Helper function to print rates for the samples with LLH. More...

unsigned int	GetNEvents () const
	Return total number of events. More...

virtual int	GetNOscChannels (const int iSample) const =0
	Get number of oscillation channels for a single sample. More...

virtual std::string	GetKinVarName (const int iSample, const int Dimension) const =0
	Return Kinematic Variable name for specified sample and dimension for example "Reconstructed_Neutrino_Energy". More...

virtual const TH1 *	GetDataHist (const int Sample)=0
	Get Data histogram. More...

virtual const TH1 *	GetMCHist (const int Sample)=0
	Get MC histogram. More...

virtual const TH1 *	GetW2Hist (const int Sample)=0
	Get W2 histogram. More...

virtual int	GetNDim (const int Sample) const =0
	DB Get what dimensionality binning for given sample has. More...

virtual std::string	GetFlavourName (const int iSample, const int iChannel) const =0
	Get the flavour name for a given sample and oscillation channel. More...

virtual std::vector< double >	ReturnKinematicParameterBinning (const int Sample, const std::string &KinematicParameter) const =0
	Return the binning used to draw a kinematic parameter. More...

virtual std::unique_ptr< TH1 >	Get1DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_Str, const int kModeToFill=-1, const int kChannelToFill=-1, const int WeightStyle=0)=0
	Build a 1D histogram for a given variable, optionally filtered by mode and channel. More...

virtual std::unique_ptr< TH2 >	Get2DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_StrX, const std::string &ProjectionVar_StrY, int kModeToFill=-1, const int kChannelToFill=-1, const int WeightStyle=0)=0
	Build a 2D histogram for given variables, optionally filtered by mode and channel. More...

virtual std::unique_ptr< TH1 >	Get1DVarHist (const int iSample, const std::string &ProjectionVar, const std::vector< KinematicCut > &EventSelectionVec={}, int WeightStyle=0, const std::vector< KinematicCut > &SubEventSelectionVec={})=0
	Return 1D projection of MC into given 1D variable (doesn't have to be variable used in the fit) More...

virtual std::unique_ptr< TH2 >	Get2DVarHist (const int iSample, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &EventSelectionVec={}, const int WeightStyle=0, const std::vector< KinematicCut > &SubEventSelectionVec={})=0
	Build a 2D projection of MC events into specified variables. More...

double	GetPoissonLLH (const double data, const double mc) const
	Calculate test statistic for a single bin using Poisson. More...

double	GetTestStatLLH (const double data, const double mc, const double w2) const
	Calculate test statistic for a single bin. Calculation depends on setting of `fTestStatistic`. Data and mc -> 0 cut-offs are defined in M3::LOW_MC_BOUND. More...

void	SetTestStatistic (TestStatistic testStat)
	Set the test statistic to be used when calculating the binned likelihoods. More...

TestStatistic	GetTestStatistic () const
	Get the test statistic used when calculating the binned likelihoods. More...

Protected Member Functions
void	QuietPlease ()
	CW: Redirect std::cout to silence some experiment specific libraries. More...

void	NowTalk ()
	CW: Redirect std::cout to silence some experiment specific libraries. More...

template<typename T >
bool	MatchCondition (const std::vector< T > &allowedValues, const T &value)
	check if event is affected by following conditions, for example pdg, or modes etc More...

Protected Attributes
TestStatistic	fTestStatistic
	Test statistic tells what kind of likelihood sample is using. More...

std::streambuf *	buf
	Keep the cout buffer. More...

std::streambuf *	errbuf
	Keep the cerr buffer. More...

M3::int_t	nSamples
	Contains how many samples we've got. More...

unsigned int	nEvents
	Number of MC events are there. More...

std::unique_ptr< MaCh3Modes >	Modes
	Holds information about used Generator and MaCh3 modes. More...

Detailed Description

Class responsible for handling implementation of samples used in analysis, reweighting and returning LLH.

It serves as basic interface for fit running, as well as other fucnitonalsities liek llh scan, sigma var or event posterior predictive disturibon. Concrete implementations of this interface are responsible for defining the specific structure of samples, event selections, and histogram filling.

Author: Asher Kaboth; Richard Calland

Definition at line 34 of file SampleHandlerInterface.h.

Constructor & Destructor Documentation

◆ SampleHandlerInterface()

SampleHandlerInterface::SampleHandlerInterface ( )

The main constructor.

Definition at line 4 of file SampleHandlerInterface.cpp.

                                                {
 // ***************************************************************************
   nEvents = 0;
   nSamples = 0;
 }

◆ ~SampleHandlerInterface()

SampleHandlerInterface::~SampleHandlerInterface ( )

virtual

destructor

Definition at line 11 of file SampleHandlerInterface.cpp.

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

Member Function Documentation

◆ CleanMemoryBeforeFit()

virtual void SampleHandlerInterface::CleanMemoryBeforeFit ( )

pure virtual

Allow to clean not used memory before fit starts.

Implemented in PySampleHandlerBase, PySampleHandlerInterface, and SampleHandlerNuDockBase.

◆ Get1DVarHist()

virtual std::unique_ptr<TH1> SampleHandlerInterface::Get1DVarHist	(	const int	iSample,
		const std::string &	ProjectionVar,
		const std::vector< KinematicCut > &	EventSelectionVec = `{}`,
		int	WeightStyle = `0`,
		const std::vector< KinematicCut > &	SubEventSelectionVec = `{}`
	)

pure virtual

Return 1D projection of MC into given 1D variable (doesn't have to be variable used in the fit)

Parameters

iSample	Index of the sample.
ProjectionVar	name of variable
EventSelectionVec	Vector of additional cuts like cut on interaction mode
WeightStyle	Alow to modify weight for example if equal to 1 all weights are set to 1
SubEventSelectionVec	Vector of additional cuts for sub event (particle, ring etc.)

Implemented in SampleHandlerNuDockBase, SampleHandlerBase, and PySampleHandlerInterface.

◆ Get1DVarHistByModeAndChannel()

virtual std::unique_ptr<TH1> SampleHandlerInterface::Get1DVarHistByModeAndChannel	(	const int	iSample,
		const std::string &	ProjectionVar_Str,
		const int	kModeToFill = `-1`,
		const int	kChannelToFill = `-1`,
		const int	WeightStyle = `0`
	)

pure virtual

Build a 1D histogram for a given variable, optionally filtered by mode and channel.

Parameters

iSample	Index of the sample.
ProjectionVar_Str	Name of the variable to project onto.
kModeToFill	Interaction mode to select (-1 means all modes).
kChannelToFill	Oscillation channel to select (-1 means all channels).
WeightStyle	Weighting scheme (e.g. 0 = nominal weights, 1 = unit weights).

Implemented in SampleHandlerNuDockBase, PySampleHandlerInterface, and SampleHandlerBase.

◆ Get2DVarHist()

virtual std::unique_ptr<TH2> SampleHandlerInterface::Get2DVarHist	(	const int	iSample,
		const std::string &	ProjectionVarX,
		const std::string &	ProjectionVarY,
		const std::vector< KinematicCut > &	EventSelectionVec = `{}`,
		const int	WeightStyle = `0`,
		const std::vector< KinematicCut > &	SubEventSelectionVec = `{}`
	)

pure virtual

Build a 2D projection of MC events into specified variables.

Parameters

iSample	Index of the sample.
ProjectionVarX	Name of the variable for the X axis.
ProjectionVarY	Name of the variable for the Y axis.
EventSelectionVec	Vector of event-level selection cuts.
WeightStyle	Weighting scheme (e.g. 0 = nominal weights, 1 = unit weights).
SubEventSelectionVec	Vector of sub-event selection cuts.

Implemented in SampleHandlerNuDockBase, SampleHandlerBase, and PySampleHandlerInterface.

◆ Get2DVarHistByModeAndChannel()

virtual std::unique_ptr<TH2> SampleHandlerInterface::Get2DVarHistByModeAndChannel	(	const int	iSample,
		const std::string &	ProjectionVar_StrX,
		const std::string &	ProjectionVar_StrY,
		int	kModeToFill = `-1`,
		const int	kChannelToFill = `-1`,
		const int	WeightStyle = `0`
	)

pure virtual

Build a 2D histogram for given variables, optionally filtered by mode and channel.

Parameters

iSample	Index of the sample.
ProjectionVar_StrX	Name of the variable for the X axis.
ProjectionVar_StrY	Name of the variable for the Y axis.
kModeToFill	Interaction mode to select (-1 means all modes).
kChannelToFill	Oscillation channel to select (-1 means all channels).
WeightStyle	Weighting scheme (e.g. 0 = nominal weights, 1 = unit weights).

Implemented in SampleHandlerNuDockBase, PySampleHandlerInterface, and SampleHandlerBase.

◆ GetDataHist()

virtual const TH1* SampleHandlerInterface::GetDataHist ( const int Sample )

pure virtual

Get Data histogram.

Parameters

Sample Sample enumerator

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetFlavourName()

virtual std::string SampleHandlerInterface::GetFlavourName	(	const int	iSample,
		const int	iChannel
	)		const

pure virtual

Get the flavour name for a given sample and oscillation channel.

Parameters

iSample	Index of the sample.
iChannel	Index of the oscillation channel within the sample.

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetKinVarName()

virtual std::string SampleHandlerInterface::GetKinVarName	(	const int	iSample,
		const int	Dimension
	)		const

pure virtual

Return Kinematic Variable name for specified sample and dimension for example "Reconstructed_Neutrino_Energy".

Parameters

iSample	Sample index
Dimension	Dimension index

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetLikelihood()

virtual double SampleHandlerInterface::GetLikelihood ( ) const

pure virtual

Return likelihood (-logL) for all samples.

Implemented in SampleHandlerBase, PySampleHandlerInterface, and SampleHandlerNuDockBase.

◆ GetMaCh3Modes()

MaCh3Modes* SampleHandlerInterface::GetMaCh3Modes ( ) const

inline

Return pointer to MaCh3 modes.

Definition at line 58 of file SampleHandlerInterface.h.

58 { return Modes.get(); }

SampleHandlerInterface::Modes

std::unique_ptr< MaCh3Modes > Modes

Holds information about used Generator and MaCh3 modes.

Definition: SampleHandlerInterface.h:269

◆ GetMCHist()

virtual const TH1* SampleHandlerInterface::GetMCHist ( const int Sample )

pure virtual

Get MC histogram.

Parameters

Sample Sample enumerator

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetName()

virtual std::string SampleHandlerInterface::GetName ( ) const

pure virtual

Get name for Sample Handler.

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetNDim()

virtual int SampleHandlerInterface::GetNDim ( const int Sample ) const

pure virtual

DB Get what dimensionality binning for given sample has.

Parameters

Sample Number of sample

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetNEvents()

unsigned int SampleHandlerInterface::GetNEvents ( ) const

inline

Return total number of events.

Definition at line 69 of file SampleHandlerInterface.h.

69 {return nEvents;}

◆ GetNOscChannels()

virtual int SampleHandlerInterface::GetNOscChannels ( const int iSample ) const

pure virtual

Get number of oscillation channels for a single sample.

Parameters

iSample Sample enumerator

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ GetNSamples()

virtual M3::int_t SampleHandlerInterface::GetNSamples ( )

inlinevirtual

returns total number of samples

Definition at line 43 of file SampleHandlerInterface.h.

43 { return nSamples; };

◆ GetPoissonLLH()

double SampleHandlerInterface::GetPoissonLLH	(	const double	data,
		const double	mc
	)		const

Calculate test statistic for a single bin using Poisson.

Parameters

data	is data
mc	is mc

Definition at line 17 of file SampleHandlerInterface.cpp.

                                                                                      {
 // ***************************************************************************
   // Return MC if there are no data, returns 0 for data == 0 && mc == 0  
   if ( data == 0 ) return mc;
  
   // If there are some data, but the prediction falls below the MC bound => return Poisson LogL for the low MC bound
   if ( mc < M3::_LOW_MC_BOUND_ ) {
     if ( data > M3::_LOW_MC_BOUND_ ) return ( M3::_LOW_MC_BOUND_ - data + data * std::log( data/M3::_LOW_MC_BOUND_ ) );
     else if ( data >= mc ) return 0.;
   }
   
   // Otherwise, just return usual Poisson LogL using Stirling's approximation
   // http://hyperphysics.phy-astr.gsu.edu/hbase/math/stirling.html
   return ( mc - data + data * std::log( data / mc ) );
 }

◆ GetSampleLikelihood()

virtual double SampleHandlerInterface::GetSampleLikelihood ( const int iSample ) const

pure virtual

Get likelihood (-logL) for a single sample.

Parameters

iSample Sample enumerator

Implemented in SampleHandlerBase, PySampleHandlerInterface, and SampleHandlerNuDockBase.

◆ GetSampleTitle()

virtual std::string SampleHandlerInterface::GetSampleTitle ( const int iSample ) const

pure virtual

Get fancy title for specified samples.

Parameters

iSample Sample enumerator

Implemented in SampleHandlerNuDockBase, SampleHandlerBase, and PySampleHandlerInterface.

◆ GetTestStatistic()

TestStatistic SampleHandlerInterface::GetTestStatistic ( ) const

inline

Get the test statistic used when calculating the binned likelihoods.

Definition at line 237 of file SampleHandlerInterface.h.

237 { return fTestStatistic; }

SampleHandlerInterface::fTestStatistic

TestStatistic fTestStatistic

Test statistic tells what kind of likelihood sample is using.

Definition: SampleHandlerInterface.h:255

◆ GetTestStatLLH()

double SampleHandlerInterface::GetTestStatLLH	(	const double	data,
		const double	mc,
		const double	w2
	)		const

Calculate test statistic for a single bin. Calculation depends on setting of fTestStatistic. Data and mc -> 0 cut-offs are defined in M3::LOW_MC_BOUND.

Poisson

Standard Poisson log-likelihood (Stirling approximation) [2]

\[ - \log \mathcal{L}_\mathrm{Poisson} = \sum_i N_i^\mathrm{MC} - N_i^\mathrm{data} + N_i^\mathrm{data} \ln \frac{N_i^\mathrm{data}}{N_i^\mathrm{MC}}, \]

Pearson

Standard Pearson likelihood [27] (assumes Gaussian approximation of bin counts):

\[ - \log \mathcal{L}_\mathrm{Pearson} = \sum_i \frac{(N_i^\mathrm{data} - N_i^\mathrm{MC})^2}{2 \, N_i^\mathrm{MC}} \]

Barlow-Beeston

Based on [3] and following Conway approximation ([5]) The generation of MC is a stochastic process, so even identical settings can lead to different outputs (assuming that the seeds of the random number generator are different). This introduces uncertainty in MC distributions, especially in bins with low statistics.

\[ - \log \mathcal{L}_\mathrm{BB} = - \log \mathcal{L}_\mathrm{Poisson} - \log \mathcal{L}_\mathrm{MC_{stat}} = \sum_i \Biggl[ N_i^\mathrm{MC}(\vec{\theta}) - N_i^\mathrm{data} + N_i^\mathrm{data} \ln \frac{N_i^\mathrm{data}}{N_i^\mathrm{MC}(\vec{\theta})} + \frac{(\beta_i - 1)^2}{2 \sigma_{\beta_i}^2} \Biggr], \]

where \(\beta_i\) is a scaling parameter between ideal ("true") and generated MC in a bin ( \(N^\mathrm{true}_{\mathrm{MC},i} = \beta_i N_i^\mathrm{MC}\)), and \(\sigma^2_{\beta_i} = \frac{\sum_i w_i^2}{N_i^\mathrm{MC}}\), with \(\sum_i w_i^2\) being the sum of the squares of weights in bin \(i\). Assuming \(\beta_i\) follows a Gaussian, its mean can be found by solving the quadratic equation derived by Conway:

\[ \beta_i^2 + (N_i^\mathrm{MC} \sigma_{\beta_i}^2 - 1)\beta_i - N_i^\mathrm{data} \sigma_{\beta_i}^2 = 0 \]

Dembinski-Abdelmotteleb

Alternative treatment of MC statistical uncertainty following Hans Dembinski and Ahmed Abdelmotteleb [6]

This approach extends the Barlow-Beeston method. For each bin:

\[ - \log \mathcal{L}_\mathrm{DA} = (N_i^{\mathrm{MC},\prime} - N_i^\mathrm{data} + N_i^\mathrm{data} \ln \frac{N_i^\mathrm{data}}{N_i^{\mathrm{MC},\prime}}) + k \beta - k + k \ln \frac{k}{k \beta} \]

where

\[ k = \frac{(N_i^\mathrm{MC})^2}{\sum_i w_i^2} \]

and

\[ \beta = \frac{N_i^\mathrm{data} + k}{N_i^\mathrm{MC} + k}, \quad N_i^{\mathrm{MC},\prime} = N_i^\mathrm{MC} \cdot \beta \]

IceCube

Alternative likelihood definition described by the IceCube collaboration [1]

\[ - \log \mathcal{L} = - \sum_i \Biggl( a_i \log(b_i) + \log[\Gamma(N_i^{\mathrm{data}}+a_i)] - (N_i^{\mathrm{data}}+a_i)\log(b_i+1) - \log[\Gamma(a_i)] \Biggr), \]

where the auxiliary variables are

\[ a_i = N^{\mathrm{gen}}_{\mathrm{MC},i} \, b_i + 1, \quad b_i = \frac{N^{\mathrm{gen}}_{\mathrm{MC},i}}{\sum_i w_i^2}. \]

Treatment of low data/mc

Implemented fTestStatistic are kPoisson (with Stirling's approx.), kBarlowBeeston (arXiv:1103.0354), kDembinskiAbdelmotteleb (arXiv:2206.12346), kIceCube (arxiv:1901.04645), and kPearson. Test statistics require mc > 0, therefore low mc and data values are treated with cut-offs based on M3::LOW_MC_BOUND = .00001 by default. For kPoisson, kBarlowBeeston, kDembinskiAbdelmotteleb, kPearson: data > LOW_MC_BOUND & mc <= LOW_MC_BOUND: returns GetTestStatLLH(data, LOW_MC_BOUND, w2), with Poisson(data,LOW_MC_BOUND) limit for mc->0, w2->0. mc < data <= LOW_MC_BOUND: returns 0 (as if any data <= LOW_MC_BOUND were effectively consistent with 0 data count), with a limit of 0 for mc->0. data = 0: returns mc (or mc/2. for kPearson), with a limit of 0 for mc->0. For kIceCube: mc < data returns the lower of IceCube(data,mc,w2) and Poisson(data,mc) penalties, with a Poisson(data,LOW_MC_BOUND) limit for mc->0, w2->0.

Parameters

data	is data
mc	is mc
w2	is \(\sum_{i} w_{i}^2\) (sum of weights squared), which is \(\sigma^2_{\text{MC stats}}\)

Definition at line 35 of file SampleHandlerInterface.cpp.

                                                                                                        {
 // *************************
   switch (fTestStatistic)
   {
     //CW: Not full Barlow-Beeston or what is referred to as "light": we're not introducing any more parameters
     // Assume the MC has a Gaussian distribution around generated
     // As in https://arxiv.org/abs/1103.0354 eq 10, 11
     //CW: Calculate the Barlow-Beeston likelihood contribution from MC statistics
     // Assumes the beta scaling parameters are Gaussian distributed
     // Follows arXiv:1103.0354 section 5 and equation 8, 9, 10, 11 on page 4/5
     // Essentially solves equation 11
     case (kBarlowBeeston):
     {
       // The MC used in the likelihood calculation is allowed to be changed by Barlow Beeston beta parameters
       double newmc = mc;
  
       // If MC falls below the low MC bound, use low MC bound for newmc
       if ( mc < M3::_LOW_MC_BOUND_ ) {
         if ( data > M3::_LOW_MC_BOUND_ ) newmc = M3::_LOW_MC_BOUND_;
         else if ( data >= mc ) return 0.;
       }
         
       // Barlow-Beeston uses fractional uncertainty on MC, so sqrt(sum[w^2])/mc
       const double fractional = std::sqrt( w2 ) / newmc;
       // fractional^2 to avoid doing same operation several times
       const double fractional2 = fractional * fractional;
       // b in quadratic equation
       const double temp = newmc * fractional2 - 1;
       // b^2 - 4ac in quadratic equation
       const double temp2 = temp * temp + 4 * data * fractional2;
       if ( temp2 < 0 ) {
         MACH3LOG_ERROR("Negative square root in Barlow Beeston coefficient calculation!");
         throw MaCh3Exception(__FILE__ , __LINE__ );
       }
       // Solve for the positive beta
       const double beta = ( -1 * temp + sqrt( temp2 ) ) / 2.;
  
       // If there is no data, test-stat shall return only MC*beta
       double stat = mc * beta;
       // With data, test-stat shall return LogL for newMC*beta which includes the low MC bound
       if ( data > 0 ) {
         newmc *= beta;
         stat = newmc - data + data * std::log( data / newmc );
       }
  
       // Now, MC stat penalty
       // The penalty from MC statistics using Conways approach (https://cds.cern.ch/record/1333496?)
       double penalty = 0;
       if ( fractional > 0 )  penalty = ( beta - 1 ) * ( beta - 1 ) / ( 2 * fractional2 );
  
       // Returns test-stat plus the MC stat penalty 
       return stat+penalty;
     }
     break;
     //KS: Alternative calculation of Barlow-Beeston following Hans Dembinski and Ahmed Abdelmottele arXiv:2206.12346v2
     case (kDembinskiAbdelmotteleb):
     {
       //KS: code follows authors implementation from:
       //https://github.com/scikit-hep/iminuit/blob/059d06b00cae097ebf340b218b4eb57357111df8/src/iminuit/cost.py#L274-L300
       
       // If there is no MC stat error for any reason, return Poisson LogL
       if ( w2 == 0 ) return GetPoissonLLH(data,mc);
       
       // The MC can be changed
       double newmc = mc;
       
       // If MC falls below the low MC bound, use low MC bound for newmc
       if ( mc < M3::_LOW_MC_BOUND_ ) {
         if ( data > M3::_LOW_MC_BOUND_ ) newmc = M3::_LOW_MC_BOUND_;
         else if ( data >= mc ) return 0.;
       }
       
       //the so-called effective count
       const double k = newmc * newmc / w2;
       //Calculate beta which is scaling factor between true and generated MC
       const double beta = ( data + k ) / ( newmc + k );
  
       newmc *= beta;
       
       // And penalise the movement in beta relative the mc uncertainty
       const double penalty = k * beta - k + k * std::log( k / ( k * beta ) );
  
       // If there are no data, this shall return newmc
       double stat = newmc;
       // All likelihood calculations may use the bare Poisson likelihood, so calculate here
       // Only if there are some data
       if ( data > 0 ) stat = newmc - data + data * std::log( data / newmc );
  
       // Return the statistical contribution and penalty
       return stat+penalty;
     }
     break;
     //CW: Also try the IceCube likelihood
     // It does not modify the MC content
     // https://arxiv.org/abs/1901.04645
     // Argüelles, C.A., Schneider, A. & Yuan, T. J. High Energ. Phys. (2019) 2019: 30. https://doi.org/10.1007/JHEP06(2019)030
     // We essentially construct eq 3.16 and take the logarithm
     // in eq 3.16, mu is MC, sigma2 is w2, k is data
     case (kIceCube):
     {
       // IceCube low MC bound is implemented to return Poisson(data, _LOW_MC_BOUND_)
       // up until the IceCube(data, mc) test-statistic is less than Poisson(data, _LOW_MC_BOUND_)
       // The 0 MC limit is set to Poisson(data, 0.) as there is no way to get a non-diverging and reasonable guess on w2
  
       // If there is 0 MC uncertainty (technically also when MC is 0) => Return Poisson(data, mc)
       if ( w2 == 0 ) return GetPoissonLLH(data,mc);
       
       // Auxiliary variables
       const long double b = mc / w2;
       const long double a = mc * b + 1;
       
       // Use C99's implementation of log of gamma function to not be C++11 dependent
       const double stat = double( -1 * ( a * logl( b ) + lgammal( data + a ) - lgammal( data + 1 ) - ( ( data + a ) * log1pl( b ) ) - lgammal( a ) ) );
  
       // Check whether the stat is more than Poisson-like bound for low mc (mc < data)
       // TN: I believe this might get some extra optimization
       if ( mc <= data ) {
         if ( data <= M3::_LOW_MC_BOUND_ ) return 0.;
         const double poisson = GetPoissonLLH(data, M3::_LOW_MC_BOUND_);
         if ( stat > poisson ) return poisson;
       }
  
       // Otherwise, return IceCube test-stat
       return stat;
     }
     break;
     //KS: Pearson works on assumption that event distribution in each bin is described by a Gaussian which in our case is not fulfilled for all bins, hence use it at your own risk
     case (kPearson):
     {
       //KS: 2 is because this function returns -LLH not -2LLH
       // With no data return the MC/2.
       if ( data == 0 ) return mc/2.;
       
       // If MC is lower than the low MC bound, return the test-stat at the bound
       if ( mc < M3::_LOW_MC_BOUND_ ) {
         if ( data > M3::_LOW_MC_BOUND_ ) return ( data - M3::_LOW_MC_BOUND_ ) * ( data - M3::_LOW_MC_BOUND_ ) / ( 2. * M3::_LOW_MC_BOUND_ ); 
         else if ( data >= mc ) return 0.;
       }
       
       // Return the Pearson metric
       return ( data - mc ) * ( data - mc ) / ( 2 * mc );
     }
     break;
     case (kPoisson):
     {
       //Just call GetPoissonLLH which doesn't take in weights
       //and is a Poisson likelihood comparison.
       return GetPoissonLLH(data, mc);
     }
     break;
     case TestStatistic::kNTestStatistics:
       MACH3LOG_ERROR("kNTestStatistics is not a valid TestStatistic!");
       throw MaCh3Exception(__FILE__, __LINE__);
     default:
     MACH3LOG_ERROR("Couldn't find TestStatistic {} exiting!", static_cast<int>(fTestStatistic));
     throw MaCh3Exception(__FILE__ , __LINE__ );
   } // end switch
 }

◆ GetW2Hist()

virtual const TH1* SampleHandlerInterface::GetW2Hist ( const int Sample )

pure virtual

Get W2 histogram.

Parameters

Sample Sample enumerator

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ MatchCondition()

template<typename T >

bool SampleHandlerInterface::MatchCondition	(	const std::vector< T > &	allowedValues,
		const T &	value
	)

inlineprotected

check if event is affected by following conditions, for example pdg, or modes etc

Definition at line 247 of file SampleHandlerInterface.h.

                                                                          {
     if (allowedValues.empty()) {
       return true;  // Apply to all if no specific values are specified
     }
     return std::find(allowedValues.begin(), allowedValues.end(), value) != allowedValues.end();
   }

◆ NowTalk()

void SampleHandlerInterface::NowTalk ( )

protected

CW: Redirect std::cout to silence some experiment specific libraries.

Definition at line 210 of file SampleHandlerInterface.cpp.

                                      {
 // ***************************************************************************
   #if MACH3_DEBUG > 0
   return;
   #else
   std::cout.rdbuf(buf);
   std::cerr.rdbuf(errbuf);
   #endif
 }

◆ PrintRates()

virtual void SampleHandlerInterface::PrintRates ( const bool DataOnly = false )

pure virtual

Helper function to print rates for the samples with LLH.

Parameters

DataOnly whether to print data only rates

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ QuietPlease()

void SampleHandlerInterface::QuietPlease ( )

protected

CW: Redirect std::cout to silence some experiment specific libraries.

Definition at line 196 of file SampleHandlerInterface.cpp.

                                          {
 // ***************************************************************************
   #if MACH3_DEBUG > 0
   return;
   #else
   buf = std::cout.rdbuf();
   errbuf = std::cerr.rdbuf();
   std::cout.rdbuf( nullptr );
   std::cerr.rdbuf( nullptr );
   #endif
 }

◆ ReturnKinematicParameterBinning()

virtual std::vector<double> SampleHandlerInterface::ReturnKinematicParameterBinning	(	const int	Sample,
		const std::string &	KinematicParameter
	)		const

pure virtual

Return the binning used to draw a kinematic parameter.

Parameters

iSample	Index of the sample.
KinematicParameter	name of variable

Implemented in PySampleHandlerInterface, SampleHandlerNuDockBase, and SampleHandlerBase.

◆ Reweight()

virtual void SampleHandlerInterface::Reweight ( )

pure virtual

main routine modifying MC prediction based on proposed parameter values

Implemented in SampleHandlerBase, PySampleHandlerInterface, and SampleHandlerNuDockBase.

◆ SaveAdditionalInfo()

virtual void SampleHandlerInterface::SaveAdditionalInfo ( TDirectory * Dir )

inlinevirtual

Store additional info in a chain.

Parameters

Dir	directory to which we save additional info

Reimplemented in SampleHandlerBase.

Definition at line 56 of file SampleHandlerInterface.h.

56 {(void) Dir;};

◆ SetTestStatistic()

void SampleHandlerInterface::SetTestStatistic ( TestStatistic testStat )

inline

Set the test statistic to be used when calculating the binned likelihoods.

Parameters

testStat The test statistic to use.

Definition at line 235 of file SampleHandlerInterface.h.

235 { fTestStatistic = testStat; }

Member Data Documentation

◆ buf

std::streambuf* SampleHandlerInterface::buf

protected

Keep the cout buffer.

Definition at line 258 of file SampleHandlerInterface.h.

◆ errbuf

std::streambuf* SampleHandlerInterface::errbuf

protected

Keep the cerr buffer.

Definition at line 260 of file SampleHandlerInterface.h.

◆ fTestStatistic

TestStatistic SampleHandlerInterface::fTestStatistic

protected

Test statistic tells what kind of likelihood sample is using.

Definition at line 255 of file SampleHandlerInterface.h.

◆ Modes

std::unique_ptr<MaCh3Modes> SampleHandlerInterface::Modes

protected

Holds information about used Generator and MaCh3 modes.

Definition at line 269 of file SampleHandlerInterface.h.

◆ nEvents

unsigned int SampleHandlerInterface::nEvents

protected

Number of MC events are there.

Definition at line 266 of file SampleHandlerInterface.h.

◆ nSamples

M3::int_t SampleHandlerInterface::nSamples

protected

Contains how many samples we've got.

Definition at line 263 of file SampleHandlerInterface.h.

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

Samples/SampleHandlerInterface.h
Samples/SampleHandlerInterface.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ SampleHandlerInterface()

◆ ~SampleHandlerInterface()

Member Function Documentation

◆ CleanMemoryBeforeFit()

◆ Get1DVarHist()

◆ Get1DVarHistByModeAndChannel()

◆ Get2DVarHist()

◆ Get2DVarHistByModeAndChannel()

◆ GetDataHist()

◆ GetFlavourName()

◆ GetKinVarName()

◆ GetLikelihood()

◆ GetMaCh3Modes()

◆ GetMCHist()

◆ GetName()

◆ GetNDim()

◆ GetNEvents()

◆ GetNOscChannels()

◆ GetNSamples()

◆ GetPoissonLLH()

◆ GetSampleLikelihood()

◆ GetSampleTitle()

◆ GetTestStatistic()

◆ GetTestStatLLH()

Poisson

Pearson

Barlow-Beeston

Dembinski-Abdelmotteleb

IceCube

Treatment of low data/mc

◆ GetW2Hist()

◆ MatchCondition()

◆ NowTalk()

◆ PrintRates()

◆ QuietPlease()

◆ ReturnKinematicParameterBinning()

◆ Reweight()

◆ SaveAdditionalInfo()

◆ SetTestStatistic()

Member Data Documentation

◆ buf

◆ errbuf

◆ fTestStatistic

◆ Modes

◆ nEvents

◆ nSamples