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

virtual	~SampleHandlerFD ()
	destructor More...

int	GetNDim () const
	DB Function to differentiate 1D or 2D binning. More...

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

std::string	GetTitle () const override

std::string	GetXBinVarName () const

std::string	GetYBinVarName () const

void	PrintIntegral (const TString &OutputName="/dev/null", const int WeightStyle=0, const 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. More...

TH1 *	GetMCHist (const int Dimension)
	Get MC histogram. More...

TH1 *	GetW2Hist (const int Dimension)
	Get W2 histogram. More...

TH1 *	GetDataHist (const int Dimension)
	Get Data histogram. More...

void	Reweight () override

M3::float_t	GetEventWeight (const int iEntry) const

void	InitialiseNuOscillatorObjects ()
	including Dan's magic NuOscillator More...

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

std::string	ReturnStringFromKinematicParameter (const int KinematicVariable) const
	ETA function to generically convert a kinematic type from xsec cov to a string. More...

void	SaveAdditionalInfo (TDirectory *Dir) override
	Store additional info in a chan. More...

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

Public Member Functions inherited from SampleHandlerBase
	SampleHandlerBase ()
	The main constructor. More...

virtual	~SampleHandlerBase ()
	destructor More...

virtual M3::int_t	GetNsamples ()

virtual double	GetSampleLikelihood (const int isample)

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

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

unsigned int	GetNEvents () const

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	GetPoissonLLH (const double data, const double mc) const
	Calculate test statistic for a single bin using Poisson. More...

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

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

Protected Member Functions
virtual void	SetupWeightPointers ()=0
	DB Function to determine which weights apply to which types of samples pure virtual!! More...

void	SetupKinematicMap ()
	Ensure Kinematic Map is setup and make sure it is initialised correctly. More...

virtual void	SetupSplines ()=0
	initialise your splineXX object and then use InitialiseSplineObject to conviently setup everything up More...

virtual void	Init ()=0
	Initialise any variables that your experiment specific SampleHandler needs. More...

virtual int	SetupExperimentMC ()=0
	Experiment specific setup, returns the number of events which were loaded. More...

virtual void	SetupFDMC ()=0
	Function which translates experiment struct into core struct. More...

void	Initialise ()
	Function which does a lot of the lifting regarding the workflow in creating different MC objects. More...

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

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

void	Set2DBinning (size_t nbins1, double boundaries1, size_t nbins2, double boundaries2)
	set the binning for 2D sample used for the likelihood calculation More...

void	Set1DBinning (std::vector< double > &XVec)
	set the binning for 1D sample used for the likelihood calculation More...

void	Set2DBinning (std::vector< double > &XVec, std::vector< double > &YVec)
	set the binning for 2D sample used for the likelihood calculation More...

void	SetupSampleBinning ()
	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. More...

void	SetupReweightArrays (const size_t numberXBins, const size_t numberYBins)
	Initialise data, MC and W2 histograms. More...

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

void	RegisterIndividualFunctionalParameter (const std::string &fpName, int fpEnum, FuncParFuncType fpFunc)
	HH - a helper function for RegisterFunctionalParameter. More...

virtual void	RegisterFunctionalParameters ()=0
	HH - a experiment-specific function where the maps to actual functions are set up. More...

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

virtual void	ApplyShifts (int iEvent)
	ETA - generic function applying shifts. More...

bool	IsEventSelected (const int iEvent)
	DB Function which determines if an event is selected, where Selection double looks like {{ND280KinematicTypes Var1, douuble LowBound}. More...

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

virtual void	resetShifts (int iEvent)
	HH - reset the shifted values to the original values. More...

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

M3::float_t	CalcWeightSpline (const FarDetectorCoreInfo *MCEvent) const
	Calculate the spline weight for a given event. More...

M3::float_t	CalcWeightNorm (const FarDetectorCoreInfo *MCEvent) const
	Calculate the norm weight for a given event. More...

virtual void	CalcWeightFunc (int iEvent)
	Calculate weights for function parameters. More...

virtual double	ReturnKinematicParameter (std::string KinematicParamter, int iEvent)=0
	Return the value of an associated kinematic parameter for an event. More...

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

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

void	SetupNormParameters ()
	Setup the norm parameters by assigning each event with bin. More...

void	FillMCHist (const int Dimension)
	Fill a histogram with the event-level information used in the fit. More...

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

void	ResetHistograms ()
	Helper function to reset histograms. More...

void	InitialiseSingleFDMCObject ()
	function to create the member of the FarDetectorInfo struct so they are the appropriate size. More...

void	InitialiseSplineObject ()

NuPDG	GetInitPDGFromFileName (const std::string &FileName) const
	Retrieve the initial neutrino PDG code associated with a given input file name. More...

NuPDG	GetFinalPDGFromFileName (const std::string &FileName) const
	Retrieve the final neutrino PDG code associated with a given input file name. More...

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

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

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

Protected Attributes
std::unique_ptr< BinnedSplineHandler >	SplineHandler
	Contains all your binned splines and handles the setup and the returning of weights from spline evaluations. More...

std::shared_ptr< OscillationHandler >	Oscillator
	Contains all your binned splines and handles the setup and the returning of weights from spline evaluations. More...

std::vector< FunctionalParameter >	funcParsVec
	HH - a vector that stores all the FuncPars struct. More...

std::unordered_map< std::string, int >	funcParsNamesMap
	HH - a map that relates the name of the functional parameter to funcpar enum. More...

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

std::unordered_map< int, FuncParFuncType >	funcParsFuncMap
	HH - a map that relates the funcpar enum to pointer of the actual function. More...

std::vector< std::vector< int > >	funcParsGrid
	HH - a grid of vectors of enums for each sample and event. More...

std::vector< std::string >	funcParsNamesVec = {}
	HH - a vector of string names for each functional parameter. More...

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

double *	SampleHandlerFD_array
	DB Array to be filled after reweighting. More...

double *	SampleHandlerFD_array_w2
	KS Array used for MC stat. More...

double *	SampleHandlerFD_data
	DB Array to be filled in AddData. More...

std::vector< FarDetectorCoreInfo >	MCSamples
	Stores information about every MC event. More...

SampleInfo	SampleDetails
	Stores info about currently initialised sample. More...

std::vector< OscChannelInfo >	OscChannels
	Stores info about oscillation channel for a single sample. More...

ParameterHandlerGeneric *	ParHandler = nullptr

std::string	SampleName
	A unique ID for each sample based on which we can define what systematic should be applied. More...

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

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

const std::unordered_map< std::string, int > *	KinematicParameters
	Mapping between string and kinematic enum. More...

const std::unordered_map< int, std::string > *	ReversedKinematicParameters
	Mapping between kinematic enum and string. More...

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

std::vector< std::string >	mc_files
	names of mc files associated associated with this object More...

std::vector< std::string >	spline_files
	names of spline files associated associated with this object More...

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

TLegend *	THStackLeg = nullptr
	DB Miscellaneous Variables. More...

bool	FirstTimeW2
	KS:Super hacky to update W2 or not. More...

bool	UpdateW2
	KS:Super hacky to update W2 or not. More...

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

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

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

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

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

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

Private Types
enum	FDPlotType { kModePlot = 0 , kOscChannelPlot = 1 }

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.

Member Enumeration Documentation

◆ FDPlotType

enum SampleHandlerFD::FDPlotType

private

Enumerator
kModePlot
kOscChannelPlot

Definition at line 391 of file SampleHandlerFD.h.

                   {
     kModePlot = 0,
     kOscChannelPlot = 1
   };

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_;
   }
  
   KinematicParameters = nullptr;
   ReversedKinematicParameters = nullptr;
   KinematicVectors = nullptr;
   ReversedKinematicVectors = nullptr;
  
   SampleHandlerFD_array = nullptr;
   SampleHandlerFD_data = nullptr;
   SampleHandlerFD_array_w2 = nullptr;
   SampleName = "";
   SampleManager = std::make_unique<manager>(ConfigFileName.c_str());
  
   // Variables related to MC stat
   FirstTimeW2 = true;
   UpdateW2 = false;
 }

