Sample Handler Getters

Functions
virtual M3::int_t	SampleHandlerBase::GetNsamples ()
virtual std::string	SampleHandlerBase::GetTitle () const
virtual std::string	SampleHandlerBase::GetSampleName (int Sample) const =0
virtual double	SampleHandlerBase::GetSampleLikelihood (const int isample)
MaCh3Modes *	SampleHandlerBase::GetMaCh3Modes () const
	Return pointer to MaCh3 modes. More...
virtual double	SampleHandlerBase::GetLikelihood ()=0
unsigned int	SampleHandlerBase::GetNEvents () const
virtual int	SampleHandlerBase::GetNOscChannels ()
double	SampleHandlerBase::GetPoissonLLH (const double data, const double mc) const
	Calculate test statistic for a single bin using Poisson. More...
double	SampleHandlerBase::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	SampleHandlerBase::SetTestStatistic (TestStatistic testStat)
	Set the test statistic to be used when calculating the binned likelihoods. More...
int	SampleHandlerFD::GetNDim () const
	DB Function to differentiate 1D or 2D binning. More...
std::string	SampleHandlerFD::GetSampleName (int iSample=0) const override
std::string	SampleHandlerFD::GetTitle () const override
std::string	SampleHandlerFD::GetXBinVarName () const
std::string	SampleHandlerFD::GetYBinVarName () const
TH1 *	SampleHandlerFD::GetMCHist (const int Dimension)
	Get MC histogram. More...
TH1 *	SampleHandlerFD::GetW2Hist (const int Dimension)
	Get W2 histogram. More...
TH1 *	SampleHandlerFD::GetDataHist (const int Dimension)
	Get Data histogram. More...
int	SampleHandlerFD::GetNOscChannels () override
std::string	SampleHandlerFD::GetFlavourName (const int iChannel)
TH1 *	SampleHandlerFD::Get1DVarHist (const std::string &ProjectionVar, const std::vector< KinematicCut > &EventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0, TAxis *Axis=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >())
void	SampleHandlerFD::Fill1DSubEventHist (TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
TH1 *	SampleHandlerFD::Get1DVarHistByModeAndChannel (const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr)
TH2 *	SampleHandlerFD::Get2DVarHistByModeAndChannel (const std::string &ProjectionVar_StrX, const std::string &ProjectionVar_StrY, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr)
TH1 *	SampleHandlerFD::GetModeHist1D (int s, int m, int style=0)
TH2 *	SampleHandlerFD::GetModeHist2D (int s, int m, int style=0)
std::vector< TH1 * >	SampleHandlerFD::ReturnHistsBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
std::vector< TH2 * >	SampleHandlerFD::ReturnHistsBySelection2D (std::string KinematicProjectionX, std::string KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *XAxis=0, TAxis *YAxis=0)
THStack *	SampleHandlerFD::ReturnStackedHistBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
TLegend *	SampleHandlerFD::ReturnStackHistLegend ()
int	SampleHandlerFD::ReturnKinematicParameterFromString (const std::string &KinematicStr) const
	ETA function to generically convert a string from xsec cov to a kinematic type. More...
std::string	SampleHandlerFD::ReturnStringFromKinematicParameter (const int KinematicVariable) const
	ETA function to generically convert a kinematic type from xsec cov to a string. More...

Detailed Description

Group of functions to get various parameters, names, and values.

Function Documentation

Fill1DSubEventHist()

void SampleHandlerFD::Fill1DSubEventHist	(	TH1D *	_h1DVar,
		const std::string &	ProjectionVar,
		const std::vector< KinematicCut > &	SubEventSelectionVec = std::vector< KinematicCut >(),
		int	WeightStyle = 0
	)

Definition at line 1491 of file SampleHandlerFD.cpp.

                                                                                                                                                                  {
  int ProjectionVar_Int = ReturnKinematicVectorFromString(ProjectionVar_Str);
 
  //JM Loop over all events
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    if (IsEventSelected(iEvent)) {
      double Weight = GetEventWeight(iEvent);
      if (WeightStyle == 1) {
        Weight = 1.;
      }
      std::vector<double> Vec = ReturnKinematicVector(ProjectionVar_Int,iEvent);
      size_t nsubevents = Vec.size();
      //JM Loop over all subevents in event
      for (unsigned int iSubEvent = 0; iSubEvent < nsubevents; iSubEvent++) {
        if (IsSubEventSelected(SubEventSelectionVec, iEvent, iSubEvent, nsubevents)) {
          double Var = Vec[iSubEvent];
          _h1DVar->Fill(Var,Weight);
        }
      }
    }
  }
}

