SampleHandlerFD Class Reference

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

#include <Samples/SampleHandlerFD.h>

Inheritance diagram for SampleHandlerFD:

Collaboration diagram for SampleHandlerFD:

Public Member Functions
	SampleHandlerFD (std::string ConfigFileName, ParameterHandlerGeneric *xsec_cov, const std::shared_ptr< OscillationHandler > &OscillatorObj_=nullptr)
	Constructor.
virtual	~SampleHandlerFD ()
	destructor
int	GetNDim () const
	DB Function to differentiate 1D or 2D binning.
std::string	GetSampleName (int iSample=0) const override
std::string	GetTitle () const override
std::string	GetXBinVarName ()
std::string	GetYBinVarName ()
void	PrintIntegral (TString OutputName="/dev/null", int WeightStyle=0, TString OutputCSVName="/dev/null")
void	AddData (TH1D *Data)
void	AddData (TH2D *Data)
void	AddData (std::vector< double > &data)
void	AddData (std::vector< std::vector< double > > &data)
double	GetLikelihood () override
	DB Multi-threaded GetLikelihood.
TH1 *	GetMCHist (const int Dimension)
	Get MC histogram.
TH1 *	GetW2Hist (const int Dimension)
	Get W2 histogram.
TH1 *	GetDataHist (const int Dimension)
	Get Data histogram.
void	Reweight () override
M3::float_t	GetEventWeight (const int iEntry) const
void	InitialiseNuOscillatorObjects ()
	including Dan's magic NuOscillator