◆ ~SampleHandlerFD()

SampleHandlerFD::~SampleHandlerFD ( )

virtual

destructor

Definition at line 46 of file SampleHandlerFD.cpp.

                                   {
   MACH3LOG_DEBUG("I'm deleting SampleHandlerFD");
   
   if (SampleHandlerFD_array != nullptr) delete[] SampleHandlerFD_array;
   if (SampleHandlerFD_array_w2 != nullptr) 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;
 }

Member Function Documentation

◆ AddData() [1/4]

void SampleHandlerFD::AddData ( std::vector< double > & data )

Definition at line 1061 of file SampleHandlerFD.cpp.

                                                      {
   if (SampleDetails.dathist == nullptr) {
     MACH3LOG_ERROR("Data hist hasn't been initialised yet");
     throw MaCh3Exception(__FILE__, __LINE__);
   }
  
   SampleDetails.dathist2d = nullptr;
   SampleDetails.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){
     SampleDetails.dathist->Fill(data_point);
   }
  
   if(SampleHandlerFD_data == nullptr) {
     MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
     throw MaCh3Exception(__FILE__, __LINE__);
   }
   // Assuming nBins == nXBins here, because you have 1D data
   for (size_t bin = 0; bin < Binning.nXBins; ++bin) {
     // ROOT histograms are 1-based, so bin index + 1
     SampleHandlerFD_data[bin] = SampleDetails.dathist->GetBinContent(static_cast<int>(bin + 1));
   }
 }