Get1DVarHist()

TH1 * SampleHandlerFD::Get1DVarHist	(	const std::string &	ProjectionVar,
		const std::vector< KinematicCut > &	EventSelectionVec = std::vector< KinematicCut >(),
		int	WeightStyle = 0,
		TAxis *	Axis = nullptr,
		const std::vector< KinematicCut > &	SubEventSelectionVec = std::vector< KinematicCut >()
	)

Definition at line 1437 of file SampleHandlerFD.cpp.

                                                                                         {
  //DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
  //   Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
  std::vector< KinematicCut > tmp_Selection = Selection;
  std::vector< KinematicCut > SelectionVecToApply;
 
  //DB Add all the predefined selections to the selection vector which will be applied
  for (size_t iSelec=0;iSelec<Selection.size();iSelec++) {
    SelectionVecToApply.emplace_back(Selection[iSelec]);
  }
 
  //DB Add all requested cuts from the argument to the selection vector which will be applied
  for (size_t iSelec=0;iSelec<EventSelectionVec.size();iSelec++) {
    SelectionVecToApply.emplace_back(EventSelectionVec[iSelec]);
  }
 
  //DB Set the member variable to be the cuts to apply
  Selection = SelectionVecToApply;
 
  //DB Define the histogram which will be returned
  TH1D* _h1DVar;
  if (Axis) {
    _h1DVar = new TH1D("","",Axis->GetNbins(),Axis->GetXbins()->GetArray());
  } else {
    std::vector<double> xBinEdges = ReturnKinematicParameterBinning(ProjectionVar_Str);
    _h1DVar = new TH1D("", "", int(xBinEdges.size())-1, xBinEdges.data());
  }
  _h1DVar->GetXaxis()->SetTitle(ProjectionVar_Str.c_str());
 
  if (IsSubEventVarString(ProjectionVar_Str)) {
    Fill1DSubEventHist(_h1DVar, ProjectionVar_Str, SubEventSelectionVec, WeightStyle);
  } else {
    //DB Grab the associated enum with the argument string
    int ProjectionVar_Int = ReturnKinematicParameterFromString(ProjectionVar_Str);
 
    //DB Loop over all events
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
      if (IsEventSelected(iEvent)) {
        double Weight = GetEventWeight(iEvent);
        if (WeightStyle == 1) {
          Weight = 1.;
        }
        double Var = ReturnKinematicParameter(ProjectionVar_Int,iEvent);
        _h1DVar->Fill(Var,Weight);
      }
    }
  }
  //DB Reset the saved selection
  Selection = tmp_Selection;
 
  return _h1DVar;
}

Get1DVarHistByModeAndChannel()

TH1 * SampleHandlerFD::Get1DVarHistByModeAndChannel	(	const std::string &	ProjectionVar_Str,
		int	kModeToFill = -1,
		int	kChannelToFill = -1,
		int	WeightStyle = 0,
		TAxis *	Axis = nullptr
	)