void	SetupNuOscillatorPointers ()
void	ReadSampleConfig ()
int	GetNOscChannels () override
std::string	GetFlavourName (const int iChannel)
TH1 *	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 >())
TH2 *	Get2DVarHist (const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &EventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >())
void	Fill1DSubEventHist (TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
void	Fill2DSubEventHist (TH2D *_h2DVar, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
TH1 *	Get1DVarHistByModeAndChannel (const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr)
TH2 *	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 *	GetModeHist1D (int s, int m, int style=0)
TH2 *	GetModeHist2D (int s, int m, int style=0)
std::vector< TH1 * >	ReturnHistsBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
std::vector< TH2 * >	ReturnHistsBySelection2D (std::string KinematicProjectionX, std::string KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *XAxis=0, TAxis *YAxis=0)
THStack *	ReturnStackedHistBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
TLegend *	ReturnStackHistLegend ()
int	ReturnKinematicParameterFromString (const std::string &KinematicStr) const
	ETA function to generically convert a string from xsec cov to a kinematic type.
std::string	ReturnStringFromKinematicParameter (const int KinematicVariable) const
	ETA function to generically convert a kinematic type from xsec cov to a string.
void	SaveAdditionalInfo (TDirectory *Dir) override
	Store additional info in a chan.
int	ReturnKinematicVectorFromString (const std::string &KinematicStr) const
std::string	ReturnStringFromKinematicVector (const int KinematicVariable) const
bool	IsSubEventVarString (const std::string &VarStr)
	JM Check if a kinematic parameter string corresponds to a subevent-level variable.
Public Member Functions inherited from SampleHandlerBase
	SampleHandlerBase ()
	The main constructor.
virtual	~SampleHandlerBase ()
	destructor
virtual M3::int_t	GetNsamples ()
virtual std::string	GetTitle () const
virtual std::string	GetSampleName (int Sample) const =0
virtual double	GetSampleLikelihood (const int isample)
virtual void	CleanMemoryBeforeFit ()=0
	Allow to clean not used memory before fit starts.
virtual void	SaveAdditionalInfo (TDirectory *Dir)
	Store additional info in a chan.
MaCh3Modes *	GetMaCh3Modes () const
	Return pointer to MaCh3 modes.
virtual void	Reweight ()=0
virtual double	GetLikelihood ()=0
unsigned int	GetNEvents () const
virtual int	GetNMCSamples ()
virtual int	GetNOscChannels ()
virtual void	SetupBinning (const M3::int_t Selection, std::vector< double > &BinningX, std::vector< double > &BinningY)
virtual TH1 *	GetData (const int Selection)
virtual TH2Poly *	GetW2 (const int Selection)
virtual TH1 *	GetPDF (const int Selection)
virtual TH1 *	GetPDFMode (const int Selection, const int Mode)
virtual std::string	GetKinVarLabel (const int sample, const int Dimension)
double	GetTestStatLLH (double data, double mc) const
	Calculate test statistic for a single bin using Poisson.
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.
void	SetTestStatistic (TestStatistic testStat)
	Set the test statistic to be used when calculating the binned likelihoods.

Protected Member Functions
virtual void	SetupWeightPointers ()=0
	DB Function to determine which weights apply to which types of samples pure virtual!!
void	SetupKinematicMap ()
	Ensure Kinematic Map is setup and make sure it is initialised correctly.
virtual void	SetupSplines ()=0
	initialise your splineXX object and then use InitialiseSplineObject to conviently setup everything up
virtual void	Init ()=0
	Initialise any variables that your experiment specific SampleHandler needs.
virtual int	SetupExperimentMC ()=0
	Experiment specific setup, returns the number of events which were loaded.
virtual void	SetupFDMC ()=0
	Function which translates experiment struct into core struct.
void	Initialise ()
	Function which does a lot of the lifting regarding the workflow in creating different MC objects.
void	FillSplineBins ()
	Finds the binned spline that an event should apply to and stored them in a a vector for easy evaluation in the fillArray() function.
void	FindNominalBinAndEdges1D ()
void	FindNominalBinAndEdges2D ()
void	Set1DBinning (size_t nbins, double *boundaries)
	sets the binning used for the likelihood calculation, used for both data and MC
void	Set2DBinning (size_t nbins1, double *boundaries1, size_t nbins2, double *boundaries2)
	set the binning for 2D sample used for the likelihood calculation
void	Set1DBinning (std::vector< double > &XVec)
	set the binning for 1D sample used for the likelihood calculation
void	Set2DBinning (std::vector< double > &XVec, std::vector< double > &YVec)
	set the binning for 2D sample used for the likelihood calculation
void	SetupSampleBinning ()
	wrapper to call set binning functions based on sample config info
void	SetupReweightArrays (const size_t numberXBins, const size_t numberYBins)
	Initialise data, MC and W2 histograms.
virtual void	SetupFunctionalParameters ()
	ETA - a function to setup and pass values to functional parameters where you need to pass a value to some custom reweight calc or engine.
void	RegisterIndividualFunctionalParameter (const std::string &fpName, int fpEnum, FuncParFuncType fpFunc)
	HH - a helper function for RegisterFunctionalParameter.
virtual void	RegisterFunctionalParameters ()=0
	HH - a experiment-specific function where the maps to actual functions are set up.
virtual void	PrepFunctionalParameters ()
	Update the functional parameter values to the latest proposed values. Needs to be called before every new reweight so is called in fillArray.
virtual void	ApplyShifts (int iEvent)
	ETA - generic function applying shifts.
bool	IsEventSelected (const int iEvent)
	DB Function which determines if an event is selected, where Selection double looks like {{ND280KinematicTypes Var1, douuble LowBound}.
bool	IsSubEventSelected (const std::vector< KinematicCut > &SubEventCuts, const int iEvent, unsigned const int iSubEvent, size_t nsubevents)
	JM Function which determines if a subevent is selected.
virtual void	resetShifts (int iEvent)
	HH - reset the shifted values to the original values.
void	CalcNormsBins (std::vector< NormParameter > &norm_parameters, std::vector< std::vector< int > > &xsec_norms_bins)
	Check whether a normalisation systematic affects an event or not.
M3::float_t	CalcWeightSpline (const FarDetectorCoreInfo *MCEvent) const
	Calculate the spline weight for a given event.
M3::float_t	CalcWeightNorm (const FarDetectorCoreInfo *MCEvent) const
	Calculate the norm weight for a given event.
virtual void	CalcWeightFunc (int iEvent)
	Calculate weights for function parameters.
virtual double	ReturnKinematicParameter (std::string KinematicParamter, int iEvent)=0
	Return the value of an assocaited kinematic parameter for an event.
virtual double	ReturnKinematicParameter (int KinematicVariable, int iEvent)=0
virtual std::vector< double >	ReturnKinematicVector (std::string KinematicParameter, int iEvent)
virtual std::vector< double >	ReturnKinematicVector (int KinematicVariable, int iEvent)
std::vector< double >	ReturnKinematicParameterBinning (const std::string &KinematicParameter)
	Return the binning used to draw a kinematic parameter.
virtual const double *	GetPointerToKinematicParameter (std::string KinematicParamter, int iEvent)=0
virtual const double *	GetPointerToKinematicParameter (double KinematicVariable, int iEvent)=0
const double *	GetPointerToOscChannel (const int iEvent) const
	Get pointer to oscillation channel associated with given event. Osc channel is const.
void	SetupNormParameters ()
	Setup the norm parameters by assigning each event with bin.
void	Fill1DHist ()
	Fill a 1D histogram with the event-level information used in the fit.
void	Fill2DHist ()
	Fill a 2D histogram with the event-level information used in the fit.
void	FillArray ()
	DB Nice new multi-threaded function which calculates the event weights and fills the relevant bins of an array.
void	ResetHistograms ()
	Helper function to reset histograms.
void	InitialiseSingleFDMCObject ()
	function to create the member of the FarDetectorInfo struct so they are the appropriate size.
void	InitialiseSplineObject ()
NuPDG	GetInitPDGFromFileName (const std::string &FileName) const
	Retrieve the initial neutrino PDG code associated with a given input file name.
NuPDG	GetFinalPDGFromFileName (const std::string &FileName) const
	Retrieve the final neutrino PDG code associated with a given input file name.
Protected Member Functions inherited from SampleHandlerBase
void	QuietPlease ()
	CW: Redirect std::cout to silence some experiment specific libraries.
void	NowTalk ()
	CW: Redirect std::cout to silence some experiment specific libraries.
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

Protected Attributes
std::unique_ptr< BinnedSplineHandler >	SplineHandler
	Contains all your binned splines and handles the setup and the returning of weights from spline evaluations.
std::shared_ptr< OscillationHandler >	Oscillator
	Contains all your binned splines and handles the setup and the returning of weights from spline evaluations.
std::string	XVarStr
	the strings associated with the variables used for the binning e.g. "RecoNeutrinoEnergy"
std::string	YVarStr
std::vector< FunctionalParameter >	funcParsVec
	HH - a vector that stores all the FuncPars struct.
std::unordered_map< std::string, int >	funcParsNamesMap
	HH - a map that relates the name of the functional parameter to funcpar enum.
std::vector< FunctionalParameter * >	funcParsMap
	HH - a map that relates the funcpar enum to pointer of FuncPars struct HH - Changed to a vector of pointers since it's faster than unordered_map and we are using ints as keys.
std::unordered_map< int, FuncParFuncType >	funcParsFuncMap
	HH - a map that relates the funcpar enum to pointer of the actual function.
std::vector< std::vector< int > >	funcParsGrid
	HH - a grid of vectors of enums for each sample and event.
std::vector< std::string >	funcParsNamesVec = {}
	HH - a vector of string names for each functional parameter.
SampleBinningInfo	Binning
	KS: This stores binning information, in future could be come vector to store binning for every used sample.
double **	SampleHandlerFD_array
	DB Array to be filled after reweighting.
double **	SampleHandlerFD_array_w2
	KS Array used for MC stat.
double **	SampleHandlerFD_data
	DB Array to be filled in AddData.
std::vector< FarDetectorCoreInfo >	MCSamples
std::vector< OscChannelInfo >	OscChannels
ParameterHandlerGeneric *	ParHandler = nullptr
int	nDimensions = M3::_BAD_INT_
	Keep track of the dimensions of the sample binning.
std::string	SampleName
	A unique ID for each sample based on powers of two for quick binary operator comparisons.
std::string	SampleTitle
	the name of this sample e.g."muon-like"
std::vector< KinematicCut >	StoredSelection
	What gets pulled from config options, these are constant after loading in this is of length 3: 0th index is the value, 1st is lower bound, 2nd is upper bound.
std::vector< KinematicCut >	Selection
	a way to store selection cuts which you may push back in the get1DVar functions most of the time this is just the same as StoredSelection
const std::unordered_map< std::string, int > *	KinematicParameters
	Mapping between string and kinematic enum.
const std::unordered_map< int, std::string > *	ReversedKinematicParameters
	Mapping between kinematic enum and string.
const std::unordered_map< std::string, int > *	KinematicVectors
const std::unordered_map< int, std::string > *	ReversedKinematicVectors
std::unique_ptr< manager >	SampleManager
	The manager object used to read the sample yaml file.
std::vector< std::string >	mc_files
	names of mc files associated associated with this object
std::vector< std::string >	spline_files
	names of spline files associated associated with this object
std::unordered_map< std::string, double >	_modeNomWeightMap
TLegend *	THStackLeg = nullptr
	DB Miscellaneous Variables.
bool	FirstTimeW2
	KS:Super hacky to update W2 or not.
bool	UpdateW2
	KS:Super hacky to update W2 or not.
TH1D *	dathist
TH2D *	dathist2d
TH1D *	_hPDF1D
TH2D *	_hPDF2D
Protected Attributes inherited from SampleHandlerBase
TestStatistic	fTestStatistic
	Test statistic tells what kind of likelihood sample is using.
std::streambuf *	buf
	Keep the cout buffer.
std::streambuf *	errbuf
	Keep the cerr buffer.
M3::int_t	nSamples
	Contains how many samples we've got.
unsigned int	nEvents
	Number of MC events are there.
std::unique_ptr< MaCh3Modes >	Modes
	Holds information about used Generator and MaCh3 modes.

Private Attributes
std::unordered_map< std::string, NuPDG >	FileToInitPDGMap
std::unordered_map< std::string, NuPDG >	FileToFinalPDGMap

Detailed Description

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

Author

Dan Barrow

Ed Atkin

Definition at line 18 of file SampleHandlerFD.h.

Constructor & Destructor Documentation

SampleHandlerFD()

SampleHandlerFD::SampleHandlerFD	(	std::string	ConfigFileName,
		ParameterHandlerGeneric *	xsec_cov,
		const std::shared_ptr< OscillationHandler > &	OscillatorObj_ = nullptr
	)

Constructor.

Parameters

ConfigFileName

Name of config to initialise the sample object

Definition at line 10 of file SampleHandlerFD.cpp.

                                                                                          : SampleHandlerBase() {
// ************************************************
  MACH3LOG_INFO("-------------------------------------------------------------------");
  MACH3LOG_INFO("Creating SampleHandlerFD object");
 
  //ETA - safety feature so you can't pass a NULL xsec_cov
  if(!xsec_cov){
    MACH3LOG_ERROR("You've passed me a nullptr to a SystematicHandlerGeneric... I need this to setup splines!");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  ParHandler = xsec_cov;
 
  nSamples = 1;
 
  if (OscillatorObj_ != nullptr) {
    MACH3LOG_WARN("You have passed an Oscillator object through the constructor of a SampleHandlerFD object - this will be used for all oscillation channels");
    Oscillator = OscillatorObj_;
  }
  dathist = nullptr;
  dathist2d = nullptr;
  _hPDF1D = nullptr;
  _hPDF2D = nullptr;
 
  KinematicParameters = nullptr;
  ReversedKinematicParameters = nullptr;
  KinematicVectors = nullptr;
  ReversedKinematicVectors = nullptr;
 
  SampleHandlerFD_array = nullptr;
  SampleHandlerFD_data = nullptr;
  SampleHandlerFD_array_w2 = nullptr;
  SampleName = "";
  SampleManager = std::unique_ptr<manager>(new manager(ConfigFileName.c_str()));
 
  // Variables related to MC stat
  FirstTimeW2 = true;
  UpdateW2 = false;
}

~SampleHandlerFD()

SampleHandlerFD::~SampleHandlerFD ( )

virtual

destructor

Definition at line 50 of file SampleHandlerFD.cpp.

                                  {
  MACH3LOG_DEBUG("I'm deleting SampleHandlerFD");
  
  for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
    if(SampleHandlerFD_array != nullptr){delete[] SampleHandlerFD_array[yBin];}
    delete[] SampleHandlerFD_array_w2[yBin];
    //ETA - there is a chance that you haven't added any data...
    if(SampleHandlerFD_data != nullptr){delete[] SampleHandlerFD_data[yBin];}
  }
 
  if(SampleHandlerFD_array != nullptr){delete[] SampleHandlerFD_array;}
  delete[] SampleHandlerFD_array_w2;
  //ETA - there is a chance that you haven't added any data...
  if(SampleHandlerFD_data != nullptr){delete[] SampleHandlerFD_data;}
 
  if(THStackLeg != nullptr) delete THStackLeg;
 
  if(dathist != nullptr) delete dathist;
  if(dathist2d != nullptr) delete dathist2d;
  if(_hPDF1D != nullptr) delete _hPDF1D;
  if(_hPDF2D != nullptr) delete _hPDF2D;
}

Member Function Documentation

AddData() [1/4]

void SampleHandlerFD::AddData

(

std::vector< double > &

data

)

Definition at line 1148 of file SampleHandlerFD.cpp.

                                                     {
  dathist2d = nullptr;
  dathist->Reset();
 
  if (GetNDim()!=1) {
    MACH3LOG_ERROR("Trying to set a 1D 'data' histogram when the number of dimensions for this sample is {}", GetNDim());
    MACH3LOG_ERROR("This won't work, please specify the correct dimensions in your sample config with the X and Y variables");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  for (auto const& data_point : data){
    dathist->Fill(data_point);
  }
 
  if(SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
    for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
      SampleHandlerFD_data[yBin][xBin] = dathist->GetBinContent(static_cast<int>(xBin+1));
    }
  }
}

AddData() [2/4]

void SampleHandlerFD::AddData

(

std::vector< std::vector< double > > &

data

)

Definition at line 1173 of file SampleHandlerFD.cpp.

                                                                   {
  dathist = nullptr;
  dathist2d->Reset();
 
  if (GetNDim()!=2) {
    MACH3LOG_ERROR("Trying to set a 2D 'data' histogram when the number of dimensions for this sample is {}", GetNDim());
    MACH3LOG_ERROR("This won't work, please specify the correct dimensions in your sample config with the X and Y variables");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  //TODO: this assumes that std::vector is of length 2 and then both data.at(0) and data
  //ETA: I think this might just be wrong? We should probably just make this AddData(std::vector<double> data_x, std::vector<double> data_y)
  // or maybe something like AddData(std::vector<std::pair<double, double>> data)?
  for (int i = 0; i < int(data.size()); i++) {
    dathist2d->Fill(data.at(0)[i],data.at(1)[i]);
  }
 
  if(SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
    for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
      //Need to cast to an int (Int_t) for ROOT
      //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
      SampleHandlerFD_data[yBin][xBin] = dathist2d->GetBinContent(static_cast<int>(xBin+1),static_cast<int>(yBin+1));
    }
  }
}

AddData() [3/4]

void SampleHandlerFD::AddData

(

TH1D *

Data

)

Definition at line 1203 of file SampleHandlerFD.cpp.

                                        {
  MACH3LOG_INFO("Adding 1D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
  dathist2d = nullptr;
  dathist = static_cast<TH1D*>(Data->Clone());
 
  if (GetNDim() != 1) {
    MACH3LOG_ERROR("Trying to set a 1D 'data' histogram in a 2D sample - Quitting"); 
    MACH3LOG_ERROR("The number of dimensions for this sample is {}", GetNDim());
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
    
  if(SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
    for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
      //Need to cast to an int (Int_t) for ROOT
      //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
      SampleHandlerFD_data[yBin][xBin] = Data->GetBinContent(static_cast<int>(xBin+1));
    }
  }
}

AddData() [4/4]

void SampleHandlerFD::AddData

(

TH2D *

Data

)

Definition at line 1227 of file SampleHandlerFD.cpp.

                                        {
  MACH3LOG_INFO("Adding 2D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
  dathist2d = static_cast<TH2D*>(Data->Clone());
  dathist = nullptr;
 
  if (GetNDim() != 2) {
    MACH3LOG_ERROR("Trying to set a 2D 'data' histogram in a 1D sample - Quitting"); 
    throw MaCh3Exception(__FILE__ , __LINE__ );}
   
  if(SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
    for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
      //Need to cast to an int (Int_t) for ROOT
      //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
      SampleHandlerFD_data[yBin][xBin] = dathist2d->GetBinContent(static_cast<int>(xBin+1),static_cast<int>(yBin+1));
    }
  }
}

ApplyShifts()

void SampleHandlerFD::ApplyShifts ( int  iEvent )

protectedvirtual

ETA - generic function applying shifts.

Definition at line 714 of file SampleHandlerFD.cpp.

                                            {
  // Given a sample and event, apply the shifts to the event based on the vector of functional parameter enums
  // First reset shifted array back to nominal values
  resetShifts(iEvent);
 
  const std::vector<int>& fpEnums = funcParsGrid[iEvent];
  const std::size_t nShifts = fpEnums.size();
 
  for (std::size_t i = 0; i < nShifts; ++i) {
    const int fpEnum = fpEnums[i];
    FunctionalParameter *fp = funcParsMap[static_cast<std::size_t>(fpEnum)];
    // if (fp->funcPtr) {
    //   (*fp->funcPtr)(fp->valuePtr, iEvent);
    // } else {
    //   MACH3LOG_ERROR("Functional parameter function pointer for {} is null for event {} in sample {}", fp->name, iEvent);
    //   throw MaCh3Exception(__FILE__, __LINE__);
    // }
    (*fp->funcPtr)(fp->valuePtr, iEvent);
  }
}

CalcNormsBins()

void SampleHandlerFD::CalcNormsBins ( std::vector< NormParameter > &  norm_parameters, std::vector< std::vector< int > > &  xsec_norms_bins  )

protected

Check whether a normalisation systematic affects an event or not.

Definition at line 806 of file SampleHandlerFD.cpp.

                                                                                                                             {
// ************************************************
  #ifdef DEBUG
  std::vector<int> VerboseCounter(norm_parameters.size(), 0);
  #endif
  for(unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent){
    std::vector< int > NormBins = {};
    if (ParHandler) {
      // Skip oscillated NC events
      // Not strictly needed, but these events don't get included in oscillated predictions, so
      // no need to waste our time calculating and storing information about xsec parameters
      // that will never be used.
      if (MCSamples[iEvent].isNC && (*MCSamples[iEvent].nupdg != *MCSamples[iEvent].nupdgUnosc) ) {
        MACH3LOG_TRACE("Event {}, missed NC/signal check", iEvent);
        continue;
      } //DB Abstract check on MaCh3Modes to determine which apply to neutral current
      for (std::vector<NormParameter>::iterator it = norm_parameters.begin(); it != norm_parameters.end(); ++it) {
        //Now check that the target of an interaction matches with the normalisation parameters
        bool TargetMatch = MatchCondition((*it).targets, *(MCSamples[iEvent].Target));
        if (!TargetMatch) {
          MACH3LOG_TRACE("Event {}, missed target check ({}) for dial {}", iEvent, *(MCSamples[iEvent].Target), (*it).name);
          continue;
        }
 
        //Now check that the neutrino flavour in an interaction matches with the normalisation parameters
        bool FlavourMatch = MatchCondition((*it).pdgs, *(MCSamples[iEvent].nupdg));
        if (!FlavourMatch) {
          MACH3LOG_TRACE("Event {}, missed PDG check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdg), (*it).name);
          continue;
        }
 
        //Now check that the unoscillated neutrino flavour in an interaction matches with the normalisation parameters
        bool FlavourUnoscMatch = MatchCondition((*it).preoscpdgs, *(MCSamples[iEvent].nupdgUnosc));
        if (!FlavourUnoscMatch){
          MACH3LOG_TRACE("Event {}, missed FlavourUnosc check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdgUnosc), (*it).name);
          continue;
        }
 
        //Now check that the mode of an interaction matches with the normalisation parameters
        bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
        if (!ModeMatch) {
          MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
          continue;
        }
 
        //Now check whether the norm has kinematic bounds
        //i.e. does it only apply to events in a particular kinematic region?
        // Now check whether within kinematic bounds
        bool IsSelected = true;
        if ((*it).hasKinBounds) {
          const auto& kinVars = (*it).KinematicVarStr;
          const auto& selection = (*it).Selection;
 
          for (std::size_t iKinPar = 0; iKinPar < kinVars.size(); ++iKinPar) {
            const double kinVal = ReturnKinematicParameter(kinVars[iKinPar], static_cast<int>(iEvent));
 
            bool passedAnyBound = false;
            const auto& boundsList = selection[iKinPar];
 
            for (const auto& bounds : boundsList) {
              if (kinVal > bounds[0] && kinVal <= bounds[1]) {
                passedAnyBound = true;
                break;
              }
            }
 
            if (!passedAnyBound) {
              MACH3LOG_TRACE("Event {}, missed kinematic check ({}) for dial {}",
                             iEvent, kinVars[iKinPar], (*it).name);
              IsSelected = false;
              break;
            }
          }
        }
        // Need to then break the event loop
        if(!IsSelected){
          MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", iEvent, (*it).name);
          continue;
        }
        // Now set 'index bin' for each normalisation parameter
        // All normalisations are just 1 bin for 2015, so bin = index (where index is just the bin for that normalisation)
        int bin = (*it).index;
 
        NormBins.push_back(bin);
        MACH3LOG_TRACE("Event {}, will be affected by dial {}", iEvent, (*it).name);
        #ifdef DEBUG
        VerboseCounter[std::distance(norm_parameters.begin(), it)]++;
        #endif
        //}
      } // end iteration over norm_parameters
    } // end if (ParHandler)
    xsec_norms_bins[iEvent] = NormBins;
  }//end loop over events
  #ifdef DEBUG
  MACH3LOG_DEBUG("┌──────────────────────────────────────────────────────────┐");
  for (std::size_t i = 0; i < norm_parameters.size(); ++i) {
    const auto& norm = norm_parameters[i];
    double eventRatio = static_cast<double>(VerboseCounter[i]) / static_cast<double>(GetNEvents());
 
    MACH3LOG_DEBUG("│ Param {:<15}, affects {:<8} events ({:>6.2f}%) │",
                  ParHandler->GetParFancyName(norm.index), VerboseCounter[i], eventRatio);
  }
  MACH3LOG_DEBUG("└──────────────────────────────────────────────────────────┘");
  #endif
}

CalcWeightFunc()

virtual void SampleHandlerFD::CalcWeightFunc ( int  iEvent )

inlineprotectedvirtual

Calculate weights for function parameters.

First you need to setup additional pointers in you experiment code in SetupWeightPointers Then in this function you can calculate whatever fancy function you want by filling weight to which you have pointer This way func weight shall be used in GetEventWeight

Definition at line 259 of file SampleHandlerFD.h.

{return; (void)iEvent;};

CalcWeightNorm()

M3::float_t SampleHandlerFD::CalcWeightNorm ( const FarDetectorCoreInfo *  MCEvent ) const

protected

Calculate the norm weight for a given event.

Definition at line 754 of file SampleHandlerFD.cpp.

                                                                                  {
// ***************************************************************************
  M3::float_t xsecw = 1.0;
  const int nNorms = static_cast<int>(MCEvent->xsec_norm_pointers.size());
  //Loop over stored normalisation and function pointers
  #ifdef MULTITHREAD
  #pragma omp simd
  #endif
  for (int iParam = 0; iParam < nNorms; ++iParam)
  {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wuseless-cast"
    xsecw *= static_cast<M3::float_t>(*(MCEvent->xsec_norm_pointers[iParam]));
#pragma GCC diagnostic pop
    #ifdef DEBUG
    if (std::isnan(xsecw)) MACH3LOG_WARN("iParam= {} xsecweight=nan from norms", iParam);
    #endif
  }
  return xsecw;
}

CalcWeightSpline()

M3::float_t SampleHandlerFD::CalcWeightSpline ( const FarDetectorCoreInfo *  MCEvent ) const

protected

Calculate the spline weight for a given event.

Definition at line 737 of file SampleHandlerFD.cpp.

                                                                                    {
// ***************************************************************************
  M3::float_t spline_weight = 1.0;
  const int nSplines = static_cast<int>(MCEvent->xsec_spline_pointers.size());
  //DB Xsec syst
  //Loop over stored spline pointers
  #ifdef MULTITHREAD
  #pragma omp simd
  #endif
  for (int iSpline = 0; iSpline < nSplines; ++iSpline) {
    spline_weight *= *(MCEvent->xsec_spline_pointers[iSpline]);
  }
  return spline_weight;
}

Fill1DHist()

void SampleHandlerFD::Fill1DHist ( )

protected

Fill a 1D histogram with the event-level information used in the fit.

Definition at line 277 of file SampleHandlerFD.cpp.

{
  // DB Commented out by default - Code heading towards GetLikelihood using arrays instead of root objects
  // Wouldn't actually need this for GetLikelihood as TH objects wouldn't be filled
  _hPDF1D->Reset();
  for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
    for (unsigned int xBin=0;xBin<(Binning.XBinEdges.size()-1);xBin++) {
      _hPDF1D->AddBinContent(xBin+1,SampleHandlerFD_array[yBin][xBin]);
    }
  }
}

Fill2DHist()

void SampleHandlerFD::Fill2DHist ( )

protected

Fill a 2D histogram with the event-level information used in the fit.

Definition at line 289 of file SampleHandlerFD.cpp.

{
  // DB Commented out by default - Code heading towards GetLikelihood using arrays instead of root objects
  // Wouldn't actually need this for GetLikelihood as TH objects wouldn't be filled
  _hPDF2D->Reset();
  for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
    for (unsigned int xBin=0;xBin<(Binning.XBinEdges.size()-1);xBin++) {
      _hPDF2D->SetBinContent(xBin+1,yBin+1,SampleHandlerFD_array[yBin][xBin]);
    }
  }
}

Fill2DSubEventHist()

void SampleHandlerFD::Fill2DSubEventHist	(	TH2D *	_h2DVar,
		const std::string &	ProjectionVarX,
		const std::string &	ProjectionVarY,
		const std::vector< KinematicCut > &	SubEventSelectionVec = std::vector< KinematicCut >(),
		int	WeightStyle = 0
	)

Definition at line 1634 of file SampleHandlerFD.cpp.

                                                                            {
  bool IsSubEventVarX = IsSubEventVarString(ProjectionVar_StrX);
  bool IsSubEventVarY = IsSubEventVarString(ProjectionVar_StrY);   
 
  int ProjectionVar_IntX, ProjectionVar_IntY;
  if (IsSubEventVarX) ProjectionVar_IntX = ReturnKinematicVectorFromString(ProjectionVar_StrX);
  else ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
  if (IsSubEventVarY) ProjectionVar_IntY = ReturnKinematicVectorFromString(ProjectionVar_StrY);
  else ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY); 
 
  //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> VecX = {}, VecY = {};
      double VarX = M3::_BAD_DOUBLE_, VarY = M3::_BAD_DOUBLE_;
      size_t nsubevents = 0;
      // JM Three cases: subeventX vs eventY || eventX vs subeventY || subeventX vs subeventY
      if (IsSubEventVarX && !IsSubEventVarY) {
        VecX = ReturnKinematicVector(ProjectionVar_IntX, iEvent);
        VarY = ReturnKinematicParameter(ProjectionVar_IntY, iEvent);
        nsubevents = VecX.size();
      }
      else if (!IsSubEventVarX && IsSubEventVarY) {
        VecY = ReturnKinematicVector(ProjectionVar_IntY, iEvent);
        VarX = ReturnKinematicParameter(ProjectionVar_IntX, iEvent);
        nsubevents = VecY.size();
      }
      else {
        VecX = ReturnKinematicVector(ProjectionVar_IntX, iEvent);
        VecY = ReturnKinematicVector(ProjectionVar_IntY, iEvent);
        if (VecX.size() != VecY.size()) {
          MACH3LOG_ERROR("Cannot plot {} of size {} against {} of size {}", ProjectionVar_StrX, VecX.size(), ProjectionVar_StrY, VecY.size());
          throw MaCh3Exception(__FILE__, __LINE__);
        }
        nsubevents = VecX.size();
      }
      //JM Loop over all subevents in event
      for (unsigned int iSubEvent = 0; iSubEvent < nsubevents; iSubEvent++) {
        if (IsSubEventSelected(SubEventSelectionVec, iEvent, iSubEvent, nsubevents)) {
          if (IsSubEventVarX) VarX = VecX[iSubEvent];
          if (IsSubEventVarY) VarY = VecY[iSubEvent];
          _h2DVar->Fill(VarX,VarY,Weight);
        }
      } 
    }
  }
}

FillArray()

void SampleHandlerFD::FillArray ( )

protected

DB Nice new multi-threaded function which calculates the event weights and fills the relevant bins of an array.

fills the SampleHandlerFD_array vector with the weight calculated from reweighting

@function SampleHandlerFD::fillArray() Function which does the core reweighting. This assumes that oscillation weights have already been calculated and stored in SampleHandlerFD.osc_w[iEvent]. This function takes advantage of most of the things called in setupSKMC to reduce reweighting time. It also follows the ND code reweighting pretty closely. This function fills the SampleHandlerFD array array which is binned to match the sample binning, such that bin[1][1] is the equivalent of _hPDF2D->GetBinContent(2,2) {Noticing the offset}

Definition at line 427 of file SampleHandlerFD.cpp.

                                {
//************************************************
  //DB Reset which cuts to apply
  Selection = StoredSelection;
  
  PrepFunctionalParameters();
 
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    ApplyShifts(iEvent);
 
    if (!IsEventSelected(iEvent)) {
      continue;
    }
 
    FarDetectorCoreInfo* MCEvent = &MCSamples[iEvent];
    M3::float_t splineweight = CalcWeightSpline(MCEvent);
    //DB Catch negative spline weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient. Do this on a spline-by-spline basis
    if (splineweight <= 0.){
      MCEvent->xsec_w = 0.;
      continue;
    }
 
    //Loop over stored normalisation and function pointers
    M3::float_t normweight = CalcWeightNorm(MCEvent);
 
    // Virtual by default does nothing
    CalcWeightFunc(iEvent);
 
    MCEvent->xsec_w = splineweight*normweight;
 
    //DB Total weight
    M3::float_t totalweight = GetEventWeight(iEvent);
    //DB Catch negative weights and skip any event with a negative event
    if (totalweight <= 0.){
      MCEvent->xsec_w = 0.;
      continue;
    }
    //DB Switch on BinningOpt to allow different binning options to be implemented
    //The alternative would be to have inheritance based on BinningOpt
    const double XVar = *(MCEvent->x_var);
 
    //DB Find the relevant bin in the PDF for each event
    const int XBinToFill = Binning.FindXBin(XVar, MCEvent->NomXBin);
    const int YBinToFill = MCEvent->NomYBin;
 
    //DB Fill relevant part of thread array
    if (XBinToFill != -1 && YBinToFill != -1) {
      SampleHandlerFD_array[YBinToFill][XBinToFill] += totalweight;
      if (FirstTimeW2) SampleHandlerFD_array_w2[YBinToFill][XBinToFill] += totalweight*totalweight;
    }
  }
}