◆ AddData() [2/4]

void SampleHandlerFD::AddData ( std::vector< std::vector< double > > & data )

Definition at line 1093 of file SampleHandlerFD.cpp.

                                                                    {
   if (SampleDetails.dathist2d == nullptr) {
     MACH3LOG_ERROR("Data hist hasn't been initialised yet");
     throw MaCh3Exception(__FILE__, __LINE__);
   }
  
   SampleDetails.dathist = nullptr;
   SampleDetails.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++) {
     SampleDetails.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
       const int idx = Binning.GetBinSafe(xBin, yBin);
       SampleHandlerFD_data[idx] = SampleDetails.dathist2d->GetBinContent(static_cast<int>(xBin + 1), static_cast<int>(yBin + 1));
     }
   }
 }

◆ AddData() [3/4]

void SampleHandlerFD::AddData ( TH1D * Data )

Definition at line 1129 of file SampleHandlerFD.cpp.

                                         {
   MACH3LOG_INFO("Adding 1D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
   if (SampleDetails.dathist != nullptr) {
     delete SampleDetails.dathist;
   }
   SampleDetails.dathist2d = nullptr;
   SampleDetails.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 bin = 0; bin < Binning.nXBins; ++bin) {
     // ROOT histograms are 1-based, so bin index + 1
     SampleHandlerFD_data[bin] = SampleDetails.dathist->GetBinContent(static_cast<int>(bin + 1));
   }
 }

◆ AddData() [4/4]

void SampleHandlerFD::AddData ( TH2D * Data )

Definition at line 1154 of file SampleHandlerFD.cpp.

                                         {
   MACH3LOG_INFO("Adding 2D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
   if (SampleDetails.dathist2d != nullptr) {
     delete SampleDetails.dathist2d;
   }
   SampleDetails.dathist2d = static_cast<TH2D*>(Data->Clone());
   SampleDetails.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
       const int idx = Binning.GetBinSafe(xBin, yBin);
       SampleHandlerFD_data[idx] = SampleDetails.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 675 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 767 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.

259 {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 715 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 698 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;
 }

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

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 SampleDetails._hPDF2D->GetBinContent(2,2) {Noticing the offset}

Definition at line 398 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) {
       const int GlobalBin = Binning.GetBin(XBinToFill, YBinToFill);
       SampleHandlerFD_array[GlobalBin] += totalweight;
       if (FirstTimeW2) SampleHandlerFD_array_w2[GlobalBin] += totalweight*totalweight;
     }
   }
 }

◆ FillMCHist()

void SampleHandlerFD::FillMCHist ( const int Dimension )

protected

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

Definition at line 263 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
   if(Dimension == 1){
     SampleDetails._hPDF1D->Reset();
     for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
       for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
         const int idx = Binning.GetBinSafe(xBin, yBin);
         SampleDetails._hPDF1D->SetBinContent(idx + 1, SampleHandlerFD_array[idx]);
       }
     }
   } else if (Dimension == 2) {
     SampleDetails._hPDF2D->Reset();
     for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
       for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
         const int idx = Binning.GetBinSafe(xBin, yBin);
         SampleDetails._hPDF2D->SetBinContent(static_cast<int>(xBin + 1), static_cast<int>(yBin + 1), SampleHandlerFD_array[idx]);
       }
     }
   } else {
     MACH3LOG_ERROR("Asking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
     throw MaCh3Exception(__FILE__, __LINE__);
   }
 }

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