Definition at line 1728 of file SampleHandlerFD.cpp.

                                                                       {
  bool fChannel;
  bool fMode;
 
  if (kChannelToFill != -1) {
    if (kChannelToFill > GetNOscChannels()) {
      MACH3LOG_ERROR("Required channel is not available. kChannelToFill should be between 0 and {}", GetNOscChannels());
      MACH3LOG_ERROR("kChannelToFill given:{}", kChannelToFill);
      MACH3LOG_ERROR("Exiting.");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    fChannel = true;
  } else {
    fChannel = false;
  }
 
  if (kModeToFill!=-1) {
    if (!(kModeToFill >= 0) && (kModeToFill < Modes->GetNModes())) {
      MACH3LOG_ERROR("Required mode is not available. kModeToFill should be between 0 and {}", Modes->GetNModes());
      MACH3LOG_ERROR("kModeToFill given:{}", kModeToFill);
      MACH3LOG_ERROR("Exiting..");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    fMode = true;
  } else {
    fMode = false;
  }
 
  std::vector< KinematicCut > SelectionVec;
 
  if (fMode) {
    KinematicCut SelecMode;
    SelecMode.ParamToCutOnIt = ReturnKinematicParameterFromString("Mode");
    SelecMode.LowerBound = kModeToFill;
    SelecMode.UpperBound = kModeToFill+1;
    SelectionVec.push_back(SelecMode);
  }
 
  if (fChannel) {
    KinematicCut SelecChannel;
    SelecChannel.ParamToCutOnIt = ReturnKinematicParameterFromString("OscillationChannel");
    SelecChannel.LowerBound = kChannelToFill;
    SelecChannel.UpperBound = kChannelToFill+1;
    SelectionVec.push_back(SelecChannel);
  }
 
  return Get1DVarHist(ProjectionVar_Str,SelectionVec,WeightStyle,Axis);
}

Get2DVarHistByModeAndChannel()

TH2 * SampleHandlerFD::Get2DVarHistByModeAndChannel	(	const std::string &	ProjectionVar_StrX,
		const std::string &	ProjectionVar_StrY,
		int	kModeToFill = -1,
		int	kChannelToFill = -1,
		int	WeightStyle = 0,
		TAxis *	AxisX = nullptr,
		TAxis *	AxisY = nullptr
	)

Definition at line 1778 of file SampleHandlerFD.cpp.

                                                                                      {
  bool fChannel;
  bool fMode;
 
  if (kChannelToFill!=-1) {
    if (kChannelToFill > GetNOscChannels()) {
      MACH3LOG_ERROR("Required channel is not available. kChannelToFill should be between 0 and {}", GetNOscChannels());
      MACH3LOG_ERROR("kChannelToFill given:{}", kChannelToFill);
      MACH3LOG_ERROR("Exiting.");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    fChannel = true;
  } else {
    fChannel = false;
  }
 
  if (kModeToFill!=-1) {
    if (!(kModeToFill >= 0) && (kModeToFill < Modes->GetNModes())) {
      MACH3LOG_ERROR("Required mode is not available. kModeToFill should be between 0 and {}", Modes->GetNModes());
      MACH3LOG_ERROR("kModeToFill given:{}", kModeToFill);
      MACH3LOG_ERROR("Exiting..");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    fMode = true;
  } else {
    fMode = false;
  }
 
  std::vector< KinematicCut > SelectionVec;
 
  if (fMode) {
    KinematicCut SelecMode;
    SelecMode.ParamToCutOnIt = ReturnKinematicParameterFromString("Mode");
    SelecMode.LowerBound = kModeToFill;
    SelecMode.UpperBound = kModeToFill+1;
    SelectionVec.push_back(SelecMode);
  }
 
  if (fChannel) {
    KinematicCut SelecChannel;
    SelecChannel.ParamToCutOnIt = ReturnKinematicParameterFromString("OscillationChannel");
    SelecChannel.LowerBound = kChannelToFill;
    SelecChannel.UpperBound = kChannelToFill+1;
    SelectionVec.push_back(SelecChannel);
  }
 
  return Get2DVarHist(ProjectionVar_StrX,ProjectionVar_StrY,SelectionVec,WeightStyle,AxisX,AxisY);
}

GetDataHist()

TH1 * SampleHandlerFD::GetDataHist

(

const int

Dimension

)

Get Data histogram.

Definition at line 1049 of file SampleHandlerFD.cpp.

                                                     {
// ************************************************
  if(Dimension == 1) {
    return SampleDetails.dathist;
  } else if(Dimension == 2) {
    return SampleDetails.dathist2d;
  } else{
    MACH3LOG_ERROR("Asdking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

GetFlavourName()

std::string SampleHandlerFD::GetFlavourName ( const int  iChannel )

inline

Definition at line 83 of file SampleHandlerFD.h.

                                               {
    if (iChannel < 0 || iChannel > GetNOscChannels()) {
      MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
      throw MaCh3Exception(__FILE__ , __LINE__);      
    }
    return OscChannels[iChannel].flavourName;
  }

GetLikelihood()

virtual double SampleHandlerBase::GetLikelihood ( )

pure virtual

Implemented in SampleHandlerFD, and PySampleHandlerBase.

GetMaCh3Modes()

MaCh3Modes* SampleHandlerBase::GetMaCh3Modes ( ) const

inline

Return pointer to MaCh3 modes.

Definition at line 53 of file SampleHandlerBase.h.

{ return Modes.get(); }

GetMCHist()

TH1 * SampleHandlerFD::GetMCHist

(

const int

Dimension

)

Get MC histogram.

Definition at line 1034 of file SampleHandlerFD.cpp.

                                                   {
// ************************************************
  FillMCHist(Dimension);
 
  if(Dimension == 1) {
    return SampleDetails._hPDF1D;
  } else if(Dimension == 2) {
    return SampleDetails._hPDF2D;
  } else{
    MACH3LOG_ERROR("Asking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

GetModeHist1D()

TH1* SampleHandlerFD::GetModeHist1D ( int  s, int  m, int  style = 0  )

inline

Definition at line 111 of file SampleHandlerFD.h.

                                                  {
    return Get1DVarHistByModeAndChannel(GetXBinVarName(),m,s,style);
  }

GetModeHist2D()

TH2* SampleHandlerFD::GetModeHist2D ( int  s, int  m, int  style = 0  )

inline

Definition at line 115 of file SampleHandlerFD.h.

                                                  {
    return Get2DVarHistByModeAndChannel(GetXBinVarName(),GetYBinVarName(),m,s,style);
  }

GetNDim()

int SampleHandlerFD::GetNDim ( ) const

inline

DB Function to differentiate 1D or 2D binning.

Definition at line 31 of file SampleHandlerFD.h.

{ return SampleDetails.nDimensions; }

GetNEvents()

unsigned int SampleHandlerBase::GetNEvents ( ) const

inline

Definition at line 60 of file SampleHandlerBase.h.

{return nEvents;}

GetNOscChannels() [1/2]

virtual int SampleHandlerBase::GetNOscChannels ( )

inlinevirtual

Reimplemented in SampleHandlerFD.

Definition at line 62 of file SampleHandlerBase.h.

{ return 1; }

GetNOscChannels() [2/2]

int SampleHandlerFD::GetNOscChannels ( )

inlineoverridevirtual

Reimplemented from SampleHandlerBase.

Definition at line 80 of file SampleHandlerFD.h.

{return static_cast<int>(OscChannels.size());}

GetNsamples()

virtual M3::int_t SampleHandlerBase::GetNsamples ( )

inlinevirtual

Definition at line 40 of file SampleHandlerBase.h.

{ return nSamples; };

GetPoissonLLH()

double SampleHandlerBase::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 18 of file SampleHandlerBase.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 SampleHandlerBase::GetSampleLikelihood ( const int  isample )

inlinevirtual

Definition at line 46 of file SampleHandlerBase.h.

{(void) isample; return GetLikelihood();};

GetSampleName() [1/2]

std::string SampleHandlerFD::GetSampleName ( int  iSample = 0 ) const

overridevirtual

Implements SampleHandlerBase.

Definition at line 1273 of file SampleHandlerFD.cpp.

                                                          {
  //ETA - this is just to suppress a warning for an unused variable
  (void)iSample;
 
  //ETA - extra safety to make sure SampleName is actually set
  // probably unnecessary due to the requirement for it to be in the yaml config
  if(SampleName.length() == 0){
    MACH3LOG_ERROR("No sample name provided");
    MACH3LOG_ERROR("Please provide a SampleName in your configuration file: {}", SampleManager->GetFileName());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  return SampleName;
}

GetSampleName() [2/2]

virtual std::string SampleHandlerBase::GetSampleName ( int  Sample ) const

pure virtual

Implemented in SampleHandlerFD, and PySampleHandlerBase.

GetTestStatLLH()

double SampleHandlerBase::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 [24] (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 36 of file SampleHandlerBase.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
}

GetTitle() [1/2]

virtual std::string SampleHandlerBase::GetTitle ( ) const

inlinevirtual

Reimplemented in SampleHandlerFD.

Definition at line 42 of file SampleHandlerBase.h.

{return "SampleHandler";};

GetTitle() [2/2]

std::string SampleHandlerFD::GetTitle ( ) const

inlineoverridevirtual

Reimplemented from SampleHandlerBase.

Definition at line 35 of file SampleHandlerFD.h.

{return SampleDetails.SampleTitle;}

GetW2Hist()

TH1 * SampleHandlerFD::GetW2Hist

(

const int

Dimension

)

Get W2 histogram.

Definition at line 996 of file SampleHandlerFD.cpp.

                                                   {
// ************************************************
  if(Dimension == 1) {
    TH1D* W2Hist = dynamic_cast<TH1D*>(SampleDetails._hPDF1D->Clone((SampleDetails._hPDF1D->GetName() + std::string("_W2")).c_str()));
    if (!W2Hist) {
      MACH3LOG_ERROR("Failed to cast");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    W2Hist->Reset();
    for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
      for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
        const int idx = Binning.GetBinSafe(xBin, yBin);
        W2Hist->AddBinContent(idx + 1, SampleHandlerFD_array_w2[idx]);
      }
    }
    return W2Hist;
  } else if(Dimension == 2) {
    TH2D* W2Hist = dynamic_cast<TH2D*>(SampleDetails._hPDF2D->Clone((SampleDetails._hPDF2D->GetName() + std::string("_W2")).c_str()));
    if (!W2Hist) {
      MACH3LOG_ERROR("Failed to cast");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    W2Hist->Reset();
    for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
      for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
        const int idx = Binning.GetBinSafe(xBin, yBin);
        W2Hist->SetBinContent(static_cast<int>(xBin + 1), static_cast<int>(yBin + 1), SampleHandlerFD_array_w2[idx]);
      }
    }
    return W2Hist;
  } else{
    MACH3LOG_ERROR("Asking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

GetXBinVarName()

std::string SampleHandlerFD::GetXBinVarName ( ) const

inline

Definition at line 38 of file SampleHandlerFD.h.

{return SampleDetails.XVarStr;}

GetYBinVarName()

std::string SampleHandlerFD::GetYBinVarName ( ) const

inline

Definition at line 40 of file SampleHandlerFD.h.

{return SampleDetails.YVarStr;}

ReturnHistsBySelection1D()

std::vector< TH1 * > SampleHandlerFD::ReturnHistsBySelection1D	(	std::string	KinematicProjection,
		int	Selection1,
		int	Selection2 = -1,
		int	WeightStyle = 0,
		TAxis *	Axis = 0
	)

Definition at line 1970 of file SampleHandlerFD.cpp.

                                                                                                                                                      {
  std::vector<TH1*> hHistList;
  std::string legendEntry;
 
  if (THStackLeg != nullptr) {delete THStackLeg;}
  THStackLeg = new TLegend(0.1,0.1,0.9,0.9);
 
  int iMax = -1;
  if (Selection1 == FDPlotType::kModePlot) {
    iMax = Modes->GetNModes();
  }
  if (Selection1 == FDPlotType::kOscChannelPlot) {
    iMax = GetNOscChannels();
  }
  if (iMax == -1) {
    MACH3LOG_ERROR("You've passed me a Selection1 which was not implemented in ReturnHistsBySelection1D. Selection1 and Selection2 are counters for different indexable quantities");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  for (int i=0;i<iMax;i++) {
    if (Selection1 == FDPlotType::kModePlot) {
      hHistList.push_back(Get1DVarHistByModeAndChannel(KinematicProjection,i,Selection2,WeightStyle,XAxis));
      THStackLeg->AddEntry(hHistList[i],(Modes->GetMaCh3ModeName(i)+Form(" : (%4.2f)",hHistList[i]->Integral())).c_str(),"f");
 
      hHistList[i]->SetFillColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
      hHistList[i]->SetLineColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
    }
    if (Selection1 == FDPlotType::kOscChannelPlot) {
      hHistList.push_back(Get1DVarHistByModeAndChannel(KinematicProjection,Selection2,i,WeightStyle,XAxis));
      THStackLeg->AddEntry(hHistList[i],(GetFlavourName(i)+Form(" | %4.2f",hHistList[i]->Integral())).c_str(),"f");
    }
  }
 
  return hHistList;
}

ReturnHistsBySelection2D()

std::vector< TH2 * > SampleHandlerFD::ReturnHistsBySelection2D	(	std::string	KinematicProjectionX,
		std::string	KinematicProjectionY,
		int	Selection1,
		int	Selection2 = -1,
		int	WeightStyle = 0,
		TAxis *	XAxis = 0,
		TAxis *	YAxis = 0
	)

Definition at line 2006 of file SampleHandlerFD.cpp.

                                                                                        {
// ************************************************
  std::vector<TH2*> hHistList;
 
  int iMax = -1;
  if (Selection1 == FDPlotType::kModePlot) {
    iMax = Modes->GetNModes();
  }
  if (Selection1 == FDPlotType::kOscChannelPlot) {
    iMax = GetNOscChannels();
  }
  if (iMax == -1) {
    MACH3LOG_ERROR("You've passed me a Selection1 which was not implemented in ReturnHistsBySelection1D. Selection1 and Selection2 are counters for different indexable quantities");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  for (int i=0;i<iMax;i++) {
    if (Selection1 == FDPlotType::kModePlot) {
      hHistList.push_back(Get2DVarHistByModeAndChannel(KinematicProjectionX,KinematicProjectionY,i,Selection2,WeightStyle,XAxis,YAxis));
    }
    if (Selection1 == FDPlotType::kOscChannelPlot) {
      hHistList.push_back(Get2DVarHistByModeAndChannel(KinematicProjectionX,KinematicProjectionY,Selection2,i,WeightStyle,XAxis,YAxis));
    }
  }
 
  return hHistList;
}

ReturnKinematicParameterFromString()

int SampleHandlerFD::ReturnKinematicParameterFromString

(

const std::string &

KinematicStr

)

const

ETA function to generically convert a string from xsec cov to a kinematic type.

Definition at line 1627 of file SampleHandlerFD.cpp.

                                                                                                    {
// ************************************************
  auto it = KinematicParameters->find(KinematicParameterStr);
  if (it != KinematicParameters->end()) return it->second;
 
  MACH3LOG_ERROR("Did not recognise Kinematic Parameter type: {}", KinematicParameterStr);
  throw MaCh3Exception(__FILE__, __LINE__);
 
  return M3::_BAD_INT_;
}

ReturnStackedHistBySelection1D()

THStack * SampleHandlerFD::ReturnStackedHistBySelection1D	(	std::string	KinematicProjection,
		int	Selection1,
		int	Selection2 = -1,
		int	WeightStyle = 0,
		TAxis *	Axis = 0
	)

Definition at line 2036 of file SampleHandlerFD.cpp.

                                                                                                                                                     {
  std::vector<TH1*> HistList = ReturnHistsBySelection1D(KinematicProjection, Selection1, Selection2, WeightStyle, XAxis);
  THStack* StackHist = new THStack((GetTitle()+"_"+KinematicProjection+"_Stack").c_str(),"");
  for (unsigned int i=0;i<HistList.size();i++) {
    StackHist->Add(HistList[i]);
  }
  return StackHist;
}

ReturnStackHistLegend()

TLegend* SampleHandlerFD::ReturnStackHistLegend ( )

inline

Definition at line 126 of file SampleHandlerFD.h.

{return THStackLeg;}

ReturnStringFromKinematicParameter()

std::string SampleHandlerFD::ReturnStringFromKinematicParameter

(

const int

KinematicVariable

)

const

ETA function to generically convert a kinematic type from xsec cov to a string.

Definition at line 1639 of file SampleHandlerFD.cpp.

                                                                                                {
// ************************************************
  auto it = ReversedKinematicParameters->find(KinematicParameter);
  if (it != ReversedKinematicParameters->end()) {
    return it->second;
  }
 
  MACH3LOG_ERROR("Did not recognise Kinematic Parameter type: {}", KinematicParameter);
  throw MaCh3Exception(__FILE__, __LINE__);
 
  return "";
}

SetTestStatistic()

void SampleHandlerBase::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 172 of file SampleHandlerBase.h.

{ fTestStatistic = testStat; }