FillSplineBins()

void SampleHandlerFD::FillSplineBins ( )

protected

Finds the binned spline that an event should apply to and stored them in a a vector for easy evaluation in the fillArray() function.

@func fillSplineBins()

Definition at line 1371 of file SampleHandlerFD.cpp.

                                     {
  //Now loop over events and get the spline bin for each event
  for (unsigned int j = 0; j < GetNEvents(); ++j) {
    const int OscIndex = GetOscChannel(OscChannels, (*MCSamples[j].nupdgUnosc), (*MCSamples[j].nupdg));
 
    std::vector< std::vector<int> > EventSplines;
    switch(nDimensions){
      case 1:
        EventSplines = SplineHandler->GetEventSplines(GetSampleName(), OscIndex, int(*(MCSamples[j].mode)), *(MCSamples[j].rw_etru), *(MCSamples[j].x_var), 0.);
        break;
      case 2:
        EventSplines = SplineHandler->GetEventSplines(GetSampleName(), OscIndex, int(*(MCSamples[j].mode)), *(MCSamples[j].rw_etru), *(MCSamples[j].x_var), *(MCSamples[j].y_var));
        break;
      default:
        MACH3LOG_ERROR("Error in assigning spline bins because nDimensions = {}", nDimensions);
        MACH3LOG_ERROR("MaCh3 only supports splines binned in Etrue + the sample binning");
        MACH3LOG_ERROR("Please check the sample binning you specified in your sample config ");
        throw MaCh3Exception(__FILE__, __LINE__);
        break;
    }
    int NSplines = int(EventSplines.size());
    if(NSplines < 0){
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    MCSamples[j].xsec_spline_pointers.resize(NSplines);
    for(size_t spline = 0; spline < MCSamples[j].xsec_spline_pointers.size(); spline++) {
      //Event Splines indexed as: sample name, oscillation channel, syst, mode, etrue, var1, var2 (var2 is a dummy 0 for 1D splines)
      MCSamples[j].xsec_spline_pointers[spline] = SplineHandler->retPointer(EventSplines[spline][0], EventSplines[spline][1],
                                                                                EventSplines[spline][2], EventSplines[spline][3],
                                                                                EventSplines[spline][4], EventSplines[spline][5],
                                                                                EventSplines[spline][6]);
    }
  }
}

FindNominalBinAndEdges1D()

void SampleHandlerFD::FindNominalBinAndEdges1D ( )

protected

Definition at line 963 of file SampleHandlerFD.cpp.

                                               {
  for(unsigned int event_i = 0; event_i < GetNEvents(); event_i++){
    //Set x_var and y_var values based on XVarStr and YVarStr
    MCSamples[event_i].x_var = GetPointerToKinematicParameter(XVarStr, event_i);
    if (std::isnan(*MCSamples[event_i].x_var) || std::isinf(*MCSamples[event_i].x_var)) {
      MACH3LOG_ERROR("X var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].x_var);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
 
    //Give y_var M3::_BAD_DOUBLE_ value for the 1D case since this won't be used
    MCSamples[event_i].y_var = &(M3::_BAD_DOUBLE_);
    int bin = _hPDF1D->FindBin(*(MCSamples[event_i].x_var));
 
    if ((bin-1) >= 0 && (bin-1) < int(Binning.XBinEdges.size()-1)) {
      MCSamples[event_i].NomXBin = bin-1;
    } else {
      MCSamples[event_i].NomXBin = -1;
    }
    MCSamples[event_i].NomYBin = 0;
  }
 
  Binning.rw_lower_xbinedge.resize(Binning.nXBins);
  Binning.rw_lower_lower_xbinedge.resize(Binning.nXBins);
  Binning.rw_upper_xbinedge.resize(Binning.nXBins);
  Binning.rw_upper_upper_xbinedge.resize(Binning.nXBins);
  //Set rw_pdf_bin and rw_upper_xbinedge and rw_lower_xbinedge for each skmc_base
  for(size_t bin_x = 0; bin_x < Binning.nXBins; bin_x++){
    double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
    double low_edge = Binning.XBinEdges[bin_x];
    double upper_edge = Binning.XBinEdges[bin_x+1];
    double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
 
    if (bin_x == 0) {
      low_lower_edge = Binning.XBinEdges[0];
    } else {
      low_lower_edge = Binning.XBinEdges[bin_x-1];
    }
 
    if (bin_x + 2 < Binning.nXBins) {
      upper_upper_edge = Binning.XBinEdges[bin_x + 2];
    } else if (bin_x + 1 < Binning.nXBins) {
      upper_upper_edge = Binning.XBinEdges[bin_x + 1];
    }
 
    Binning.rw_lower_xbinedge[bin_x] = low_edge;
    Binning.rw_upper_xbinedge[bin_x] = upper_edge;
    Binning.rw_lower_lower_xbinedge[bin_x] = low_lower_edge;
    Binning.rw_upper_upper_xbinedge[bin_x] = upper_upper_edge;
  }
}

FindNominalBinAndEdges2D()

void SampleHandlerFD::FindNominalBinAndEdges2D ( )

protected

Definition at line 1031 of file SampleHandlerFD.cpp.

                                               {
// ************************************************
  for(unsigned int event_i = 0 ; event_i < GetNEvents(); event_i++) {
    //Set x_var and y_var values based on XVarStr and YVarStr
    MCSamples[event_i].x_var = GetPointerToKinematicParameter(XVarStr, event_i);
    MCSamples[event_i].y_var = GetPointerToKinematicParameter(YVarStr, event_i);
 
    if (std::isnan(*MCSamples[event_i].x_var) || std::isinf(*MCSamples[event_i].x_var)) {
      MACH3LOG_ERROR("X var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].x_var);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    if (std::isnan(*MCSamples[event_i].y_var) || std::isinf(*MCSamples[event_i].y_var)) {
      MACH3LOG_ERROR("Y var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].y_var);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    //Global bin number
    int bin = _hPDF2D->FindBin(*(MCSamples[event_i].x_var), *(MCSamples[event_i].y_var));
 
    int bin_x = M3::_BAD_INT_;
    int bin_y = M3::_BAD_INT_;
    int bin_z = M3::_BAD_INT_;
    _hPDF2D->GetBinXYZ(bin, bin_x, bin_y, bin_z);
 
    if ((bin_x-1) >= 0 && (bin_x-1) < int(Binning.XBinEdges.size()-1)) {
      MCSamples[event_i].NomXBin = bin_x-1;
    }
    MCSamples[event_i].NomYBin = bin_y-1;
    if(MCSamples[event_i].NomYBin < 0){
      MACH3LOG_WARN("Nominal YBin PROBLEM, y-bin is {}", MCSamples[event_i].NomYBin);
    }
  }
 
  Binning.rw_lower_xbinedge.resize(Binning.nXBins);
  Binning.rw_lower_lower_xbinedge.resize(Binning.nXBins);
  Binning.rw_upper_xbinedge.resize(Binning.nXBins);
  Binning.rw_upper_upper_xbinedge.resize(Binning.nXBins);
  //Set rw_pdf_bin and rw_upper_xbinedge and rw_lower_xbinedge for each skmc_base
  for(size_t bin_x = 0; bin_x < Binning.nXBins; bin_x++){
    double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
    double low_edge = Binning.XBinEdges[bin_x];
    double upper_edge = Binning.XBinEdges[bin_x+1];
    double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
 
    if (bin_x == 0) {
      low_lower_edge = Binning.XBinEdges[0];
    } else {
      low_lower_edge = Binning.XBinEdges[bin_x-1];
    }
 
    if (bin_x + 2 < Binning.nXBins) {
      upper_upper_edge = Binning.XBinEdges[bin_x + 2];
    } else if (bin_x + 1 < Binning.nXBins) {
      upper_upper_edge = Binning.XBinEdges[bin_x + 1];
    }
 
    Binning.rw_lower_xbinedge[bin_x] = low_edge;
    Binning.rw_upper_xbinedge[bin_x] = upper_edge;
    Binning.rw_lower_lower_xbinedge[bin_x] = low_lower_edge;
    Binning.rw_upper_upper_xbinedge[bin_x] = upper_upper_edge;
  }
}

Get2DVarHist()

TH2 * SampleHandlerFD::Get2DVarHist	(	const std::string &	ProjectionVarX,
		const std::string &	ProjectionVarY,
		const std::vector< KinematicCut > &	EventSelectionVec = std::vector< KinematicCut >(),
		int	WeightStyle = 0,
		TAxis *	AxisX = nullptr,
		TAxis *	AxisY = nullptr,
		const std::vector< KinematicCut > &	SubEventSelectionVec = std::vector< KinematicCut >()
	)

Definition at line 1577 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
  TH2D* _h2DVar;
  if (AxisX && AxisY) {
    _h2DVar = new TH2D("","",AxisX->GetNbins(),AxisX->GetXbins()->GetArray(),AxisY->GetNbins(),AxisY->GetXbins()->GetArray());
  } else {
    std::vector<double> xBinEdges = ReturnKinematicParameterBinning(ProjectionVar_StrX);
    std::vector<double> yBinEdges = ReturnKinematicParameterBinning(ProjectionVar_StrY);
    _h2DVar = new TH2D("", "", int(xBinEdges.size())-1, xBinEdges.data(), int(yBinEdges.size())-1, yBinEdges.data());
  }
 
  bool IsSubEventHist = IsSubEventVarString(ProjectionVar_StrX) || IsSubEventVarString(ProjectionVar_StrY);
  if (IsSubEventHist) Fill2DSubEventHist(_h2DVar, ProjectionVar_StrX, ProjectionVar_StrY, SubEventSelectionVec, WeightStyle);
  else {
    //DB Grab the associated enum with the argument string
    int ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
    int ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY);
 
    //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 VarX = ReturnKinematicParameter(ProjectionVar_IntX, iEvent);
        double VarY = ReturnKinematicParameter(ProjectionVar_IntY, iEvent);
        _h2DVar->Fill(VarX,VarY,Weight);
      }
    }
  }
  //DB Reset the saved selection
  Selection = tmp_Selection;
 
  return _h2DVar;
}

GetEventWeight()

M3::float_t SampleHandlerFD::GetEventWeight

(

const int

iEntry

)

const

Definition at line 1355 of file SampleHandlerFD.cpp.

                                                                {
  M3::float_t totalweight = 1.0;
  const int nParams = static_cast<int>(MCSamples[iEntry].total_weight_pointers.size());
  #ifdef MULTITHREAD
  #pragma omp simd
  #endif
  for (int iParam = 0; iParam < nParams; ++iParam) {
    totalweight *= *(MCSamples[iEntry].total_weight_pointers[iParam]);
  }
 
  return totalweight;
}

GetFinalPDGFromFileName()

NuPDG SampleHandlerFD::GetFinalPDGFromFileName ( const std::string &  FileName ) const

inlineprotected

Retrieve the final neutrino PDG code associated with a given input file name.

Definition at line 392 of file SampleHandlerFD.h.

{return FileToFinalPDGMap.at(FileName);}

GetInitPDGFromFileName()

NuPDG SampleHandlerFD::GetInitPDGFromFileName ( const std::string &  FileName ) const

inlineprotected

Retrieve the initial neutrino PDG code associated with a given input file name.

Definition at line 390 of file SampleHandlerFD.h.

{return FileToInitPDGMap.at(FileName);}

GetLikelihood()

double SampleHandlerFD::GetLikelihood ( )

overridevirtual

DB Multi-threaded GetLikelihood.

Implements SampleHandlerBase.

Definition at line 1407 of file SampleHandlerFD.cpp.

                                      {
// ************************************************
  if (SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("Data sample is empty! Can't calculate a likelihood!");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
    
  double negLogL = 0.;
  #ifdef MULTITHREAD
  #pragma omp parallel for collapse(2) reduction(+:negLogL)
  #endif
  for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin)
  {
    for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin)
    {
      const double DataVal = SampleHandlerFD_data[yBin][xBin];
      const double MCPred = SampleHandlerFD_array[yBin][xBin];
      const double w2 = SampleHandlerFD_array_w2[yBin][xBin];
      
      //KS: Calculate likelihood using Barlow-Beeston Poisson or even IceCube
      negLogL += GetTestStatLLH(DataVal, MCPred, w2);
    }
  }
  return negLogL;
}

GetPointerToKinematicParameter() [1/2]

virtual const double * SampleHandlerFD::GetPointerToKinematicParameter ( double  KinematicVariable, int  iEvent  )

protectedpure virtual

Implemented in PySampleHandlerFD.

GetPointerToKinematicParameter() [2/2]

virtual const double * SampleHandlerFD::GetPointerToKinematicParameter ( std::string  KinematicParamter, int  iEvent  )

protectedpure virtual

Implemented in PySampleHandlerFD.

GetPointerToOscChannel()

const double * SampleHandlerFD::GetPointerToOscChannel ( const int  iEvent ) const

protected

Get pointer to oscillation channel associated with given event. Osc channel is const.

Definition at line 2106 of file SampleHandlerFD.cpp.

                                                                      {
// ************************************************
  const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
 
  return &(OscChannels[Channel].ChannelIndex);
}

Init()

virtual void SampleHandlerFD::Init ( )

protectedpure virtual

Initialise any variables that your experiment specific SampleHandler needs.

Implemented in PySampleHandlerFD.

Initialise()

void SampleHandlerFD::Initialise ( )

protected

Function which does a lot of the lifting regarding the workflow in creating different MC objects.

Definition at line 202 of file SampleHandlerFD.cpp.

                                 {
  //First grab all the information from your sample config via your manager
  ReadSampleConfig();
 
  //Now initialise all the variables you will need
  Init();
 
  nEvents = SetupExperimentMC();
 
  InitialiseSingleFDMCObject();
  SetupFDMC();
 
  MACH3LOG_INFO("=============================================");
  MACH3LOG_INFO("Total number of events is: {}", GetNEvents());
 
  auto OscParams = ParHandler->GetOscParsFromSampleName(SampleName);
  if (OscParams.size() > 0) {
    MACH3LOG_INFO("Setting up NuOscillator..");
    if (Oscillator != nullptr) {
      MACH3LOG_INFO("You have passed an OscillatorBase object through the constructor of a SampleHandlerFD object - this will be used for all oscillation channels");
      if(Oscillator->isEqualBinningPerOscChannel() != true) {
        MACH3LOG_ERROR("Trying to run shared NuOscillator without EqualBinningPerOscChannel, this will not work");
        throw MaCh3Exception(__FILE__, __LINE__);
      }
 
      if(OscParams.size() != Oscillator->GetOscParamsSize()){
        MACH3LOG_ERROR("Sample {} with {} has {} osc params, while shared NuOsc has {} osc params", GetTitle(), GetSampleName(),
                       OscParams.size(), Oscillator->GetOscParamsSize());
        MACH3LOG_ERROR("This indicate misconfiguration in your Osc yaml");
        throw MaCh3Exception(__FILE__, __LINE__);
      }
    } else {
      InitialiseNuOscillatorObjects();
    }
    SetupNuOscillatorPointers();
  } else{
    MACH3LOG_WARN("Didn't find any oscillation params, thus will not enable oscillations");
  }
 
  MACH3LOG_INFO("Setting up Sample Binning..");
  SetupSampleBinning();
  MACH3LOG_INFO("Setting up Splines..");
  SetupSplines();
  MACH3LOG_INFO("Setting up Normalisation Pointers..");
  SetupNormParameters();
  MACH3LOG_INFO("Setting up Functional Pointers..");
  SetupFunctionalParameters();
  MACH3LOG_INFO("Setting up Weight Pointers..");
  SetupWeightPointers();
  MACH3LOG_INFO("Setting up Kinematic Map..");
  SetupKinematicMap();
  MACH3LOG_INFO("=======================================================");
}

InitialiseNuOscillatorObjects()

void SampleHandlerFD::InitialiseNuOscillatorObjects

(

)

including Dan's magic NuOscillator

Definition at line 1250 of file SampleHandlerFD.cpp.

                                                    {
// ************************************************
  auto NuOscillatorConfigFile = Get<std::string>(SampleManager->raw()["NuOsc"]["NuOscConfigFile"], __FILE__ , __LINE__);
  auto EqualBinningPerOscChannel = Get<bool>(SampleManager->raw()["NuOsc"]["EqualBinningPerOscChannel"], __FILE__ , __LINE__);
 
  // TN override the sample setting if not using binned oscillation
  if (EqualBinningPerOscChannel) {
    if (YAML::LoadFile(NuOscillatorConfigFile)["General"]["CalculationType"].as<std::string>() == "Unbinned") {
      MACH3LOG_WARN("Tried using EqualBinningPerOscChannel while using Unbinned oscillation calculation, changing EqualBinningPerOscChannel to false");
      EqualBinningPerOscChannel = false;
    }
  }
  std::vector<const double*> OscParams = ParHandler->GetOscParsFromSampleName(SampleName);
  if (OscParams.empty()) {
    MACH3LOG_ERROR("OscParams is empty for sample '{}'.", GetTitle());
    MACH3LOG_ERROR("This likely indicates an error in your oscillation YAML configuration.");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  Oscillator = std::make_shared<OscillationHandler>(NuOscillatorConfigFile, EqualBinningPerOscChannel, OscParams, static_cast<int>(OscChannels.size()));
 
  if (!EqualBinningPerOscChannel) {
    for(int iChannel = 0; iChannel < static_cast<int>(OscChannels.size()); iChannel++) {
      std::vector<M3::float_t> EnergyArray;
      std::vector<M3::float_t> CosineZArray;
 
      for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
        // KS: This is bit weird but we basically loop over all events and push to vector only these which are part of a given OscChannel
        const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
        //DB Remove NC events from the arrays which are handed to the NuOscillator objects
        if (!MCSamples[iEvent].isNC && Channel == iChannel) {
          EnergyArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_etru)));
        }
      }
      std::sort(EnergyArray.begin(),EnergyArray.end());
 
      //============================================================================
      //DB Atmospheric only part, can only happen if truecz has been initialised within the experiment specific code
      if (*(MCSamples[0].rw_truecz) != M3::_BAD_DOUBLE_) {
        for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
          // KS: This is bit weird but we basically loop over all events and push to vector only these which are part of a given OscChannel
          const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
          //DB Remove NC events from the arrays which are handed to the NuOscillator objects
          if (!MCSamples[iEvent].isNC && Channel == iChannel) {
            CosineZArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_truecz)));
          }
        }
        std::sort(CosineZArray.begin(),CosineZArray.end());
      }
      Oscillator->SetOscillatorBinning(iChannel, EnergyArray, CosineZArray);
    } // end loop over channels
  }
}