Definition at line 1302 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(GetNDim()){
       case 1:
         EventSplines = SplineHandler->GetEventSplines(GetTitle(), OscIndex, int(*(MCSamples[j].mode)), *(MCSamples[j].rw_etru), *(MCSamples[j].x_var), 0.);
         break;
       case 2:
         EventSplines = SplineHandler->GetEventSplines(GetTitle(), 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 = {}", GetNDim());
         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 917 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(GetXBinVarName(), 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 = SampleDetails._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.InitialiseBinMigrationLookUp();
 }

◆ FindNominalBinAndEdges2D()

void SampleHandlerFD::FindNominalBinAndEdges2D ( )

protected

Definition at line 958 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(GetXBinVarName(), event_i);
     MCSamples[event_i].y_var = GetPointerToKinematicParameter(GetYBinVarName(), 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 = SampleDetails._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_;
     SampleDetails._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.InitialiseBinMigrationLookUp();
 }

◆ 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 1515 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 1287 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 385 of file SampleHandlerFD.h.

385 {return FileToFinalPDGMap.at(FileName);}

SampleHandlerFD::FileToFinalPDGMap

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

Definition: SampleHandlerFD.h:389

◆ 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 383 of file SampleHandlerFD.h.

383 {return FileToInitPDGMap.at(FileName);}

SampleHandlerFD::FileToInitPDGMap

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

Definition: SampleHandlerFD.h:388

◆ GetLikelihood()

double SampleHandlerFD::GetLikelihood ( )

overridevirtual

DB Multi-threaded GetLikelihood.

Implements SampleHandlerBase.

Definition at line 1338 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 reduction(+:negLogL)
   #endif
   for (size_t idx = 0; idx < Binning.nBins; ++idx)
   {
     const double DataVal = SampleHandlerFD_data[idx];
     const double MCPred = SampleHandlerFD_array[idx];
     const double w2 = SampleHandlerFD_array_w2[idx];
  
     //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 2047 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 185 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 1182 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, GetNOscChannels());
  
   if (!EqualBinningPerOscChannel) {
     for(int iChannel = 0; iChannel < GetNOscChannels(); 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 1361 of file SampleHandlerFD.cpp.

                                                  {
   MCSamples.resize(nEvents);
 }

◆ InitialiseSplineObject()

void SampleHandlerFD::InitialiseSplineObject ( )

protected

Definition at line 1400 of file SampleHandlerFD.cpp.

                                              {
   std::vector<std::string> spline_filepaths;
   for(int iChannel = 0 ; iChannel < GetNOscChannels(); iChannel++){
     spline_filepaths.push_back(spline_files[iChannel]);
   }
  
   //Keep a track of the spline variables
   std::vector<std::string> SplineVarNames = {"TrueNeutrinoEnergy"};
   if(GetXBinVarName().length() > 0){
     SplineVarNames.push_back(GetXBinVarName());
   }
   if(GetYBinVarName().length() > 0){
     SplineVarNames.push_back(GetYBinVarName());
   }
   
   bool LoadSplineFile = GetFromManager<bool>(SampleManager->raw()["InputFiles"]["LoadSplineFile"], false, __FILE__, __LINE__);
   bool PrepSplineFile = GetFromManager<bool>(SampleManager->raw()["InputFiles"]["PrepSplineFile"], false, __FILE__, __LINE__);
   auto SplineFileName = GetFromManager<std::string>(SampleManager->raw()["InputFiles"]["SplineFileName"],
                                                     (SampleName + "_SplineFile.root"), __FILE__, __LINE__);
   if(!LoadSplineFile) {
     SplineHandler->AddSample(SampleName, GetTitle(), spline_filepaths, SplineVarNames);
     SplineHandler->CountNumberOfLoadedSplines(false, 1);
     SplineHandler->TransferToMonolith();
     if(PrepSplineFile) SplineHandler->PrepareSplineFile(SplineFileName);
   } else {
     // KS: Skip default spline loading and use flattened spline format allowing to read stuff much faster
     SplineHandler->LoadSplineFile(SplineFileName);
   }
   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 334 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 349 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 1714 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.

219 {};

◆ PrintIntegral()

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

Definition at line 1828 of file SampleHandlerFD.cpp.

                                                                                                                           {
   constexpr 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(GetXBinVarName(),1,i,WeightStyle);
     } else {
       IntegralList[i] = CastVector<TH2, TH1>(ReturnHistsBySelection2D(GetXBinVarName(),GetYBinVarName(),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 (int i = 0;i < GetNOscChannels(); i++) {
     table_headings += fmt::format(" {:<17} |", GetFlavourName(i));
     table_footline += "--------------------";
     if (printToFile) {outfile << "&" << std::setw(space) << OscChannels[i].flavourName_Latex << " ";}
     if (printToCSV)  {outcsv << GetFlavourName(i) << ",";}
   }
   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();
   }
   // KS: Clean memory we could use smart pointers in future
   CleanContainer(IntegralList);
 }

◆ ReadSampleConfig()

void SampleHandlerFD::ReadSampleConfig ( )

Definition at line 57 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
   SampleDetails.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
   SampleDetails.nDimensions = 0;
   SampleDetails.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(GetXBinVarName().length() > 0){
     SampleDetails.nDimensions++;
   } else{
     MACH3LOG_ERROR("Please specify an X-variable string in sample config {}", SampleManager->GetFileName());
     throw MaCh3Exception(__FILE__, __LINE__);
   }
   
   SampleDetails.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(GetYBinVarName().length() > 0){
     if(GetXBinVarName().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__);
     }
     SampleDetails.nDimensions++;
   }
   
   if(GetNDim() == 0){
     MACH3LOG_ERROR("Error setting up the sample binning");
     MACH3LOG_ERROR("Number of dimensions is {}", GetNDim());
     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", GetNDim());
   }
   
   //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      = GetNOscChannels();
  
     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 581 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 568 of file SampleHandlerFD.cpp.

                                       {
 // **************************************************  
   //DB Reset values stored in PDF array to 0.
   // Don't openMP this; no significant gain
   #ifdef MULTITHREAD
   #pragma omp simd
   #endif
   for (size_t i = 0; i < Binning.nBins; ++i) {
     SampleHandlerFD_array[i] = 0.;
     if (FirstTimeW2) SampleHandlerFD_array_w2[i] = 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.

228 {(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 associated 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 1680 of file SampleHandlerFD.cpp.

                                                                                                       {
 // ************************************************
   // If x or y variable return used binning
   if(KinematicParameter == GetXBinVarName()) {
     return Binning.XBinEdges;
   } else if (KinematicParameter == GetYBinVarName()) {
     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.

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

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

◆ ReturnKinematicVectorFromString()

int SampleHandlerFD::ReturnKinematicVectorFromString ( const std::string & KinematicStr ) const

Definition at line 1654 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 1666 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 368 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 1367 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(GetXBinVarName().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(GetXBinVarName().c_str());
     data_hist->GetYaxis()->SetTitle(GetYBinVarName().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 899 of file SampleHandlerFD.cpp.

 {
   SampleDetails._hPDF1D->Reset();
   SampleDetails._hPDF1D->SetBins(static_cast<int>(nbins),boundaries);
   SampleDetails.dathist->SetBins(static_cast<int>(nbins),boundaries);
  
   Binning.YBinEdges = {-1e8, 1e8};
  
   SampleDetails._hPDF2D->Reset();
   SampleDetails._hPDF2D->SetBins(static_cast<int>(nbins),boundaries,1, Binning.YBinEdges.data());
   SampleDetails.dathist2d->SetBins(static_cast<int>(nbins),boundaries,1, Binning.YBinEdges.data());
  
   //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.

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

SampleHandlerFD::Set1DBinning

void Set1DBinning(size_t nbins, double *boundaries)

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

Definition: SampleHandlerFD.cpp:899

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

 {
   SampleDetails._hPDF1D->Reset();
   SampleDetails._hPDF1D->SetBins(static_cast<int>(nbins1),boundaries1);
   SampleDetails.dathist->SetBins(static_cast<int>(nbins1),boundaries1);
  
   SampleDetails._hPDF2D->Reset();
   SampleDetails._hPDF2D->SetBins(static_cast<int>(nbins1),boundaries1,static_cast<int>(nbins2),boundaries2);
   SampleDetails.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.

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

SampleHandlerFD::Set2DBinning

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

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

Definition: SampleHandlerFD.cpp:941

◆ 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 600 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 240 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 738 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 1235 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 874 of file SampleHandlerFD.cpp.

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

◆ SetupSampleBinning()

void SampleHandlerFD::SetupSampleBinning ( )

protected

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.

Definition at line 294 of file SampleHandlerFD.cpp.

                                         {
 // ************************************************
   MACH3LOG_INFO("Setting up Sample Binning");
   SampleDetails.InitialiseHistograms();
  
   //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(GetNDim() == 1){
     MACH3LOG_INFO("Setting up {}D binning with {}", GetNDim(), GetXBinVarName());
     Set1DBinning(Binning.XBinEdges);
   }
   else if(GetNDim() == 2){
     MACH3LOG_INFO("Setting up {}D binning with {} and {}", GetNDim(), GetXBinVarName(), GetYBinVarName());
     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

◆ _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 305 of file SampleHandlerFD.h.

◆ FileToFinalPDGMap

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

private

Definition at line 389 of file SampleHandlerFD.h.

◆ FileToInitPDGMap

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

private

Definition at line 388 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

Stores information about every MC event.

Definition at line 317 of file SampleHandlerFD.h.

◆ OscChannels

std::vector<OscChannelInfo> SampleHandlerFD::OscChannels

protected

Stores info about oscillation channel for a single sample.

Definition at line 321 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 328 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.

◆ SampleDetails

SampleInfo SampleHandlerFD::SampleDetails

protected

Stores info about currently initialised sample.

Definition at line 319 of file SampleHandlerFD.h.

◆ SampleHandlerFD_array

double* SampleHandlerFD::SampleHandlerFD_array

protected

DB Array to be filled after reweighting.

Definition at line 308 of file SampleHandlerFD.h.

◆ SampleHandlerFD_array_w2

double* SampleHandlerFD::SampleHandlerFD_array_w2

protected

KS Array used for MC stat.

Definition at line 310 of file SampleHandlerFD.h.

◆ SampleHandlerFD_data

double* SampleHandlerFD::SampleHandlerFD_data

protected

DB Array to be filled in AddData.

Definition at line 312 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 which we can define what systematic should be applied.

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.

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

Samples/SampleHandlerFD.h
Samples/SampleHandlerFD.cpp

Public Member Functions

Protected Member Functions

Protected Attributes

Private Types

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ FDPlotType

Constructor & Destructor Documentation

◆ SampleHandlerFD()

◆ ~SampleHandlerFD()

Member Function Documentation

◆ AddData() [1/4]

◆ AddData() [2/4]

◆ AddData() [3/4]

◆ AddData() [4/4]

◆ ApplyShifts()

◆ CalcNormsBins()

◆ CalcWeightFunc()

◆ CalcWeightNorm()

◆ CalcWeightSpline()

◆ Fill2DSubEventHist()

◆ FillArray()

◆ FillMCHist()

◆ FillSplineBins()

◆ FindNominalBinAndEdges1D()

◆ FindNominalBinAndEdges2D()

◆ Get2DVarHist()

◆ GetEventWeight()

◆ GetFinalPDGFromFileName()

◆ GetInitPDGFromFileName()

◆ GetLikelihood()

◆ GetPointerToKinematicParameter() [1/2]

◆ GetPointerToKinematicParameter() [2/2]

◆ GetPointerToOscChannel()

◆ Init()

◆ Initialise()

◆ InitialiseNuOscillatorObjects()

◆ InitialiseSingleFDMCObject()

◆ InitialiseSplineObject()

◆ IsEventSelected()

◆ IsSubEventSelected()

◆ IsSubEventVarString()

◆ PrepFunctionalParameters()

◆ PrintIntegral()

◆ ReadSampleConfig()

◆ RegisterFunctionalParameters()

◆ RegisterIndividualFunctionalParameter()

◆ ResetHistograms()

◆ resetShifts()

◆ ReturnKinematicParameter() [1/2]

◆ ReturnKinematicParameter() [2/2]

◆ ReturnKinematicParameterBinning()

◆ ReturnKinematicVector() [1/2]

◆ ReturnKinematicVector() [2/2]

◆ ReturnKinematicVectorFromString()

◆ ReturnStringFromKinematicVector()

◆ Reweight()

◆ SaveAdditionalInfo()

◆ Set1DBinning() [1/2]

◆ Set1DBinning() [2/2]

◆ Set2DBinning() [1/2]

◆ Set2DBinning() [2/2]

◆ SetupExperimentMC()

◆ SetupFDMC()

◆ SetupFunctionalParameters()

◆ SetupKinematicMap()

◆ SetupNormParameters()

◆ SetupNuOscillatorPointers()

◆ SetupReweightArrays()

◆ SetupSampleBinning()

◆ SetupSplines()

◆ SetupWeightPointers()

Member Data Documentation

◆ _modeNomWeightMap

◆ Binning

◆ FileToFinalPDGMap

◆ FileToInitPDGMap

◆ FirstTimeW2

◆ funcParsFuncMap