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.
virtual double	SampleHandlerBase::GetLikelihood ()=0
unsigned int	SampleHandlerBase::GetNEvents () const
virtual int	SampleHandlerBase::GetNMCSamples ()
virtual int	SampleHandlerBase::GetNOscChannels ()
double	SampleHandlerBase::GetTestStatLLH (double data, double mc) const
	Calculate test statistic for a single bin using Poisson.
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.
void	SampleHandlerBase::SetTestStatistic (TestStatistic testStat)
	Set the test statistic to be used when calculating the binned likelihoods.
int	SampleHandlerFD::GetNDim () const
	DB Function to differentiate 1D or 2D binning.
std::string	SampleHandlerFD::GetSampleName (int iSample=0) const override
std::string	SampleHandlerFD::GetTitle () const override
std::string	SampleHandlerFD::GetXBinVarName ()
std::string	SampleHandlerFD::GetYBinVarName ()
TH1 *	SampleHandlerFD::GetMCHist (const int Dimension)
	Get MC histogram.
TH1 *	SampleHandlerFD::GetW2Hist (const int Dimension)
	Get W2 histogram.
TH1 *	SampleHandlerFD::GetDataHist (const int Dimension)
	Get Data histogram.
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.
std::string	SampleHandlerFD::ReturnStringFromKinematicParameter (const int KinematicVariable) const
	ETA function to generically convert a kinematic type from xsec cov to a string.

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 1553 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 1500 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());
  }
  
  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 1790 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 1840 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 1136 of file SampleHandlerFD.cpp.

                                                     {
// ************************************************
  if(Dimension == 1) {
    return dathist;
  } else if(Dimension == 2) {
    return 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 > static_cast<int>(OscChannels.size())) {
      MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
      throw MaCh3Exception(__FILE__ , __LINE__);      
    }
    return OscChannels[iChannel].flavourName;
  }

GetLikelihood()

virtual double SampleHandlerBase::GetLikelihood ( )

pure virtual

Implemented in PySampleHandlerBase, and SampleHandlerFD.

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 1121 of file SampleHandlerFD.cpp.

                                                   {
// ************************************************
  if(Dimension == 1) {
    Fill1DHist();
    return _hPDF1D;
  } else if(Dimension == 2) {
    Fill2DHist();
    return _hPDF2D;
  } else{
    MACH3LOG_ERROR("Asdking 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(XVarStr,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(XVarStr,YVarStr,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 nDimensions; }

GetNEvents()

unsigned int SampleHandlerBase::GetNEvents ( ) const

inline

Definition at line 60 of file SampleHandlerBase.h.

{return nEvents;}

GetNMCSamples()

virtual int SampleHandlerBase::GetNMCSamples ( )

inlinevirtual

Definition at line 62 of file SampleHandlerBase.h.

{ return nSamples; }

GetNOscChannels() [1/2]

virtual int SampleHandlerBase::GetNOscChannels ( )

inlinevirtual

Reimplemented in SampleHandlerFD.

Definition at line 64 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; };

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 1341 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 PySampleHandlerBase, and SampleHandlerFD.

GetTestStatLLH() [1/2]

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.

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 is Sum(w_{i}^2) (sum of weights squared), which is sigma^2_{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 GetTestStatLLH(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 GetTestStatLLH(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 = GetTestStatLLH(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 GetTestStatLLH which doesn't take in weights
      //and is a Poisson likelihood comparison.
      return GetTestStatLLH(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
}

GetTestStatLLH() [2/2]

double SampleHandlerBase::GetTestStatLLH	(	double	data,
		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 ) );
}

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 SampleTitle;}

GetW2Hist()

TH1 * SampleHandlerFD::GetW2Hist

(

const int

Dimension

)

Get W2 histogram.

Definition at line 1094 of file SampleHandlerFD.cpp.

                                                   {
// ************************************************
  if(Dimension == 1) {
    TH1D* W2Hist = dynamic_cast<TH1D*>(_hPDF1D->Clone((_hPDF1D->GetName() + std::string("_W2")).c_str()));
    W2Hist->Reset();
    for (unsigned int yBin = 0; yBin < (Binning.YBinEdges.size()-1); yBin++) {
      for (unsigned int xBin = 0; xBin < (Binning.XBinEdges.size()-1); xBin++) {
        W2Hist->AddBinContent(xBin+1, SampleHandlerFD_array_w2[yBin][xBin]);
      }
    }
    return W2Hist;
  } else if(Dimension == 2) {
    TH2D* W2Hist = dynamic_cast<TH2D*>(_hPDF2D->Clone((_hPDF2D->GetName() + std::string("_W2")).c_str()));
    W2Hist->Reset();
    for (unsigned int yBin = 0; yBin < (Binning.YBinEdges.size()-1); yBin++) {
      for (unsigned int xBin = 0; xBin < (Binning.XBinEdges.size()-1); xBin++) {
        W2Hist->SetBinContent(xBin+1, yBin+1, SampleHandlerFD_array_w2[yBin][xBin]);
      }
    }
    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 ( )

inline

Definition at line 38 of file SampleHandlerFD.h.

{return XVarStr;}

GetYBinVarName()

std::string SampleHandlerFD::GetYBinVarName ( )

inline

Definition at line 40 of file SampleHandlerFD.h.

{return YVarStr;}

ReturnHistsBySelection1D()

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

Definition at line 2030 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 == 0) {
    iMax = Modes->GetNModes();
  }
  if (Selection1 == 1) {
    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==0) {
      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==1) {
      hHistList.push_back(Get1DVarHistByModeAndChannel(KinematicProjection,Selection2,i,WeightStyle,XAxis));
      THStackLeg->AddEntry(hHistList[i],(OscChannels[i].flavourName+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 2066 of file SampleHandlerFD.cpp.

                                {
  std::vector<TH2*> hHistList;
 
  int iMax = -1;
  if (Selection1 == 0) {
    iMax = Modes->GetNModes();
  }
  if (Selection1 == 1) {
    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==0) {
      hHistList.push_back(Get2DVarHistByModeAndChannel(KinematicProjectionX,KinematicProjectionY,i,Selection2,WeightStyle,XAxis,YAxis));
    }
    if (Selection1==1) {
      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 1689 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 2095 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 1701 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 99 of file SampleHandlerBase.h.

{ fTestStatistic = testStat; }