InitialiseSingleFDMCObject()

void SampleHandlerFD::InitialiseSingleFDMCObject ( )

protected

function to create the member of the FarDetectorInfo struct so they are the appropriate size.

Definition at line 1433 of file SampleHandlerFD.cpp.

                                                 {
  MCSamples.resize(nEvents);
}

InitialiseSplineObject()

void SampleHandlerFD::InitialiseSplineObject ( )

protected

Definition at line 1472 of file SampleHandlerFD.cpp.

                                             {
  std::vector<std::string> spline_filepaths;
  for(unsigned iChannel = 0 ; iChannel < OscChannels.size() ; iChannel++){
    spline_filepaths.push_back(spline_files[iChannel]);
  }
 
  //Keep a track of the spline variables
  std::vector<std::string> SplineVarNames = {"TrueNeutrinoEnergy"};
  if(XVarStr.length() > 0){
    SplineVarNames.push_back(XVarStr);
  }
  if(YVarStr.length() > 0){
    SplineVarNames.push_back(YVarStr);
  }
  
  SplineHandler->AddSample(SampleName, spline_filepaths, SplineVarNames);
  SplineHandler->CountNumberOfLoadedSplines(false, 1);
  SplineHandler->TransferToMonolith();
 
  MACH3LOG_INFO("--------------------------------");
  MACH3LOG_INFO("Setup Far Detector splines");
 
  FillSplineBins();
 
  SplineHandler->cleanUpMemory();
}

IsEventSelected()

bool SampleHandlerFD::IsEventSelected ( const int  iEvent )

protected

DB Function which determines if an event is selected, where Selection double looks like {{ND280KinematicTypes Var1, douuble LowBound}.

Definition at line 362 of file SampleHandlerFD.cpp.

                                                      {
// ************************************************
  const int SelectionSize = static_cast<int>(Selection.size());
  for (int iSelection = 0; iSelection < SelectionSize; ++iSelection) {
    const auto& Cut = Selection[iSelection];
    const double Val = ReturnKinematicParameter(Cut.ParamToCutOnIt, iEvent);
    if ((Val < Cut.LowerBound) || (Val >= Cut.UpperBound)) {
      return false;
    }
  }
  //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check
  return true;
}

IsSubEventSelected()

bool SampleHandlerFD::IsSubEventSelected ( const std::vector< KinematicCut > &  SubEventCuts, const int  iEvent, unsigned const int  iSubEvent, size_t  nsubevents  )

protected

JM Function which determines if a subevent is selected.

Definition at line 377 of file SampleHandlerFD.cpp.

                                                                                                                                                       {
  for (unsigned int iSelection=0;iSelection < SubEventCuts.size() ;iSelection++) {
    std::vector<double> Vec = ReturnKinematicVector(SubEventCuts[iSelection].ParamToCutOnIt, iEvent);
    if (nsubevents != Vec.size()) {
      MACH3LOG_ERROR("Cannot apply kinematic cut on {} as it is of different size to plotting variable");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    const double Val = Vec[iSubEvent];
    if ((Val < SubEventCuts[iSelection].LowerBound) || (Val >= SubEventCuts[iSelection].UpperBound)) {
      return false;
    }
  }
  //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check
  return true;
}

IsSubEventVarString()

bool SampleHandlerFD::IsSubEventVarString

(

const std::string &

VarStr

)

JM Check if a kinematic parameter string corresponds to a subevent-level variable.

Definition at line 1776 of file SampleHandlerFD.cpp.

                                                                 {
  if (KinematicVectors == nullptr) return false;
 
  if (KinematicVectors->count(VarStr)) {
    if (!KinematicParameters->count(VarStr)) return true;
    else {
      MACH3LOG_ERROR("Attempted to plot kinematic variable {}, but it appears in both KinematicVectors and KinematicParameters", VarStr);
      throw MaCh3Exception(__FILE__,__LINE__);
    }
  }
  return false;
}

PrepFunctionalParameters()

virtual void SampleHandlerFD::PrepFunctionalParameters ( )

inlineprotectedvirtual

Update the functional parameter values to the latest proposed values. Needs to be called before every new reweight so is called in fillArray.

Definition at line 219 of file SampleHandlerFD.h.

{};

PrintIntegral()

void SampleHandlerFD::PrintIntegral	(	TString	OutputName = "/dev/null",
		int	WeightStyle = 0,
		TString	OutputCSVName = "/dev/null"
	)

Definition at line 1890 of file SampleHandlerFD.cpp.

                                                                                                      {
  int space = 14;
 
  bool printToFile=false;
  if (OutputFileName.CompareTo("/dev/null")) {printToFile = true;}
 
  bool printToCSV=false;
  if(OutputCSVFileName.CompareTo("/dev/null")) printToCSV=true;
 
  std::ofstream outfile;
  if (printToFile) {
    outfile.open(OutputFileName.Data(), std::ios_base::app);
    outfile.precision(7);
  }
 
  std::ofstream outcsv;
  if(printToCSV){
    outcsv.open(OutputCSVFileName, std::ios_base::app); // Appened to CSV
    outcsv.precision(7);
  }
 
  double PDFIntegral = 0.;
 
  std::vector< std::vector< TH1* > > IntegralList;
  IntegralList.resize(Modes->GetNModes());
 
  std::vector<double> ChannelIntegral;
  ChannelIntegral.resize(GetNOscChannels());
  for (unsigned int i=0;i<ChannelIntegral.size();i++) {ChannelIntegral[i] = 0.;}
 
  for (int i=0;i<Modes->GetNModes();i++) {
    if (GetNDim()==1) {
      IntegralList[i] = ReturnHistsBySelection1D(XVarStr,1,i,WeightStyle);
    } else {
      IntegralList[i] = CastVector<TH2, TH1>(ReturnHistsBySelection2D(XVarStr,YVarStr,1,i,WeightStyle));
    }
  }
 
  MACH3LOG_INFO("-------------------------------------------------");
 
  if (printToFile) {
    outfile << "\\begin{table}[ht]" << std::endl;
    outfile << "\\begin{center}" << std::endl;
    outfile << "\\caption{Integral breakdown for sample: " << GetTitle() << "}" << std::endl;
    outfile << "\\label{" << GetTitle() << "-EventRate}" << std::endl;
 
    TString nColumns;
    for (int i=0;i<GetNOscChannels();i++) {nColumns+="|c";}
    nColumns += "|c|";
    outfile << "\\begin{tabular}{|l" << nColumns.Data() << "}" << std::endl;
    outfile << "\\hline" << std::endl;
  }
 
  if(printToCSV){
    // HW Probably a better way but oh well, here I go making MaCh3 messy again
    outcsv<<"Integral Breakdown for sample :"<<GetTitle()<<"\n";
  }
 
  MACH3LOG_INFO("Integral breakdown for sample: {}", GetTitle());
  MACH3LOG_INFO("");
 
  if (printToFile) {outfile << std::setw(space) << "Mode:";}
  if(printToCSV) {outcsv<<"Mode,";}
 
  std::string table_headings = fmt::format("| {:<8} |", "Mode");
  std::string table_footline = "------------"; //Scalable table horizontal line
  for (size_t i = 0;i < OscChannels.size(); i++) {
    table_headings += fmt::format(" {:<17} |", OscChannels[i].flavourName);
    table_footline += "--------------------";
    if (printToFile) {outfile << "&" << std::setw(space) << OscChannels[i].flavourName_Latex << " ";}
    if (printToCSV)  {outcsv << OscChannels[i].flavourName << ",";}
  }
  if (printToFile) {outfile << "&" << std::setw(space) << "Total:" << "\\\\ \\hline" << std::endl;}
  if (printToCSV)  {outcsv <<"Total\n";}
  table_headings += fmt::format(" {:<10} |", "Total");
  table_footline += "-------------";
 
  MACH3LOG_INFO("{}", table_headings);
  MACH3LOG_INFO("{}", table_footline);
 
  for (unsigned int i=0;i<IntegralList.size();i++) {
    double ModeIntegral = 0;
    if (printToFile) {outfile << std::setw(space) << Modes->GetMaCh3ModeName(i);}
    if(printToCSV)   {outcsv << Modes->GetMaCh3ModeName(i) << ",";}
 
    table_headings = fmt::format("| {:<8} |", Modes->GetMaCh3ModeName(i)); //Start string with mode name
 
    for (unsigned int j=0;j<IntegralList[i].size();j++) {
      double Integral = IntegralList[i][j]->Integral();
 
      if (Integral<1e-100) {Integral=0;}
 
      ModeIntegral += Integral;
      ChannelIntegral[j] += Integral;
      PDFIntegral += Integral;
 
      if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",Integral) << " ";}
      if (printToCSV)  {outcsv << Form("%4.5f", Integral) << ",";}
 
      table_headings += fmt::format(" {:<17.4f} |", Integral);
    }
    if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ModeIntegral) <<  " \\\\ \\hline" << std::endl;}
    if (printToCSV)  {outcsv << Form("%4.5f", ModeIntegral) << "\n";}
 
    table_headings += fmt::format(" {:<10.4f} |", ModeIntegral);
 
    MACH3LOG_INFO("{}", table_headings);
  }
 
  if (printToFile) {outfile << std::setw(space) << "Total:";}
  if (printToCSV)  {outcsv << "Total,";}
 
  //Clear the table_headings to print last row of totals
  table_headings = fmt::format("| {:<8} |", "Total");
  for (unsigned int i=0;i<ChannelIntegral.size();i++) {
    if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ChannelIntegral[i]) << " ";}
    if (printToCSV)  {outcsv << Form("%4.5f", ChannelIntegral[i]) << ",";}
    table_headings += fmt::format(" {:<17.4f} |", ChannelIntegral[i]);
  }
  if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",PDFIntegral) << " \\\\ \\hline" << std::endl;}
  if (printToCSV)  {outcsv << Form("%4.5f", PDFIntegral) << "\n\n\n\n";} // Let's have a few new lines!
 
  table_headings += fmt::format(" {:<10.4f} |", PDFIntegral);
  MACH3LOG_INFO("{}", table_headings);
  MACH3LOG_INFO("{}", table_footline);
 
  if (printToFile) {
    outfile << "\\end{tabular}" << std::endl;
    outfile << "\\end{center}" << std::endl;
    outfile << "\\end{table}" << std::endl;
  }
 
  MACH3LOG_INFO("");
 
  if (printToFile) {
    outfile << std::endl;
    outfile.close();
  }
}

ReadSampleConfig()

void SampleHandlerFD::ReadSampleConfig

(

)

Definition at line 73 of file SampleHandlerFD.cpp.

{
  auto ModeName = Get<std::string>(SampleManager->raw()["MaCh3ModeConfig"], __FILE__ , __LINE__);
  Modes = std::make_unique<MaCh3Modes>(ModeName);
  //SampleTitle has to be provided in the sample yaml otherwise this will throw an exception
  SampleTitle = Get<std::string>(SampleManager->raw()["SampleTitle"], __FILE__ , __LINE__);
  //SampleName has to be provided in the sample yaml otherwise this will throw an exception
  SampleName = Get<std::string>(SampleManager->raw()["SampleName"], __FILE__ , __LINE__);
 
  fTestStatistic = static_cast<TestStatistic>(SampleManager->GetMCStatLLH());
  if (CheckNodeExists(SampleManager->raw(), "LikelihoodOptions")) {
    UpdateW2 = GetFromManager<bool>(SampleManager->raw()["LikelihoodOptions"]["UpdateW2"], false);
  }
  //Binning
  nDimensions = 0;
  XVarStr = GetFromManager(SampleManager->raw()["Binning"]["XVarStr"], std::string(""));
  Binning.XBinEdges = GetFromManager(SampleManager->raw()["Binning"]["XVarBins"], std::vector<double>());
  const auto& edgesx = Binning.XBinEdges;
  if (!std::is_sorted(edgesx.begin(), edgesx.end())) {
    MACH3LOG_ERROR("XVarBins must be in increasing order in sample config {}\n  XVarBins: [{}]",
                   GetTitle(), fmt::join(edgesx, ", "));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  if(XVarStr.length() > 0){
    nDimensions++;
  } else{
    MACH3LOG_ERROR("Please specify an X-variable string in sample config {}", SampleManager->GetFileName());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  
  YVarStr = GetFromManager(SampleManager->raw()["Binning"]["YVarStr"], std::string(""));
  Binning.YBinEdges = GetFromManager(SampleManager->raw()["Binning"]["YVarBins"], std::vector<double>());
  const auto& edgesy = Binning.YBinEdges;
  if (!std::is_sorted(edgesy.begin(), edgesy.end())) {
    MACH3LOG_ERROR("YBinEdges must be in increasing order in sample config {}\n  YBinEdges: [{}]",
                   GetTitle(), fmt::join(edgesy, ", "));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  if(YVarStr.length() > 0){
    if(XVarStr.length() == 0){
      MACH3LOG_ERROR("Please specify an X-variable string in sample config {}. I won't work only with a Y-variable", SampleManager->GetFileName());
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    nDimensions++;
  }
  
  if(nDimensions == 0){
    MACH3LOG_ERROR("Error setting up the sample binning");
    MACH3LOG_ERROR("Number of dimensions is {}", nDimensions);
    MACH3LOG_ERROR("Check that an XVarStr has been given in the sample config");
    throw MaCh3Exception(__FILE__, __LINE__);
  } else{
    MACH3LOG_INFO("Found {} dimensions for sample binning", nDimensions);
  }
  
  //Sanity check that some binning has been specified
  if(Binning.XBinEdges.size() == 0 && Binning.YBinEdges.size() == 0){
    MACH3LOG_ERROR("No binning specified for either X or Y of sample binning, please add some binning to the sample config {}", SampleManager->GetFileName());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  
  if (!CheckNodeExists(SampleManager->raw(), "BinningFile")){
    MACH3LOG_ERROR("BinningFile not given in for sample {}, ReturnKinematicParameterBinning will not work", GetTitle());
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  auto mtupleprefix  = Get<std::string>(SampleManager->raw()["InputFiles"]["mtupleprefix"], __FILE__, __LINE__);
  auto mtuplesuffix  = Get<std::string>(SampleManager->raw()["InputFiles"]["mtuplesuffix"], __FILE__, __LINE__);
  auto splineprefix  = Get<std::string>(SampleManager->raw()["InputFiles"]["splineprefix"], __FILE__, __LINE__);
  auto splinesuffix  = Get<std::string>(SampleManager->raw()["InputFiles"]["splinesuffix"], __FILE__, __LINE__);
  
  int NChannels = static_cast<M3::int_t>(SampleManager->raw()["SubSamples"].size());
  OscChannels.reserve(NChannels);
 
  for (auto const &osc_channel : SampleManager->raw()["SubSamples"]) {
    std::string MTupleFileName = mtupleprefix+osc_channel["mtuplefile"].as<std::string>()+mtuplesuffix;
    
    OscChannelInfo OscInfo;
    OscInfo.flavourName       = osc_channel["Name"].as<std::string>();
    OscInfo.flavourName_Latex = osc_channel["LatexName"].as<std::string>();
    OscInfo.InitPDG           = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
    OscInfo.FinalPDG          = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
    OscInfo.ChannelIndex      = static_cast<int>(OscChannels.size());
 
 
    OscChannels.push_back(std::move(OscInfo));
 
    FileToInitPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
    FileToFinalPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
 
    mc_files.push_back(MTupleFileName);
    spline_files.push_back(splineprefix+osc_channel["splinefile"].as<std::string>()+splinesuffix);
  }
 
  //Now grab the selection cuts from the manager
  for ( auto const &SelectionCuts : SampleManager->raw()["SelectionCuts"]) {
    auto TempBoundsVec = GetBounds(SelectionCuts["Bounds"]);
    KinematicCut CutObj;
    CutObj.LowerBound = TempBoundsVec[0];
    CutObj.UpperBound = TempBoundsVec[1];
    CutObj.ParamToCutOnIt = ReturnKinematicParameterFromString(SelectionCuts["KinematicStr"].as<std::string>());
    MACH3LOG_INFO("Adding cut on {} with bounds {} to {}", SelectionCuts["KinematicStr"].as<std::string>(), TempBoundsVec[0], TempBoundsVec[1]);
    StoredSelection.push_back(CutObj);
  }
 
  // EM: initialise the mode weight map
  for( int iMode=0; iMode < Modes->GetNModes(); iMode++ ) {
    _modeNomWeightMap[Modes->GetMaCh3ModeName(iMode)] = 1.0;
  }
 
  // EM: multiply by the nominal weight specified in the sample config file
  if ( SampleManager->raw()["NominalWeights"] ) {
    for( int iMode=0; iMode<Modes->GetNModes(); iMode++ ) {
      std::string modeStr = Modes->GetMaCh3ModeName(iMode);
      if( SampleManager->raw()["NominalWeights"][modeStr] ) {
        double modeWeight = SampleManager->raw()["NominalWeights"][modeStr].as<double>();
        _modeNomWeightMap[Modes->GetMaCh3ModeName(iMode)] *= modeWeight;
      }
    }
  }
 
  // EM: print em out
  MACH3LOG_INFO("  Nominal mode weights to apply: ");
  for(int iMode=0; iMode<Modes->GetNModes(); iMode++ ) {
    std::string modeStr = Modes->GetMaCh3ModeName(iMode);
    MACH3LOG_INFO("    - {}: {}", modeStr, _modeNomWeightMap.at(modeStr));
  }
}

RegisterFunctionalParameters()

virtual void SampleHandlerFD::RegisterFunctionalParameters ( )

protectedpure virtual

HH - a experiment-specific function where the maps to actual functions are set up.

Implemented in PySampleHandlerFD.

RegisterIndividualFunctionalParameter()

void SampleHandlerFD::RegisterIndividualFunctionalParameter ( const std::string &  fpName, int  fpEnum, FuncParFuncType  fpFunc  )

protected

HH - a helper function for RegisterFunctionalParameter.

Definition at line 620 of file SampleHandlerFD.cpp.

                                                                                                                      {
  // Add protections to not add the same functional parameter twice
  if (funcParsNamesMap.find(fpName) != funcParsNamesMap.end()) {
    MACH3LOG_ERROR("Functional parameter {} already registered in funcParsNamesMap with enum {}", fpName, funcParsNamesMap[fpName]);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  if (std::find(funcParsNamesVec.begin(), funcParsNamesVec.end(), fpName) != funcParsNamesVec.end()) {
    MACH3LOG_ERROR("Functional parameter {} already in funcParsNamesVec", fpName);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  if (funcParsFuncMap.find(fpEnum) != funcParsFuncMap.end()) {
    MACH3LOG_ERROR("Functional parameter enum {} already registered in funcParsFuncMap", fpEnum);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  funcParsNamesMap[fpName] = fpEnum;
  funcParsNamesVec.push_back(fpName);
  funcParsFuncMap[fpEnum] = fpFunc;
}

ResetHistograms()

void SampleHandlerFD::ResetHistograms ( )

inlineprotected

Helper function to reset histograms.

Definition at line 605 of file SampleHandlerFD.cpp.

                                      {
// **************************************************  
  //DB Reset values stored in PDF array to 0.
  // Don't openMP this; no significant gain
  for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
    #ifdef MULTITHREAD
    #pragma omp simd
    #endif
    for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
      SampleHandlerFD_array[yBin][xBin] = 0.;
      if(FirstTimeW2) SampleHandlerFD_array_w2[yBin][xBin] = 0.;
    }
  }
} // end function

resetShifts()

virtual void SampleHandlerFD::resetShifts ( int  iEvent )

inlineprotectedvirtual

HH - reset the shifted values to the original values.

Definition at line 228 of file SampleHandlerFD.h.

{(void)iEvent;};

ReturnKinematicParameter() [1/2]

virtual double SampleHandlerFD::ReturnKinematicParameter ( int  KinematicVariable, int  iEvent  )

protectedpure virtual

Implemented in PySampleHandlerFD.

ReturnKinematicParameter() [2/2]

virtual double SampleHandlerFD::ReturnKinematicParameter ( std::string  KinematicParamter, int  iEvent  )

protectedpure virtual

Return the value of an assocaited kinematic parameter for an event.

Implemented in PySampleHandlerFD.

ReturnKinematicParameterBinning()

std::vector< double > SampleHandlerFD::ReturnKinematicParameterBinning ( const std::string &  KinematicParameter )

protected

Return the binning used to draw a kinematic parameter.

Definition at line 1742 of file SampleHandlerFD.cpp.

                                                                                                      {
// ************************************************
  // If x or y variable return used binning
  if(KinematicParameter == XVarStr) {
    return Binning.XBinEdges;
  } else if (KinematicParameter == YVarStr) {
    return Binning.YBinEdges;
  }
 
  auto MakeBins = [](int nBins) {
    std::vector<double> bins(nBins + 1);
    for (int i = 0; i <= nBins; ++i)
      bins[i] = static_cast<double>(i) - 0.5;
    return bins;
  };
 
  if (KinematicParameter == "OscillationChannel") {
    return MakeBins(GetNOscChannels());
  } else if (KinematicParameter == "Mode") {
    return MakeBins(Modes->GetNModes());
  }
 
  // We first check if binning for a sample has been specified
  auto BinningConfig = M3OpenConfig(SampleManager->raw()["BinningFile"].as<std::string>());
  if(BinningConfig[GetTitle()] && BinningConfig[GetTitle()][KinematicParameter]){
    auto BinningVect = Get<std::vector<double>>(BinningConfig[GetTitle()][KinematicParameter], __FILE__, __LINE__);
    return BinningVect;
  } else {
    auto BinningVect = Get<std::vector<double>>(BinningConfig[KinematicParameter], __FILE__, __LINE__);
    return BinningVect;
  }
}

ReturnKinematicVector() [1/2]

virtual std::vector< double > SampleHandlerFD::ReturnKinematicVector ( int  KinematicVariable, int  iEvent  )

inlineprotectedvirtual

Definition at line 267 of file SampleHandlerFD.h.

{return {}; (void)KinematicVariable; (void)iEvent;};

ReturnKinematicVector() [2/2]

virtual std::vector< double > SampleHandlerFD::ReturnKinematicVector ( std::string  KinematicParameter, int  iEvent  )

inlineprotectedvirtual

Definition at line 266 of file SampleHandlerFD.h.

{return {}; (void)KinematicParameter; (void)iEvent;};

ReturnKinematicVectorFromString()

int SampleHandlerFD::ReturnKinematicVectorFromString

(

const std::string &

KinematicStr

)

const

Definition at line 1716 of file SampleHandlerFD.cpp.

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

ReturnStringFromKinematicVector()

std::string SampleHandlerFD::ReturnStringFromKinematicVector

(

const int

KinematicVariable

)

const

Definition at line 1728 of file SampleHandlerFD.cpp.

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

Reweight()

void SampleHandlerFD::Reweight ( )

overridevirtual

Implements SampleHandlerBase.

Definition at line 396 of file SampleHandlerFD.cpp.

                               {
//************************************************
  //KS: Reset the histograms before reweight 
  ResetHistograms();
  
  //You only need to do these things if Oscillator has been initialised
  //if not then you're not considering oscillations
  if (Oscillator) Oscillator->Evaluate();
 
  // Calculate weight coming from all splines if we initialised handler
  if(SplineHandler) SplineHandler->Evaluate();
 
  #ifdef MULTITHREAD
  // Call entirely different routine if we're running with openMP
  FillArray_MP();
  #else
  FillArray();
  #endif
 
  //KS: If you want to not update W2 wights then uncomment this line
  if(!UpdateW2) FirstTimeW2 = false;
}

SaveAdditionalInfo()

void SampleHandlerFD::SaveAdditionalInfo ( TDirectory *  Dir )

overridevirtual

Store additional info in a chan.

Reimplemented from SampleHandlerBase.

Definition at line 1439 of file SampleHandlerFD.cpp.

                                                        {
// ************************************************
  Dir->cd();
 
  YAML::Node Config = SampleManager->raw();
  TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + GetTitle()));
  ConfigSave.Write();
 
  std::unique_ptr<TH1> data_hist;
 
  if (GetNDim() == 1) {
    data_hist = M3::Clone<TH1D>(dynamic_cast<TH1D*>(GetDataHist(1)), "data_" + GetTitle());
    data_hist->GetXaxis()->SetTitle(XVarStr.c_str());
    data_hist->GetYaxis()->SetTitle("Number of Events");
  } else if (GetNDim() == 2) {
    data_hist = M3::Clone<TH2D>(dynamic_cast<TH2D*>(GetDataHist(2)), "data_" + GetTitle());
    data_hist->GetXaxis()->SetTitle(XVarStr.c_str());
    data_hist->GetYaxis()->SetTitle(YVarStr.c_str());
    data_hist->GetZaxis()->SetTitle("Number of Events");
  } else {
    MACH3LOG_ERROR("Not implemented");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  if (!data_hist) {
    MACH3LOG_ERROR("nullptr data hist :(");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  data_hist->SetTitle(("data_" + GetTitle()).c_str());
  data_hist->Write();
}

Set1DBinning() [1/2]

void SampleHandlerFD::Set1DBinning ( size_t  nbins, double *  boundaries  )

protected

sets the binning used for the likelihood calculation, used for both data and MC

Parameters

nbins	number of total bins
boundaries	the bin edges e.g. 0, 0.1, 0.2, 0.3

Definition at line 939 of file SampleHandlerFD.cpp.

{
  _hPDF1D->Reset();
  _hPDF1D->SetBins(static_cast<int>(nbins),boundaries);
  dathist->SetBins(static_cast<int>(nbins),boundaries);
 
  Binning.YBinEdges = std::vector<double>(2);
  Binning.YBinEdges[0] = -1e8;
  Binning.YBinEdges[1] = 1e8;
 
  double YBinEdges_Arr[2];
  YBinEdges_Arr[0] = Binning.YBinEdges[0];
  YBinEdges_Arr[1] = Binning.YBinEdges[1];
 
  _hPDF2D->Reset();
  _hPDF2D->SetBins(static_cast<int>(nbins),boundaries,1,YBinEdges_Arr);
  dathist2d->SetBins(static_cast<int>(nbins),boundaries,1,YBinEdges_Arr);
 
  //Set the number of X and Y bins now
  SetupReweightArrays(Binning.XBinEdges.size() - 1, Binning.YBinEdges.size() - 1);
 
  FindNominalBinAndEdges1D();
}

Set1DBinning() [2/2]

void SampleHandlerFD::Set1DBinning ( std::vector< double > &  XVec )

inlineprotected

set the binning for 1D sample used for the likelihood calculation

Parameters

XVec	vector containing the binning in axis 1 (the x-axis)

Definition at line 197 of file SampleHandlerFD.h.

{Set1DBinning(XVec.size()-1, XVec.data());};

Set2DBinning() [1/2]

void SampleHandlerFD::Set2DBinning ( size_t  nbins1, double *  boundaries1, size_t  nbins2, double *  boundaries2  )

protected

set the binning for 2D sample used for the likelihood calculation

Parameters

nbins1	number of bins in axis 1 (the x-axis)
nbins2	number of bins in axis 2 (the y-axis)
boundaries1	the bin boundaries used in axis 1 (the x-axis)
boundaries2	the bin boundaries used in axis 2 (the y-axis)

Definition at line 1014 of file SampleHandlerFD.cpp.

{
  _hPDF1D->Reset();
  _hPDF1D->SetBins(static_cast<int>(nbins1),boundaries1);
  dathist->SetBins(static_cast<int>(nbins1),boundaries1);
 
  _hPDF2D->Reset();
  _hPDF2D->SetBins(static_cast<int>(nbins1),boundaries1,static_cast<int>(nbins2),boundaries2);
  dathist2d->SetBins(static_cast<int>(nbins1),boundaries1,static_cast<int>(nbins2),boundaries2);
  
  //Set the number of X and Y bins now
  SetupReweightArrays(Binning.XBinEdges.size() - 1, Binning.YBinEdges.size() - 1);
 
  FindNominalBinAndEdges2D();
}

Set2DBinning() [2/2]

void SampleHandlerFD::Set2DBinning ( std::vector< double > &  XVec, std::vector< double > &  YVec  )

inlineprotected

set the binning for 2D sample used for the likelihood calculation

Parameters

XVec	vector containing the binning in axis 1 (the x-axis)
YVec	vector containing the binning in axis 2 (the y-axis)

Definition at line 201 of file SampleHandlerFD.h.

{Set2DBinning(XVec.size()-1, XVec.data(), YVec.size()-1, YVec.data());};

SetupExperimentMC()

virtual int SampleHandlerFD::SetupExperimentMC ( )

protectedpure virtual

Experiment specific setup, returns the number of events which were loaded.

Implemented in PySampleHandlerFD.

SetupFDMC()

virtual void SampleHandlerFD::SetupFDMC ( )

protectedpure virtual

Function which translates experiment struct into core struct.

Implemented in PySampleHandlerFD.

SetupFunctionalParameters()

void SampleHandlerFD::SetupFunctionalParameters ( )

protectedvirtual

ETA - a function to setup and pass values to functional parameters where you need to pass a value to some custom reweight calc or engine.

Definition at line 639 of file SampleHandlerFD.cpp.

                                                {
  funcParsVec = ParHandler->GetFunctionalParametersFromSampleName(SampleName);
  // RegisterFunctionalParameters is implemented in experiment-specific code, 
  // which calls RegisterIndividualFuncPar to populate funcParsNamesMap, funcParsNamesVec, and funcParsFuncMap
  RegisterFunctionalParameters();
  funcParsMap.resize(funcParsNamesMap.size());
  funcParsGrid.resize(GetNEvents());
 
  // For every functional parameter in XsecCov that matches the name in funcParsNames, add it to the map
  for (FunctionalParameter & fp : funcParsVec) {
    for (std::string name : funcParsNamesVec) {
      if (fp.name == name) {
        MACH3LOG_INFO("Adding functional parameter: {} to funcParsMap with key: {}", fp.name, funcParsNamesMap[fp.name]);
        fp.funcPtr = &funcParsFuncMap[funcParsNamesMap[fp.name]];
        funcParsMap[static_cast<std::size_t>(funcParsNamesMap[fp.name])] = &fp;
        continue;
      }
    }
    // If we don't find a match, we need to throw an error
    if (funcParsMap[static_cast<std::size_t>(funcParsNamesMap[fp.name])] == nullptr) {
      MACH3LOG_ERROR("Functional parameter {} not found, did you define it in RegisterFunctionalParameters()?", fp.name);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  // Mostly the same as CalcXsecNormsBins
  // For each event, make a vector of pointers to the functional parameters
  for (std::size_t iEvent = 0; iEvent < static_cast<std::size_t>(GetNEvents()); ++iEvent) {
    // Now loop over the functional parameters and get a vector of enums corresponding to the functional parameters
    for (std::vector<FunctionalParameter>::iterator it = funcParsVec.begin(); it != funcParsVec.end(); ++it) {
      // Check whether the interaction modes match
      bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
      if (!ModeMatch) {
        MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
        continue;
      }
      // Now check whether within kinematic bounds
      bool IsSelected = true;
      if ((*it).hasKinBounds) {
        const auto& kinVars = (*it).KinematicVarStr;
        const auto& selection = (*it).Selection;
 
        for (std::size_t iKinPar = 0; iKinPar < kinVars.size(); ++iKinPar) {
          const double kinVal = ReturnKinematicParameter(kinVars[iKinPar], static_cast<int>(iEvent));
 
          bool passedAnyBound = false;
          const auto& boundsList = selection[iKinPar];
 
          for (const auto& bounds : boundsList) {
            if (kinVal > bounds[0] && kinVal <= bounds[1]) {
              passedAnyBound = true;
              break;
            }
          }
 
          if (!passedAnyBound) {
            MACH3LOG_TRACE("Event {}, missed kinematic check ({}) for dial {}",
                           iEvent, kinVars[iKinPar], (*it).name);
            IsSelected = false;
            break;
          }
        }
      }
      // Need to then break the event loop
      if(!IsSelected){
        MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", iEvent, (*it).name);
        continue;
      }
      auto funcparenum = funcParsNamesMap[(*it).name];
      funcParsGrid.at(iEvent).push_back(funcparenum);
    }
  }
  MACH3LOG_INFO("Finished setting up functional parameters");
}

SetupKinematicMap()

void SampleHandlerFD::SetupKinematicMap ( )

protected

Ensure Kinematic Map is setup and make sure it is initialised correctly.

Definition at line 257 of file SampleHandlerFD.cpp.

                                        {
// ************************************************
  if(KinematicParameters == nullptr || ReversedKinematicParameters == nullptr) {
    MACH3LOG_INFO("Map KinematicParameters or ReversedKinematicParameters hasn't been initialised");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  // KS: Ensure maps exist correctly
  for (const auto& pair : *KinematicParameters) {
    const auto& key = pair.first;
    const auto& value = pair.second;
 
    auto it = ReversedKinematicParameters->find(value);
    if (it == ReversedKinematicParameters->end() || it->second != key) {
      MACH3LOG_ERROR("Mismatch found: {} -> {} but {} -> {}",
                     key, value, value, (it != ReversedKinematicParameters->end() ? it->second : "NOT FOUND"));
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
}

SetupNormParameters()

void SampleHandlerFD::SetupNormParameters ( )

protected

Setup the norm parameters by assigning each event with bin.

Definition at line 777 of file SampleHandlerFD.cpp.

                                          {
// ***************************************************************************
  std::vector< std::vector< int > > xsec_norms_bins(GetNEvents());
 
  std::vector<NormParameter> norm_parameters = ParHandler->GetNormParsFromSampleName(GetSampleName());
 
  if(!ParHandler){
    MACH3LOG_ERROR("ParHandler is not setup!");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  // Assign xsec norm bins in MCSamples tree
  CalcNormsBins(norm_parameters, xsec_norms_bins);
 
  //DB Attempt at reducing impact of SystematicHandlerGeneric::calcReweight()
  int counter;
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) {
    counter = 0;
 
    MCSamples[iEvent].xsec_norm_pointers.resize(xsec_norms_bins[iEvent].size());
    for(auto const & norm_bin: xsec_norms_bins[iEvent]) {
      MCSamples[iEvent].xsec_norm_pointers[counter] = ParHandler->RetPointer(norm_bin);
      counter += 1;
    }
  }
}

SetupNuOscillatorPointers()

void SampleHandlerFD::SetupNuOscillatorPointers

(

)

Definition at line 1303 of file SampleHandlerFD.cpp.

                                                {
  for (unsigned int iEvent=0;iEvent<GetNEvents();iEvent++) {
    MCSamples[iEvent].osc_w_pointer = &M3::Unity;
    if (MCSamples[iEvent].isNC) {
      if (*MCSamples[iEvent].nupdg != *MCSamples[iEvent].nupdgUnosc) {
        MCSamples[iEvent].osc_w_pointer = &M3::Zero;
      } else {
        MCSamples[iEvent].osc_w_pointer = &M3::Unity;
      }
    } else {
      int InitFlav = M3::_BAD_INT_;
      int FinalFlav = M3::_BAD_INT_;
 
      InitFlav =  MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iEvent].nupdgUnosc));
      FinalFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iEvent].nupdg));
 
      if (InitFlav == M3::_BAD_INT_ || FinalFlav == M3::_BAD_INT_) {
        MACH3LOG_ERROR("Something has gone wrong in the mapping between MCSamples.nutype and the enum used within NuOscillator");
        MACH3LOG_ERROR("MCSamples.nupdgUnosc: {}", (*MCSamples[iEvent].nupdgUnosc));
        MACH3LOG_ERROR("InitFlav: {}", InitFlav);
        MACH3LOG_ERROR("MCSamples.nupdg: {}", (*MCSamples[iEvent].nupdg));
        MACH3LOG_ERROR("FinalFlav: {}", FinalFlav);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
 
      const int OscIndex = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
      //Can only happen if truecz has been initialised within the experiment specific code
      if (*(MCSamples[iEvent].rw_truecz) != M3::_BAD_DOUBLE_) {
        //Atmospherics
        MCSamples[iEvent].osc_w_pointer = Oscillator->GetNuOscillatorPointers(OscIndex, InitFlav, FinalFlav, FLOAT_T(*(MCSamples[iEvent].rw_etru)), FLOAT_T(*(MCSamples[iEvent].rw_truecz)));
      } else {
        //Beam
        MCSamples[iEvent].osc_w_pointer = Oscillator->GetNuOscillatorPointers(OscIndex, InitFlav, FinalFlav, FLOAT_T(*(MCSamples[iEvent].rw_etru)));
      }
    } // end if NC
  } // end loop over events
}

SetupReweightArrays()

void SampleHandlerFD::SetupReweightArrays ( const size_t  numberXBins, const size_t  numberYBins  )

protected

Initialise data, MC and W2 histograms.

Definition at line 913 of file SampleHandlerFD.cpp.

                                                                                            {
// ************************************************
  //Set the number of X and Y bins now
  Binning.nXBins = numberXBins;
  Binning.nYBins = numberYBins;
 
  SampleHandlerFD_array = new double*[Binning.nYBins];
  SampleHandlerFD_array_w2 = new double*[Binning.nYBins];
  SampleHandlerFD_data = new double*[Binning.nYBins];
  for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
    SampleHandlerFD_array[yBin] = new double[Binning.nXBins];
    SampleHandlerFD_array_w2[yBin] = new double[Binning.nXBins];
    SampleHandlerFD_data[yBin] = new double[Binning.nXBins];
    for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
      SampleHandlerFD_array[yBin][xBin] = 0.;
      SampleHandlerFD_array_w2[yBin][xBin] = 0.;
      SampleHandlerFD_data[yBin][xBin] = 0.;
    }
  }
}

SetupSampleBinning()

void SampleHandlerFD::SetupSampleBinning ( )

protected

wrapper to call set binning functions based on sample config info

Function to setup the binning of your sample histograms and the underlying arrays that get handled in fillArray() and fillArray_MP(). The Binning.XBinEdges are filled in the daughter class from the sample config file. This "passing" can be removed.

@function SampleHandlerFD::SetupSampleBinning()

Definition at line 307 of file SampleHandlerFD.cpp.

                                        {
// ************************************************
  MACH3LOG_INFO("Setting up Sample Binning");
  TString histname1d = (XVarStr).c_str();
  TString histname2d = (XVarStr+"_"+YVarStr).c_str();
  TString histtitle = GetTitle();
 
  //The binning here is arbitrary, now we get info from cfg so the
  //set1DBinning and set2Dbinning calls below will make the binning
  //to be what we actually want
  _hPDF1D   = new TH1D("h" + histname1d + GetTitle(),histtitle, 1, 0, 1);
  dathist   = new TH1D("d" + histname1d + GetTitle(),histtitle, 1, 0, 1);
  _hPDF2D   = new TH2D("h" + histname2d + GetTitle(),histtitle, 1, 0, 1, 1, 0, 1);
  dathist2d = new TH2D("d" + histname2d + GetTitle(),histtitle, 1, 0, 1, 1, 0, 1);
 
  _hPDF1D->GetXaxis()->SetTitle(XVarStr.c_str());
  dathist->GetXaxis()->SetTitle(XVarStr.c_str());
  _hPDF2D->GetXaxis()->SetTitle(XVarStr.c_str());
  dathist2d->GetXaxis()->SetTitle(XVarStr.c_str());
 
  //A string to store the binning for a nice print out
  std::string XBinEdgesStr = "";
  std::string YBinEdgesStr = "";
 
  for(auto XBinEdge : Binning.XBinEdges){
    XBinEdgesStr += std::to_string(XBinEdge);
    XBinEdgesStr += ", ";
  }
  MACH3LOG_INFO("XBinning:");
  MACH3LOG_INFO("{}", XBinEdgesStr);
  
  //And now the YBin Edges
  for(auto YBinEdge : Binning.YBinEdges){
    YBinEdgesStr += std::to_string(YBinEdge);
    YBinEdgesStr += ", ";
  }
  MACH3LOG_INFO("YBinning:");
  MACH3LOG_INFO("{}", YBinEdgesStr);
  
  //Check whether you are setting up 1D or 2D binning
  if(nDimensions == 1){
    MACH3LOG_INFO("Setting up 1D binning with {}", XVarStr);
    Set1DBinning(Binning.XBinEdges);
  }
  else if(nDimensions == 2){
    MACH3LOG_INFO("Setting up 2D binning with {} and {}", XVarStr, YVarStr);
    Set2DBinning(Binning.XBinEdges, Binning.YBinEdges);
  }
  else{
    MACH3LOG_ERROR("Number of dimensions is not 1 or 2, this is unsupported at the moment");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

SetupSplines()

virtual void SampleHandlerFD::SetupSplines ( )

protectedpure virtual

initialise your splineXX object and then use InitialiseSplineObject to conviently setup everything up

Todo:

abstract the spline initialisation completely to core

Implemented in PySampleHandlerFD.

SetupWeightPointers()

virtual void SampleHandlerFD::SetupWeightPointers ( )

protectedpure virtual

DB Function to determine which weights apply to which types of samples pure virtual!!

Implemented in PySampleHandlerFD.

Member Data Documentation

_hPDF1D

TH1D* SampleHandlerFD::_hPDF1D

protected

Definition at line 386 of file SampleHandlerFD.h.

_hPDF2D

TH2D* SampleHandlerFD::_hPDF2D

protected

Definition at line 387 of file SampleHandlerFD.h.

_modeNomWeightMap

std::unordered_map<std::string, double> SampleHandlerFD::_modeNomWeightMap

protected

Definition at line 370 of file SampleHandlerFD.h.

Binning

SampleBinningInfo SampleHandlerFD::Binning

protected

KS: This stores binning information, in future could be come vector to store binning for every used sample.

Definition at line 302 of file SampleHandlerFD.h.

dathist

TH1D* SampleHandlerFD::dathist

protected

Definition at line 382 of file SampleHandlerFD.h.

dathist2d

TH2D* SampleHandlerFD::dathist2d

protected

Definition at line 383 of file SampleHandlerFD.h.

FileToFinalPDGMap

std::unordered_map<std::string, NuPDG> SampleHandlerFD::FileToFinalPDGMap

private

Definition at line 396 of file SampleHandlerFD.h.

FileToInitPDGMap

std::unordered_map<std::string, NuPDG> SampleHandlerFD::FileToInitPDGMap

private

Definition at line 395 of file SampleHandlerFD.h.

FirstTimeW2

bool SampleHandlerFD::FirstTimeW2

protected

KS:Super hacky to update W2 or not.

Definition at line 378 of file SampleHandlerFD.h.

funcParsFuncMap

std::unordered_map<int, FuncParFuncType> SampleHandlerFD::funcParsFuncMap

protected

HH - a map that relates the funcpar enum to pointer of the actual function.

Definition at line 241 of file SampleHandlerFD.h.

funcParsGrid

std::vector<std::vector<int> > SampleHandlerFD::funcParsGrid

protected

HH - a grid of vectors of enums for each sample and event.

Definition at line 243 of file SampleHandlerFD.h.

funcParsMap

std::vector<FunctionalParameter *> SampleHandlerFD::funcParsMap

protected

HH - a map that relates the funcpar enum to pointer of FuncPars struct HH - Changed to a vector of pointers since it's faster than unordered_map and we are using ints as keys.

Definition at line 238 of file SampleHandlerFD.h.

funcParsNamesMap

std::unordered_map<std::string, int> SampleHandlerFD::funcParsNamesMap

protected

HH - a map that relates the name of the functional parameter to funcpar enum.

Definition at line 233 of file SampleHandlerFD.h.

funcParsNamesVec

std::vector<std::string> SampleHandlerFD::funcParsNamesVec = {}

protected

HH - a vector of string names for each functional parameter.

Definition at line 245 of file SampleHandlerFD.h.

funcParsVec

std::vector<FunctionalParameter> SampleHandlerFD::funcParsVec

protected

HH - a vector that stores all the FuncPars struct.

Definition at line 230 of file SampleHandlerFD.h.

KinematicParameters

const std::unordered_map<std::string, int>* SampleHandlerFD::KinematicParameters

protected

Mapping between string and kinematic enum.

Definition at line 349 of file SampleHandlerFD.h.

KinematicVectors

const std::unordered_map<std::string, int>* SampleHandlerFD::KinematicVectors

protected

Definition at line 354 of file SampleHandlerFD.h.

mc_files

std::vector<std::string> SampleHandlerFD::mc_files

protected

names of mc files associated associated with this object

Definition at line 366 of file SampleHandlerFD.h.

MCSamples

std::vector<FarDetectorCoreInfo> SampleHandlerFD::MCSamples

protected

Definition at line 314 of file SampleHandlerFD.h.

nDimensions

int SampleHandlerFD::nDimensions = M3::_BAD_INT_

protected

Keep track of the dimensions of the sample binning.

Definition at line 327 of file SampleHandlerFD.h.

OscChannels

std::vector<OscChannelInfo> SampleHandlerFD::OscChannels

protected

Definition at line 315 of file SampleHandlerFD.h.

Oscillator

std::shared_ptr<OscillationHandler> SampleHandlerFD::Oscillator

protected

Contains all your binned splines and handles the setup and the returning of weights from spline evaluations.

Definition at line 173 of file SampleHandlerFD.h.

ParHandler

ParameterHandlerGeneric* SampleHandlerFD::ParHandler = nullptr

protected

Definition at line 322 of file SampleHandlerFD.h.

ReversedKinematicParameters

const std::unordered_map<int, std::string>* SampleHandlerFD::ReversedKinematicParameters

protected

Mapping between kinematic enum and string.

Definition at line 351 of file SampleHandlerFD.h.

ReversedKinematicVectors

const std::unordered_map<int, std::string>* SampleHandlerFD::ReversedKinematicVectors

protected

Definition at line 355 of file SampleHandlerFD.h.

SampleHandlerFD_array

double** SampleHandlerFD::SampleHandlerFD_array

protected

DB Array to be filled after reweighting.

Definition at line 305 of file SampleHandlerFD.h.

SampleHandlerFD_array_w2

double** SampleHandlerFD::SampleHandlerFD_array_w2

protected

KS Array used for MC stat.

Definition at line 307 of file SampleHandlerFD.h.

SampleHandlerFD_data

double** SampleHandlerFD::SampleHandlerFD_data

protected

DB Array to be filled in AddData.

Definition at line 309 of file SampleHandlerFD.h.

SampleManager

std::unique_ptr<manager> SampleHandlerFD::SampleManager

protected

The manager object used to read the sample yaml file.

Definition at line 359 of file SampleHandlerFD.h.

SampleName

std::string SampleHandlerFD::SampleName

protected

A unique ID for each sample based on powers of two for quick binary operator comparisons.

Definition at line 329 of file SampleHandlerFD.h.

SampleTitle

std::string SampleHandlerFD::SampleTitle

protected

the name of this sample e.g."muon-like"

Definition at line 332 of file SampleHandlerFD.h.

Selection

std::vector< KinematicCut > SampleHandlerFD::Selection

protected

a way to store selection cuts which you may push back in the get1DVar functions most of the time this is just the same as StoredSelection

Definition at line 345 of file SampleHandlerFD.h.

spline_files

std::vector<std::string> SampleHandlerFD::spline_files

protected

names of spline files associated associated with this object

Definition at line 368 of file SampleHandlerFD.h.

SplineHandler

std::unique_ptr<BinnedSplineHandler> SampleHandlerFD::SplineHandler

protected

Contains all your binned splines and handles the setup and the returning of weights from spline evaluations.

Definition at line 170 of file SampleHandlerFD.h.

StoredSelection

std::vector< KinematicCut > SampleHandlerFD::StoredSelection

protected

What gets pulled from config options, these are constant after loading in this is of length 3: 0th index is the value, 1st is lower bound, 2nd is upper bound.

Definition at line 342 of file SampleHandlerFD.h.

THStackLeg

TLegend* SampleHandlerFD::THStackLeg = nullptr

protected

DB Miscellaneous Variables.

Definition at line 374 of file SampleHandlerFD.h.

UpdateW2

bool SampleHandlerFD::UpdateW2

protected

KS:Super hacky to update W2 or not.

Definition at line 380 of file SampleHandlerFD.h.

XVarStr

std::string SampleHandlerFD::XVarStr

protected

the strings associated with the variables used for the binning e.g. "RecoNeutrinoEnergy"

Definition at line 208 of file SampleHandlerFD.h.

YVarStr

std::string SampleHandlerFD::YVarStr

protected

Definition at line 208 of file SampleHandlerFD.h.

---

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

• Samples/SampleHandlerFD.h
• Samples/SampleHandlerFD.cpp