SampleHandlerBase Class Reference

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

#include <Samples/SampleHandlerBase.h>

Inheritance diagram for SampleHandlerBase:

Collaboration diagram for SampleHandlerBase:

Public Member Functions
	SampleHandlerBase (std::string ConfigFileName, ParameterHandlerGeneric *xsec_cov, const std::shared_ptr< OscillationHandler > &OscillatorObj_=nullptr)
	Constructor. More...
virtual	~SampleHandlerBase ()
	destructor More...
int	GetNDim (const int Sample) const final
	DB Get what dimensionality binning for given sample has. More...
std::string	GetName () const final
	Get name for Sample Handler. More...
std::string	GetSampleTitle (const int Sample) const final
	Get fancy title for specified samples. More...
std::string	GetKinVarName (const int iSample, const int Dimension) const final
	Return Kinematic Variable name for specified sample and dimension for example "Reconstructed_Neutrino_Energy". More...
void	PrintIntegral (const int iSample, const TString &OutputName="/dev/null", const int WeightStyle=0, const TString &OutputCSVName="/dev/null")
	Computes and prints the integral breakdown of all modes and oscillation channels for a given sample. More...
void	AddData (const int Sample, TH1 *Data)
void	AddData (const int Sample, const std::vector< double > &Data_Array)
void	PrintRates (const bool DataOnly=false) final
	Helper function to print rates for the samples with LLH. More...
double	GetLikelihood () const override
	DB Multi-threaded GetLikelihood. More...
double	GetSampleLikelihood (const int isample) const override
	Get likelihood for single sample. More...
int	GetSampleIndex (const std::string &SampleTitle) const
	Get index of sample based on name. More...
const TH1 *	GetDataHist (const int Sample) final
	Get Data histogram. More...
const TH1 *	GetDataHist (const std::string &Sample)
const TH1 *	GetMCHist (const int Sample) final
	Get MC histogram. More...
const TH1 *	GetMCHist (const std::string &Sample)
const TH1 *	GetW2Hist (const int Sample) final
	Get W2 histogram. More...
const TH1 *	GetW2Hist (const std::string &Sample)
void	Reweight () override
	main routine modifying MC prediction based on proposed parameter values More...
M3::float_t	GetEventWeight (const int iEntry)
	Computes the total event weight for a given entry. More...
const M3::float_t *	GetNuOscillatorPointers (const int iEvent) const
int	GetNOscChannels (const int iSample) const final
	Get number of oscillation channels for a single sample. More...
std::string	GetFlavourName (const int iSample, const int iChannel) const final
std::unique_ptr< TH1 >	Get1DVarHist (const int iSample, const std::string &ProjectionVar, const std::vector< KinematicCut > &EventSelectionVec={}, int WeightStyle=0, const std::vector< KinematicCut > &SubEventSelectionVec={}) final
std::unique_ptr< TH2 >	Get2DVarHist (const int iSample, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &EventSelectionVec={}, int WeightStyle=0, const std::vector< KinematicCut > &SubEventSelectionVec={}) final
std::vector< KinematicCut >	BuildModeChannelSelection (const int iSample, const int kModeToFill, const int kChannelToFill) const
void	Fill1DSubEventHist (const int iSample, TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec={}, int WeightStyle=0)
void	Fill2DSubEventHist (const int iSample, TH2 *_h2DVar, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &SubEventSelectionVec={}, int WeightStyle=0)
std::unique_ptr< TH1 >	Get1DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_Str, const int kModeToFill=-1, const int kChannelToFill=-1, const int WeightStyle=0) final
std::unique_ptr< TH2 >	Get2DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_StrX, const std::string &ProjectionVar_StrY, const int kModeToFill=-1, const int kChannelToFill=-1, const int WeightStyle=0) final
std::unique_ptr< TH1 >	GetModeHist1D (const int iSample, int s, int m, int style=0)
std::unique_ptr< TH2 >	GetModeHist2D (const int iSample, int s, int m, int style=0)
std::vector< std::unique_ptr< TH1 > >	ReturnHistsBySelection1D (const int iSample, const std::string &KinematicProjection, const int Selection1, const int Selection2=-1, const int WeightStyle=0)
std::vector< std::unique_ptr< TH2 > >	ReturnHistsBySelection2D (const int iSample, const std::string &KinematicProjectionX, const std::string &KinematicProjectionY, const int Selection1, const int Selection2=-1, const int WeightStyle=0)
std::unique_ptr< THStack >	ReturnStackedHistBySelection1D (const int iSample, const std::string &KinematicProjection, const int Selection1, const int Selection2=-1, const int WeightStyle=0)
const TLegend *	ReturnStackHistLegend () const
	Return the legend used for stacked histograms with sample info. More...
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) final
	Store additional info in a chan. More...
int	ReturnKinematicVectorFromString (const std::string &KinematicStr) const
	JM: Convert a kinematic vector name to its corresponding integer ID. More...
std::string	ReturnStringFromKinematicVector (const int KinematicVariable) const
	JM: Convert a kinematic vector integer ID to its corresponding name as a string. More...
bool	IsSubEventVarString (const std::string &VarStr) const
	JM: Check if a kinematic parameter string corresponds to a subevent-level variable. More...
auto	GetDataArray () const
	Return array storing data entries for every bin. More...
auto	GetMCArray () const
	Return array storing MC entries for every bin. More...
auto	GetW2Array () const
	Return array storing W2 entries for every bin. More...
std::vector< double >	GetArrayForSample (const int Sample, std::vector< double > const &array) const
	Return a sub-array for a given sample. More...
std::vector< double >	GetDataArray (const int Sample) const
	Return array storing data entries for every bin. More...
std::vector< double >	GetMCArray (const int Sample) const
	Return array storing MC entries for every bin. More...
std::vector< double >	GetW2Array (const int Sample) const
	Return array storing W2 entries for single sample. More...
Public Member Functions inherited from SampleHandlerInterface
	SampleHandlerInterface ()
	The main constructor. More...
virtual	~SampleHandlerInterface ()
	destructor More...
virtual M3::int_t	GetNSamples ()
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
double	GetPoissonLLH (const double data, const double mc) const
	Calculate test statistic for a single bin using Poisson. More...
double	GetTestStatLLH (const double data, const double mc, const double w2) const
	Calculate test statistic for a single bin. Calculation depends on setting of fTestStatistic. Data and mc -> 0 cut-offs are defined in M3::LOW_MC_BOUND. More...
void	SetTestStatistic (TestStatistic testStat)
	Set the test statistic to be used when calculating the binned likelihoods. More...
TestStatistic	GetTestStatistic () const
	Get the test statistic used when calculating the binned likelihoods. More...

Protected Member Functions
void	InitialiseNuOscillatorObjects ()
	including Dan's magic NuOscillator More...
void	SetupNuOscillatorPointers ()
	Initialise pointer to oscillation weight to NuOscillator object. More...
void	ReadConfig ()
	Load information about sample handler and corresponding samples from config file. More...
void	LoadSingleSample (const int iSample, const YAML::Node &Settings)
	Initialise single sample from config file. More...
virtual void	AddAdditionalWeightPointers ()=0
	DB Function to determine which weights apply to which types of samples. 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	SetupMC ()=0
	Function which translates experiment struct into core struct. More...
virtual void	InititialiseData ()=0
	Function responsible for loading data from file or loading from file. More...
void	Initialise ()
	Function which does a lot of the lifting regarding the workflow in creating different MC objects. More...
void	SetSplinePointers ()
	Set pointers for each event to appropriate weights, for unbinned based on event number while for binned based on other kinematical properties. More...
std::vector< std::vector< int > >	GetSplineBins (int Event, BinnedSplineHandler *BinnedSpline, bool &ThrowCrititcal) const
	Retrieve the spline bin indices associated with a given event. More...
void	FindNominalBinAndEdges ()
void	SetBinning ()
	set the binning for 2D sample used for the likelihood calculation More...
void	SetupReweightArrays ()
	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 (const int iEvent)
	ETA - generic function applying shifts. More...
bool	IsEventSelected (const int iSample, const int iEvent) _noexcept_
	DB Function which determines if an event is selected based on KinematicCut. 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 (const int iEvent)
	HH - reset the shifted values to the original values. More...
virtual void	FinaliseShifts (const int iEvent)
	LP - Optionally calculate derived observables after all shifts have been applied. More...
void	CalcNormsBins (std::vector< NormParameter > &norm_parameters, std::vector< std::vector< int > > &norms_bins)
	Check whether a normalisation systematic affects an event or not. More...
template<typename ParT >
bool	PassesSelection (const ParT &Par, std::size_t iEvent)
M3::float_t	CalcWeightTotal (const EventInfo *_restrict_ MCEvent) const _noexcept_
	Calculate the total weight weight for a given event. More...
virtual void	CalcWeightFunc (const int iEvent)
	Calculate weights for function parameters. More...
double	ReturnKinematicParameter (const std::string &KinematicParameter, int iEvent) const
	Return the value of an associated kinematic parameter for an event. More...
virtual double	ReturnKinematicParameter (const int KinematicVariable, const int iEvent) const =0
std::vector< double >	ReturnKinematicVector (const std::string &KinematicParameter, const int iEvent) const
virtual std::vector< double >	ReturnKinematicVector (const int KinematicVariable, const int iEvent) const
std::vector< double >	ReturnKinematicParameterBinning (const int Sample, const std::string &KinematicParameter) const final
	Return the binning used to draw a kinematic parameter. More...
const double *	GetPointerToKinematicParameter (const std::string &KinematicParameter, int iEvent) const
virtual const double *	GetPointerToKinematicParameter (const int KinematicVariable, const int iEvent) const =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	SetupOscParameters ()
	Setup the osc parameters. More...
void	FillHist (const int Sample, TH1 *Hist, std::vector< double > &Array)
	Fill a histogram with the event-level information used in the fit. More...
void	FillArray_MP ()
	DB Nice new multi-threaded function which calculates the event weights and fills the relevant bins of an array. More...
void	FillArray ()
	Function which does the core reweighting, fills the SampleHandlerBase::SampleHandler_array vector with the weight calculated from reweighting. More...
void	ResetHistograms ()
	Helper function to reset histograms. 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 SampleHandlerInterface
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< SplineBase >	SplineHandler
	Contains all your splines (binned or unbinned) and handles the setup and the returning of weights from spline evaluations. More...
std::shared_ptr< OscillationHandler >	Oscillator
	Contains oscillator handling calculating oscillation probabilities. More...
std::vector< std::vector< FunctionalShifter * > >	funcParsGrid
	HH - a grid of vectors of enums for each sample and event. More...
std::vector< FunctionalShifter >	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::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::unordered_map< int, FuncParFuncType >	funcParsFuncMap
	HH - a map that relates the funcpar enum to pointer of the actual function. More...
std::vector< std::string >	funcParsNamesVec = {}
	HH - a vector of string names for each functional parameter. More...
std::unique_ptr< BinningHandler >	Binning
	KS: This stores binning information, in future could be come vector to store binning for every used sample. More...
std::vector< double >	SampleHandler_array
	DB Array to be filled after reweighting. More...
std::vector< double >	SampleHandler_array_w2
	KS Array used for MC stat. More...
std::vector< double >	SampleHandler_data
	DB Array to be filled in AddData. More...
std::vector< EventInfo >	MCEvents
	Stores information about every MC event. More...
std::vector< SampleInfo >	SampleDetails
	Stores info about currently initialised sample. More...
ParameterHandlerGeneric *	ParHandler = nullptr
	ETA - All experiments will need an xsec, det and osc cov. More...
std::string	SampleHandlerName
	A unique ID for each sample based on which we can define what systematic should be applied. More...
std::vector< 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< 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::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 SampleHandlerInterface
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 Member Functions
std::vector< std::vector< KinematicCut > >	ApplyTemporarySelection (const int iSample, const std::vector< KinematicCut > &ExtraCuts)
	Temporarily extend Selection for a given sample with additional cuts. Returns the original Selection so the caller can restore it later. More...

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 21 of file SampleHandlerBase.h.

Member Enumeration Documentation

FDPlotType

enum SampleHandlerBase::FDPlotType

private

Enumerator
kModePlot
kOscChannelPlot

Definition at line 426 of file SampleHandlerBase.h.

                  {
    kModePlot = 0,
    kOscChannelPlot = 1
  };

Constructor & Destructor Documentation

SampleHandlerBase()

SampleHandlerBase::SampleHandlerBase	(	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 12 of file SampleHandlerBase.cpp.

                                                                                              : SampleHandlerInterface() {
// ************************************************
  MACH3LOG_INFO("-------------------------------------------------------------------");
  MACH3LOG_INFO("Creating SampleHandlerBase object");
 
  //ETA - safety feature so you can't pass a NULL xsec_cov
  if(!xsec_cov) {
    MACH3LOG_WARN("You've passed me a nullptr ParameterHandler so I will not use any xsec parameters");
  }
  ParHandler = xsec_cov;
 
  nSamples = 1;
 
  if (OscillatorObj_ != nullptr) {
    MACH3LOG_WARN("You have passed an Oscillator object through the constructor of a SampleHandlerBase object - this will be used for all oscillation channels");
    Oscillator = OscillatorObj_;
    if(!ParHandler) {
      MACH3LOG_CRITICAL("You've passed me a nullptr to ParamHandler while non null to Oscillator");
      MACH3LOG_CRITICAL("Make up you mind");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
 
  KinematicParameters = nullptr;
  ReversedKinematicParameters = nullptr;
  KinematicVectors = nullptr;
  ReversedKinematicVectors = nullptr;
 
  SampleHandlerName = "";
  SampleManager = std::make_unique<Manager>(ConfigFileName.c_str());
  Binning = std::make_unique<BinningHandler>();
  // Variables related to MC stat
  FirstTimeW2 = true;
  UpdateW2 = false;
}

~SampleHandlerBase()

SampleHandlerBase::~SampleHandlerBase ( )

virtual

destructor

Definition at line 49 of file SampleHandlerBase.cpp.

                                      {
  MACH3LOG_DEBUG("I'm deleting SampleHandlerBase");
  if(THStackLeg != nullptr) delete THStackLeg;
}

Member Function Documentation

AddAdditionalWeightPointers()

virtual void SampleHandlerBase::AddAdditionalWeightPointers ( )

protectedpure virtual

DB Function to determine which weights apply to which types of samples.

Implemented in PySampleHandlerBase.

AddData() [1/2]

void SampleHandlerBase::AddData	(	const int	Sample,
		const std::vector< double > &	Data_Array
	)

Definition at line 1020 of file SampleHandlerBase.cpp.

                                                                                     {
// ************************************************
  const int Start = Binning->GetSampleStartBin(Sample);
  const int End = Binning->GetSampleEndBin(Sample);
  const int ExpectedSize = End - Start;
 
  if (static_cast<int>(Data_Array.size()) != ExpectedSize) {
    MACH3LOG_ERROR("Data_Array size mismatch for sample '{}'.", GetSampleTitle(Sample));
    MACH3LOG_ERROR("Expected size: {}, received size: {}.", ExpectedSize, Data_Array.size());
    MACH3LOG_ERROR("Start bin: {}, End bin: {}.", Start, End);
    MACH3LOG_ERROR("This likely indicates a binning or sample slicing bug.");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  std::copy_n(Data_Array.begin(), End-Start, SampleHandler_data.begin() + Start);
 
  FillHist(Sample, SampleDetails[Sample].DataHist, SampleHandler_data);
}

AddData() [2/2]

void SampleHandlerBase::AddData	(	const int	Sample,
		TH1 *	Data
	)

Definition at line 948 of file SampleHandlerBase.cpp.

                                                           {
// ************************************************
  int Dim = GetNDim(Sample);
  MACH3LOG_INFO("Adding {}D data histogram: {} with {:.2f} events", Dim, Data->GetTitle(), Data->Integral());
  // delete old histogram
  delete SampleDetails[Sample].DataHist;
  SampleDetails[Sample].DataHist = static_cast<TH1*>(Data->Clone());
 
  if(!SampleHandler_data.size()) {
    MACH3LOG_ERROR("SampleHandler_data haven't been initialised yet");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  auto ChecHistType = [&](const std::string& Type, const int Dimen, const TH1* Hist,
                          const std::string& file, const int line) {
    if (std::string(Hist->ClassName()) != Type) {
      MACH3LOG_ERROR("Expected {} for {}D sample, got {}", Type, Dimen, Hist->ClassName());
      throw MaCh3Exception(file, line);
    }
  };
 
  if (Dim == 1) {
    // Ensure we really have a TH1D
    ChecHistType("TH1D", Dim, SampleDetails[Sample].DataHist, __FILE__, __LINE__);
    M3::CheckBinningMatch(static_cast<TH1D*>(SampleDetails[Sample].DataHist),
                          static_cast<TH1D*>(SampleDetails[Sample].MCHist), __FILE__, __LINE__);
    for (int xBin = 0; xBin < Binning->GetNAxisBins(Sample, 0); ++xBin) {
      const int idx = Binning->GetGlobalBinSafe(Sample, {xBin});
      // ROOT histograms are 1-based, so bin index + 1
      SampleHandler_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(xBin + 1);
    }
    SampleDetails[Sample].DataHist->GetXaxis()->SetTitle(GetKinVarName(Sample, 0).c_str());
    SampleDetails[Sample].DataHist->GetYaxis()->SetTitle("Number of Events");
  } else if (Dim == 2) {
    if(Binning->IsUniform(Sample)) {
      ChecHistType("TH2D", Dim, SampleDetails[Sample].DataHist, __FILE__, __LINE__);
      M3::CheckBinningMatch(static_cast<TH2D*>(SampleDetails[Sample].DataHist),
                            static_cast<TH2D*>(SampleDetails[Sample].MCHist), __FILE__, __LINE__);
      for (int yBin = 0; yBin < Binning->GetNAxisBins(Sample, 1); ++yBin) {
        for (int xBin = 0; xBin < Binning->GetNAxisBins(Sample, 0); ++xBin) {
          const int idx = Binning->GetGlobalBinSafe(Sample, {xBin, yBin});
          //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
          SampleHandler_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(xBin + 1, yBin + 1);
        }
      }
    } else{
      ChecHistType("TH2Poly", Dim, SampleDetails[Sample].DataHist, __FILE__, __LINE__);
      M3::CheckBinningMatch(static_cast<TH2Poly*>(SampleDetails[Sample].DataHist),
                            static_cast<TH2Poly*>(SampleDetails[Sample].MCHist), __FILE__, __LINE__);
      for (int iBin = 0; iBin < Binning->GetNBins(Sample); ++iBin) {
        const int idx = iBin + Binning->GetSampleStartBin(Sample);
        //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
        SampleHandler_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(iBin + 1);
      }
    }
 
    SampleDetails[Sample].DataHist->GetXaxis()->SetTitle(GetKinVarName(Sample, 0).c_str());
    SampleDetails[Sample].DataHist->GetYaxis()->SetTitle(GetKinVarName(Sample, 1).c_str());
    SampleDetails[Sample].DataHist->GetZaxis()->SetTitle("Number of Events");
  } else {
    ChecHistType("TH1D", Dim, SampleDetails[Sample].DataHist, __FILE__, __LINE__);
    M3::CheckBinningMatch(static_cast<TH1D*>(SampleDetails[Sample].DataHist),
                          static_cast<TH1D*>(SampleDetails[Sample].MCHist), __FILE__, __LINE__);
    for (int iBin = 0; iBin < Binning->GetNBins(Sample); ++iBin) {
      const int idx = iBin + Binning->GetSampleStartBin(Sample);
      //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
      SampleHandler_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(iBin + 1);
    }
  }
}

ApplyShifts()

void SampleHandlerBase::ApplyShifts ( const int  iEvent )

protectedvirtual

ETA - generic function applying shifts.

Definition at line 549 of file SampleHandlerBase.cpp.

                                                    {
// ***************************************************************************
  const auto& shifts = funcParsGrid[iEvent];
  const auto nShifts = shifts.size();
  // KS: If there are no shifts then there is no point in resetting which can be costly.
  if(nShifts == 0) {
    return;
  }
 
  // 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);
 
  for (std::size_t iShift = 0; iShift < nShifts; ++iShift) {
    const auto* _restrict_ fp = shifts[iShift];
    (*fp->funcPtr)(fp->valuePtr, iEvent);
  }
 
  FinaliseShifts(iEvent);
}

ApplyTemporarySelection()

std::vector< std::vector< KinematicCut > > SampleHandlerBase::ApplyTemporarySelection ( const int  iSample, const std::vector< KinematicCut > &  ExtraCuts  )

private

Temporarily extend Selection for a given sample with additional cuts. Returns the original Selection so the caller can restore it later.

Definition at line 1404 of file SampleHandlerBase.cpp.

                                                                                   {
// ************************************************
  // Backup current selection
  auto originalSelection = Selection;
 
  //DB Add all the predefined selections to the selection vector which will be applied
  auto selectionToApply = Selection;
 
  //DB Add all requested cuts from the argument to the selection vector which will be applied
  for (const auto& cut : ExtraCuts) {
    selectionToApply[iSample].emplace_back(cut);
  }
 
  //DB Set the member variable to be the cuts to apply
  Selection = std::move(selectionToApply);
 
  return originalSelection;
}

BuildModeChannelSelection()

std::vector< KinematicCut > SampleHandlerBase::BuildModeChannelSelection	(	const int	iSample,
		const int	kModeToFill,
		const int	kChannelToFill
	)		const

Definition at line 1734 of file SampleHandlerBase.cpp.

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

CalcNormsBins()

void SampleHandlerBase::CalcNormsBins ( std::vector< NormParameter > &  norm_parameters, std::vector< std::vector< int > > &  norms_bins  )

protected

Check whether a normalisation systematic affects an event or not.

Definition at line 658 of file SampleHandlerBase.cpp.

                                                                                                                          {
// ************************************************
  #ifdef MACH3_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 (MCEvents[iEvent].isNC && (MCEvents[iEvent].nupdg != MCEvents[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
        const int Target = static_cast<int>(std::round(ReturnKinematicParameter("TargetNucleus", iEvent)));
        bool TargetMatch = MatchCondition(it->targets, Target);
        if (!TargetMatch) {
          MACH3LOG_TRACE("Event {}, missed target check ({}) for dial {}", iEvent, Target, it->name);
          continue;
        }
 
        //Now check that the neutrino flavour in an interaction matches with the normalisation parameters
        bool FlavourMatch = MatchCondition(it->pdgs, MCEvents[iEvent].nupdg);
        if (!FlavourMatch) {
          MACH3LOG_TRACE("Event {}, missed PDG check ({}) for dial {}", iEvent,MCEvents[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, MCEvents[iEvent].nupdgUnosc);
        if (!FlavourUnoscMatch){
          MACH3LOG_TRACE("Event {}, missed FlavourUnosc check ({}) for dial {}", iEvent,MCEvents[iEvent].nupdgUnosc, it->name);
          continue;
        }
 
        //Now check that the mode of an interaction matches with the normalisation parameters
        const int Mode = static_cast<int>(std::round(ReturnKinematicParameter("Mode", iEvent)));
        bool ModeMatch = MatchCondition(it->modes, Mode);
        if (!ModeMatch) {
          MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", 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 = PassesSelection((*it), iEvent);
        // 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 MACH3_DEBUG
        VerboseCounter[std::distance(norm_parameters.begin(), it)]++;
        #endif
        //}
      } // end iteration over norm_parameters
    } // end if (ParHandler)
    norms_bins[iEvent] = NormBins;
  }//end loop over events
  #ifdef MACH3_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 SampleHandlerBase::CalcWeightFunc ( const 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 294 of file SampleHandlerBase.h.

{return; (void)iEvent;};

CalcWeightTotal()

M3::float_t SampleHandlerBase::CalcWeightTotal ( const EventInfo *_restrict_  MCEvent ) const

protected

Calculate the total weight weight for a given event.

Definition at line 572 of file SampleHandlerBase.cpp.

                                                                                                 {
// ***************************************************************************
  M3::float_t TotalWeight = 1.0;
  const int TotalWeights = static_cast<int>(MCEvent->total_weight_pointers.size());
  //DB Loop over stored pointers
  #ifdef MULTITHREAD
  #pragma omp simd reduction(*:TotalWeight)
  #endif
  for (int iWeight = 0; iWeight < TotalWeights; ++iWeight) {
    TotalWeight *= *(MCEvent->total_weight_pointers[iWeight]);
  }
 
  return TotalWeight;
}

Fill1DSubEventHist()

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

Definition at line 1469 of file SampleHandlerBase.cpp.

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

Fill2DSubEventHist()

void SampleHandlerBase::Fill2DSubEventHist	(	const int	iSample,
		TH2 *	_h2DVar,
		const std::string &	ProjectionVarX,
		const std::string &	ProjectionVarY,
		const std::vector< KinematicCut > &	SubEventSelectionVec = {},
		int	WeightStyle = 0
	)

Definition at line 1552 of file SampleHandlerBase.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++) {
    const int EventSample = MCEvents[iEvent].NominalSample;
    if(EventSample != iSample) continue;
    if (IsEventSelected(EventSample, 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 SampleHandlerBase::FillArray ( )

protected

Function which does the core reweighting, fills the SampleHandlerBase::SampleHandler_array vector with the weight calculated from reweighting.

Definition at line 355 of file SampleHandlerBase.cpp.

                                  {
//************************************************
  //DB Reset which cuts to apply
  Selection = StoredSelection;
 
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    ApplyShifts(iEvent);
    const EventInfo* _restrict_ MCEvent = &MCEvents[iEvent];
 
    if (!IsEventSelected(MCEvent->NominalSample, iEvent)) {
      continue;
    }
 
    // Virtual by default does nothing, has to happen before CalcWeightTotal
    CalcWeightFunc(iEvent);
 
    const M3::float_t totalweight = CalcWeightTotal(MCEvent);
    //DB Catch negative total weights and skip any event with a negative weight. Previously we would set weight to zero and continue but that is inefficient
    if (totalweight <= 0.){
      continue;
    }
 
    //DB Find the relevant bin in the PDF for each event
    const int GlobalBin = Binning->FindGlobalBin(MCEvent->NominalSample, MCEvent->KinVar, MCEvent->NomBin);
 
    //DB Fill relevant part of thread array
    if (GlobalBin > M3::UnderOverFlowBin) {
      SampleHandler_array[GlobalBin] += totalweight;
      if (FirstTimeW2) SampleHandler_array_w2[GlobalBin] += totalweight*totalweight;
    }
  }
}

FillArray_MP()

void SampleHandlerBase::FillArray_MP ( )

protected

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

Multithreaded version of fillArray.

Function which does the core reweighting, fills the SampleHandlerBase::SampleHandler_array vector with the weight calculated from reweighting but multithreaded

See also

fillArray()

Definition at line 393 of file SampleHandlerBase.cpp.

                                     {
// ************************************************
  //DB Reset which cuts to apply
  Selection = StoredSelection;
 
  // NOTE comment below is left for historical reasons
  //DB - Brain dump of speedup ideas
  //
  //Those relevant to reweighting
  // 1. Don't bother storing and calculating NC signal events - Implemented and saves marginal s/step
  // 2. Loop over spline event weight calculation in the following event loop - Currently done in splineSKBase->calcWeight() where multi-threading won't be optimised - Implemented and saves 0.3s/step
  // 3. Inline getDiscVar or somehow include that calculation inside the multi-threading - Implemented and saves about 0.01s/step
  // 4. Include isCC inside SKMCStruct so don't have to have several 'if' statements determine if oscillation weight needs to be set to 1.0 for NC events - Implemented and saves marginal s/step
  // 5. Do explicit check on adjacent bins when finding event XBin instead of looping over all BinEdge indices - Implemented but doesn't significantly affect s/step
  //
  //Other aspects
  // 1. Order minituples in Y-axis variable as this will *hopefully* reduce cache misses inside SampleHandler_array_class[yBin][xBin]
  //
  // We will hit <0.1 s/step eventually! :D
  const auto TotalBins = Binning->GetNBins();
  const unsigned int NumberOfEvents = GetNEvents();
 
  double* _restrict_ MC_Array_for_reduction = SampleHandler_array.data();
  double* _restrict_ W2_array_for_reduction = SampleHandler_array_w2.data();
 
  #pragma omp parallel for reduction(+:MC_Array_for_reduction[:TotalBins], W2_array_for_reduction[:TotalBins])
  for (unsigned int iEvent = 0; iEvent < NumberOfEvents; ++iEvent) {
    //ETA - generic functions to apply shifts to kinematic variables
    // Apply this before IsEventSelected is called.
    ApplyShifts(iEvent);
 
    const EventInfo* _restrict_ MCEvent = &MCEvents[iEvent];
    //ETA - generic functions to apply shifts to kinematic variable
    if(!IsEventSelected(MCEvent->NominalSample, iEvent)){
      continue;
    }
 
    // Virtual by default does nothing, has to happen before CalcWeightTotal
    CalcWeightFunc(iEvent);
 
    const M3::float_t totalweight = CalcWeightTotal(MCEvent);
    //DB Catch negative total weights and skip any event with a negative weight. Previously we would set weight to zero and continue but that is inefficient
    if (totalweight <= 0.){
      continue;
    }
 
    //DB Find the relevant bin in the PDF for each event
    const int GlobalBin = Binning->FindGlobalBin(MCEvent->NominalSample, MCEvent->KinVar, MCEvent->NomBin);
 
    //ETA - we can probably remove this final if check on the -1?
    //Maybe we can add an overflow bin to the array and assign any events to this bin?
    //Might save us an extra if call?
    //DB Fill relevant part of thread array
    if (GlobalBin > M3::UnderOverFlowBin) {
      MC_Array_for_reduction[GlobalBin] += totalweight;
      if (FirstTimeW2) W2_array_for_reduction[GlobalBin] += totalweight*totalweight;
    }
  }
}

FillHist()

void SampleHandlerBase::FillHist ( const int  Sample, TH1 *  Hist, std::vector< double > &  Array  )

protected

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

DB Functions required for reweighting functions DB Replace previous implementation with reading bin contents from SampleHandler_array

Definition at line 245 of file SampleHandlerBase.cpp.

                                                                                      {
// ************************************************
  int Dimension = GetNDim(Sample);
  // 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) {
    Hist->Reset();
    for (int xBin = 0; xBin < Binning->GetNAxisBins(Sample, 0); ++xBin) {
      const int idx = Binning->GetGlobalBinSafe(Sample, {xBin});
      Hist->SetBinContent(xBin + 1, Array[idx]);
    }
  } else if (Dimension == 2) {
    Hist->Reset();
    if(Binning->IsUniform(Sample)) {
      for (int yBin = 0; yBin < Binning->GetNAxisBins(Sample, 1); ++yBin) {
        for (int xBin = 0; xBin < Binning->GetNAxisBins(Sample, 0); ++xBin) {
          const int idx = Binning->GetGlobalBinSafe(Sample, {xBin, yBin});
          Hist->SetBinContent(xBin + 1, yBin + 1, Array[idx]);
        }
      }
    } else {
      for (int iBin = 0; iBin < Binning->GetNBins(Sample); ++iBin) {
        const int idx = iBin + Binning->GetSampleStartBin(Sample);
        //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
        Hist->SetBinContent(iBin + 1, Array[idx]);
      }
    }
  } else {
    for (int iBin = 0; iBin < Binning->GetNBins(Sample); ++iBin) {
      const int idx = iBin + Binning->GetSampleStartBin(Sample);
      //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
      Hist->SetBinContent(iBin + 1, Array[idx]);
    }
  }
}

FinaliseShifts()

virtual void SampleHandlerBase::FinaliseShifts ( const int  iEvent )

inlineprotectedvirtual

LP - Optionally calculate derived observables after all shifts have been applied.

LP - For example, have shifts that varied lepton energy and hadron energy separately in a subclass implementation of this method you may add the shifted quantities together to build a shifted neutrino energy estimator

Definition at line 260 of file SampleHandlerBase.h.

{(void)iEvent;};

FindNominalBinAndEdges()

void SampleHandlerBase::FindNominalBinAndEdges ( )

protected

Definition at line 851 of file SampleHandlerBase.cpp.

                                               {
// ************************************************
  for (unsigned int event_i = 0; event_i < GetNEvents(); event_i++) {
    int Sample = MCEvents[event_i].NominalSample;
    const int dim = GetNDim(Sample);
    MCEvents[event_i].KinVar.resize(dim);
    MCEvents[event_i].NomBin.resize(dim);
 
    auto SetNominalBin = [&](int bin, int max_bins, int& out_bin) {
      if (bin >= 0 && bin < max_bins) {
        out_bin = bin;
      } else {
        out_bin = M3::UnderOverFlowBin;  // Out of bounds
      }
    };
 
    // Find nominal bin for each dimension
    for(int iDim = 0; iDim < dim; iDim++) {
      MCEvents[event_i].KinVar[iDim] = GetPointerToKinematicParameter(GetKinVarName(Sample, iDim), event_i);
      if (std::isnan(*MCEvents[event_i].KinVar[iDim]) || std::isinf(*MCEvents[event_i].KinVar[iDim])) {
        MACH3LOG_ERROR("Variable {} for sample {} and dimension {} is ill-defined and equal to {}",
                       GetKinVarName(Sample, iDim), GetSampleTitle(Sample), dim, *MCEvents[event_i].KinVar[iDim]);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      const int bin = Binning->FindNominalBin(Sample, iDim, *MCEvents[event_i].KinVar[iDim]);
      int NBins_i = static_cast<int>(Binning->GetBinEdges(Sample, iDim).size() - 1);
      SetNominalBin(bin, NBins_i, MCEvents[event_i].NomBin[iDim]);
    }
  }
}

Get1DVarHist()

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

finalvirtual

Implements SampleHandlerInterface.

Definition at line 1426 of file SampleHandlerBase.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
  auto tmp_Selection = ApplyTemporarySelection(iSample, EventSelectionVec);
 
  //DB Define the histogram which will be returned
  std::vector<double> xBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_Str);
  auto _h1DVar = std::make_unique<TH1D>("", "", int(xBinEdges.size())-1, xBinEdges.data());
  _h1DVar->SetDirectory(nullptr);
  _h1DVar->GetXaxis()->SetTitle(ProjectionVar_Str.c_str());
  _h1DVar->GetYaxis()->SetTitle("Events");
 
  if (IsSubEventVarString(ProjectionVar_Str)) {
    Fill1DSubEventHist(iSample, _h1DVar.get(), ProjectionVar_Str, SubEventSelectionVec, WeightStyle);
  } else {
    //DB Grab the associated enum with the argument string
    int ProjectionVar_Int = ReturnKinematicParameterFromString(ProjectionVar_Str);
 
    //DB Loop over all events
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
      const int EventSample = MCEvents[iEvent].NominalSample;
      if(EventSample != iSample) continue;
      if (IsEventSelected(EventSample, iEvent)) {
        double Weight = GetEventWeight(iEvent);
        if (WeightStyle == 1) {
          Weight = 1.;
        }
        double Var = ReturnKinematicParameter(ProjectionVar_Int, iEvent);
        _h1DVar->Fill(Var, Weight);
      }
    }
  }
  //DB Reset the saved selection
  Selection = tmp_Selection;
 
  return _h1DVar;
}

Get1DVarHistByModeAndChannel()

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

finalvirtual

Implements SampleHandlerInterface.

Definition at line 1785 of file SampleHandlerBase.cpp.

                                                                            {
// ************************************************
  auto SelectionVec = BuildModeChannelSelection(iSample, kModeToFill, kChannelToFill);
  return Get1DVarHist(iSample, ProjectionVar_Str, SelectionVec, WeightStyle);
}

Get2DVarHist()

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

finalvirtual

Implements SampleHandlerInterface.

Definition at line 1497 of file SampleHandlerBase.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
  auto tmp_Selection = ApplyTemporarySelection(iSample, EventSelectionVec);
 
  //DB Define the histogram which will be returned
  //KS: If we use 2D non uniform binning and wanting to plot fit variables use TH2Poly everything else uses TH2D with binning from config
  std::unique_ptr<TH2> _h2DVar;
  if(GetNDim(iSample) == 2 && !Binning->IsUniform(iSample) &&
     ProjectionVar_StrX == GetKinVarName(iSample, 0) && ProjectionVar_StrY == GetKinVarName(iSample, 1) ) {
    _h2DVar = std::unique_ptr<TH2>(static_cast<TH2*>(M3::Clone(GetMCHist(iSample)).release()));
  } else {
    std::vector<double> xBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_StrX);
    std::vector<double> yBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_StrY);
    _h2DVar = std::make_unique<TH2D>("", "", int(xBinEdges.size())-1, xBinEdges.data(), int(yBinEdges.size())-1, yBinEdges.data());
  }
  _h2DVar->SetDirectory(nullptr);
  _h2DVar->GetXaxis()->SetTitle(ProjectionVar_StrX.c_str());
  _h2DVar->GetYaxis()->SetTitle(ProjectionVar_StrY.c_str());
  _h2DVar->GetZaxis()->SetTitle("Events");
 
  bool IsSubEventHist = IsSubEventVarString(ProjectionVar_StrX) || IsSubEventVarString(ProjectionVar_StrY);
  if (IsSubEventHist) Fill2DSubEventHist(iSample, _h2DVar.get(), 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++) {
      const int EventSample = MCEvents[iEvent].NominalSample;
      if(EventSample != iSample) continue;
      if (IsEventSelected(EventSample, 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;
}

Get2DVarHistByModeAndChannel()

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

finalvirtual

Implements SampleHandlerInterface.

Definition at line 1793 of file SampleHandlerBase.cpp.

                                                                                                                      {
// ************************************************
  auto SelectionVec = BuildModeChannelSelection(iSample, kModeToFill, kChannelToFill);
  return Get2DVarHist(iSample, ProjectionVar_StrX, ProjectionVar_StrY, SelectionVec, WeightStyle);
}

GetArrayForSample()

std::vector< double > SampleHandlerBase::GetArrayForSample	(	const int	Sample,
		std::vector< double > const &	array
	)		const

Return a sub-array for a given sample.

Definition at line 2098 of file SampleHandlerBase.cpp.

                                                                                                              {
// ***************************************************************************
  const int Start = Binning->GetSampleStartBin(Sample);
  const int End   = Binning->GetSampleEndBin(Sample);
 
  return std::vector<double>(array.begin() + Start, array.begin() + End);
}

GetDataArray() [1/2]

auto SampleHandlerBase::GetDataArray ( ) const

inline

Return array storing data entries for every bin.

Definition at line 153 of file SampleHandlerBase.h.

                            {
    return SampleHandler_data;
  }

GetDataArray() [2/2]

std::vector<double> SampleHandlerBase::GetDataArray ( const int  Sample ) const

inline

Return array storing data entries for every bin.

Definition at line 168 of file SampleHandlerBase.h.

                                                         {
    return GetArrayForSample(Sample, SampleHandler_data);
  }

GetDataHist() [1/2]

const TH1 * SampleHandlerBase::GetDataHist ( const int  Sample )

finalvirtual

Get Data histogram.

Implements SampleHandlerInterface.

Definition at line 931 of file SampleHandlerBase.cpp.

                                                          {
// ************************************************
  if(SampleDetails[Sample].DataHist == nullptr) {
    MACH3LOG_ERROR("Can't access {} for {}Dimensions", __func__, GetNDim(Sample));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  return SampleDetails[Sample].DataHist;
}

GetDataHist() [2/2]

const TH1 * SampleHandlerBase::GetDataHist

(

const std::string &

Sample

)

Definition at line 941 of file SampleHandlerBase.cpp.

                                                                 {
// ************************************************
  int Index = GetSampleIndex(Sample);
  return GetDataHist(Index);
}

GetEventWeight()

M3::float_t SampleHandlerBase::GetEventWeight

(

const int

iEntry

)

Computes the total event weight for a given entry.

KS:

Warning

we have to here recalculate weight and cap because there is possibility weight wasn't calculated during FillArray because it didn't fulfil IsEventSelected

Definition at line 1176 of file SampleHandlerBase.cpp.

                                                            {
// ************************************************
 
  // Virtual by default does nothing, has to happen before CalcWeightTotal
  CalcWeightFunc(iEntry);
 
  const EventInfo* _restrict_ MCEvent = &MCEvents[iEntry];
  M3::float_t totalweight = CalcWeightTotal(MCEvent);
 
  //DB Catch negative total weights and skip any event with a negative weight. Previously we would set weight to zero and continue but that is inefficient
  if (totalweight <= 0.){
    totalweight = 0.;
  }
  return totalweight;
}

GetFinalPDGFromFileName()

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

inlineprotected

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

Definition at line 415 of file SampleHandlerBase.h.

{return FileToFinalPDGMap.at(FileName);}

GetFlavourName()

std::string SampleHandlerBase::GetFlavourName ( const int  iSample, const int  iChannel  ) const

inlinefinalvirtual

Implements SampleHandlerInterface.

Definition at line 90 of file SampleHandlerBase.h.

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

GetInitPDGFromFileName()

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

inlineprotected

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

Definition at line 413 of file SampleHandlerBase.h.

{return FileToInitPDGMap.at(FileName);}

GetKinVarName()

std::string SampleHandlerBase::GetKinVarName ( const int  iSample, const int  Dimension  ) const

finalvirtual

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

Parameters

iSample	Sample index
Dimension	Dimension index

Implements SampleHandlerInterface.

Definition at line 2088 of file SampleHandlerBase.cpp.

                                                                                       {
// ***************************************************************************
  if(Dimension > GetNDim(iSample)) {
    MACH3LOG_ERROR("Asking for dimension {}, while sample: {} only has {}", Dimension, GetSampleTitle(iSample), GetNDim(iSample));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  return SampleDetails[iSample].VarStr[Dimension];
}

GetLikelihood()

double SampleHandlerBase::GetLikelihood ( ) const

overridevirtual

DB Multi-threaded GetLikelihood.

Implements SampleHandlerInterface.

Definition at line 1291 of file SampleHandlerBase.cpp.

                                              {
// ************************************************
  double negLogL = 0.;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:negLogL)
  #endif
  for (int idx = 0; idx < Binning->GetNBins(); ++idx)
  {
    double const &DataVal = SampleHandler_data[idx];
    double const &MCPred = SampleHandler_array[idx];
    double const &w2 = SampleHandler_array_w2[idx];
 
    //KS: Calculate likelihood using Barlow-Beeston Poisson or even IceCube
    negLogL += GetTestStatLLH(DataVal, MCPred, w2);
  }
  return negLogL;
}

GetMCArray() [1/2]

auto SampleHandlerBase::GetMCArray ( ) const

inline

Return array storing MC entries for every bin.

Definition at line 157 of file SampleHandlerBase.h.

                          {
    return SampleHandler_array;
  }

GetMCArray() [2/2]

std::vector<double> SampleHandlerBase::GetMCArray ( const int  Sample ) const

inline

Return array storing MC entries for every bin.

Definition at line 172 of file SampleHandlerBase.h.

                                                       {
    return GetArrayForSample(Sample, SampleHandler_array);
  }

GetMCHist() [1/2]

const TH1 * SampleHandlerBase::GetMCHist ( const int  Sample )

finalvirtual

Get MC histogram.

Implements SampleHandlerInterface.

Definition at line 913 of file SampleHandlerBase.cpp.

                                                        {
// ************************************************
  FillHist(Sample, SampleDetails[Sample].MCHist, SampleHandler_array);
  if(SampleDetails[Sample].MCHist == nullptr) {
    MACH3LOG_ERROR("Can't access {} for {}Dimensions", __func__, GetNDim(Sample));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  return SampleDetails[Sample].MCHist;
}

GetMCHist() [2/2]

const TH1 * SampleHandlerBase::GetMCHist

(

const std::string &

Sample

)

Definition at line 924 of file SampleHandlerBase.cpp.

                                                               {
// ************************************************
  const int Index = GetSampleIndex(Sample);
  return GetMCHist(Index);
}

GetModeHist1D()

std::unique_ptr<TH1> SampleHandlerBase::GetModeHist1D ( const int  iSample, int  s, int  m, int  style = 0  )

inline

Definition at line 118 of file SampleHandlerBase.h.

                                                                                   {
    return Get1DVarHistByModeAndChannel(iSample, GetKinVarName(iSample, 0), m, s, style);
  }

GetModeHist2D()

std::unique_ptr<TH2> SampleHandlerBase::GetModeHist2D ( const int  iSample, int  s, int  m, int  style = 0  )

inline

Definition at line 121 of file SampleHandlerBase.h.

                                                                                   {
    return Get2DVarHistByModeAndChannel(iSample, GetKinVarName(iSample, 0), GetKinVarName(iSample, 1), m, s, style);
  }

GetName()

std::string SampleHandlerBase::GetName ( ) const

finalvirtual

Get name for Sample Handler.

Implements SampleHandlerInterface.

Definition at line 1163 of file SampleHandlerBase.cpp.

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

GetNDim()

int SampleHandlerBase::GetNDim ( const int  Sample ) const

inlinefinalvirtual

DB Get what dimensionality binning for given sample has.

Parameters

Sample

Number of sample

Implements SampleHandlerInterface.

Definition at line 34 of file SampleHandlerBase.h.

{ return SampleDetails[Sample].nDimensions; }

GetNOscChannels()

int SampleHandlerBase::GetNOscChannels ( const int  iSample ) const

inlinefinalvirtual

Get number of oscillation channels for a single sample.

Implements SampleHandlerInterface.

Definition at line 88 of file SampleHandlerBase.h.

{return static_cast<int>(SampleDetails[iSample].OscChannels.size());};

GetNuOscillatorPointers()

const M3::float_t * SampleHandlerBase::GetNuOscillatorPointers

(

const int

iEvent

)

const

Definition at line 1123 of file SampleHandlerBase.cpp.

                                                                                  {
// ************************************************
  const M3::float_t* osc_w_pointer = &M3::Unity;
  if (MCEvents[iEvent].isNC) {
    if (MCEvents[iEvent].nupdg != MCEvents[iEvent].nupdgUnosc) {
      osc_w_pointer = &M3::Zero;
    } else {
      osc_w_pointer = &M3::Unity;
    }
  } else {
    int InitFlav = M3::_BAD_INT_;
    int FinalFlav = M3::_BAD_INT_;
 
    InitFlav =  M3::Utils::PDGToNuOscillatorFlavour(MCEvents[iEvent].nupdgUnosc);
    FinalFlav = M3::Utils::PDGToNuOscillatorFlavour(MCEvents[iEvent].nupdg);
 
    if (InitFlav == M3::_BAD_INT_ || FinalFlav == M3::_BAD_INT_) {
      MACH3LOG_ERROR("Something has gone wrong in the mapping between MCEvents.nutype and the enum used within NuOscillator");
      MACH3LOG_ERROR("MCEvents.nupdgUnosc: {}", MCEvents[iEvent].nupdgUnosc);
      MACH3LOG_ERROR("InitFlav: {}", InitFlav);
      MACH3LOG_ERROR("MCEvents.nupdg: {}", MCEvents[iEvent].nupdg);
      MACH3LOG_ERROR("FinalFlav: {}", FinalFlav);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    const int Sample = MCEvents[iEvent].NominalSample;
 
    const int OscIndex = GetOscChannel(SampleDetails[Sample].OscChannels, MCEvents[iEvent].nupdgUnosc, MCEvents[iEvent].nupdg);
    //Can only happen if truecz has been initialised within the experiment specific code
    if (MCEvents[iEvent].coszenith_true != M3::_BAD_DOUBLE_) {
      //Atmospherics
      osc_w_pointer = Oscillator->GetNuOscillatorPointers(Sample, OscIndex, InitFlav, FinalFlav, FLOAT_T(MCEvents[iEvent].enu_true), FLOAT_T(MCEvents[iEvent].coszenith_true));
    } else {
      //Beam
      osc_w_pointer = Oscillator->GetNuOscillatorPointers(Sample, OscIndex, InitFlav, FinalFlav, FLOAT_T(MCEvents[iEvent].enu_true));
    }
  } // end if NC
  return osc_w_pointer;
}

GetPointerToKinematicParameter() [1/2]

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

protectedpure virtual

Implemented in PySampleHandlerBase.

GetPointerToKinematicParameter() [2/2]

const double* SampleHandlerBase::GetPointerToKinematicParameter ( const std::string &  KinematicParameter, int  iEvent  ) const

inlineprotected

Definition at line 313 of file SampleHandlerBase.h.

                                                                                                      {
    return GetPointerToKinematicParameter(ReturnKinematicParameterFromString(KinematicParameter), iEvent);
  }

GetPointerToOscChannel()

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

protected

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

Definition at line 2027 of file SampleHandlerBase.cpp.

                                                                              {
// ************************************************
  auto& OscillationChannels = SampleDetails[MCEvents[iEvent].NominalSample].OscChannels;
  const int Channel = GetOscChannel(OscillationChannels, MCEvents[iEvent].nupdgUnosc, MCEvents[iEvent].nupdg);
  return &(OscillationChannels[Channel].ChannelIndex);
}

GetSampleIndex()

int SampleHandlerBase::GetSampleIndex

(

const std::string &

SampleTitle

)

const

Get index of sample based on name.

Parameters

SampleTitle

The title of the sample to search for.

Definition at line 883 of file SampleHandlerBase.cpp.

                                                                        {
// ************************************************
  for (M3::int_t iSample = 0; iSample < nSamples; ++iSample) {
    if (SampleTitle == GetSampleTitle(iSample)) {
      return iSample;
    }
  }
  MACH3LOG_ERROR("Sample name not found: {}", SampleTitle);
  throw MaCh3Exception(__FILE__, __LINE__);
}

GetSampleLikelihood()

double SampleHandlerBase::GetSampleLikelihood ( const int  isample ) const

overridevirtual

Get likelihood for single sample.

Implements SampleHandlerInterface.

Definition at line 1269 of file SampleHandlerBase.cpp.

                                                                     {
// ************************************************
  const int Start = Binning->GetSampleStartBin(isample);
  const int End = Binning->GetSampleEndBin(isample);
 
  double negLogL = 0.;
  #ifdef MULTITHREAD
  #pragma omp parallel for reduction(+:negLogL)
  #endif
  for (int idx = Start; idx < End; ++idx)
  {
    double const &DataVal = SampleHandler_data[idx];
    double const &MCPred = SampleHandler_array[idx];
    double const &w2 = SampleHandler_array_w2[idx];
 
    //KS: Calculate likelihood using Barlow-Beeston Poisson or even IceCube
    negLogL += GetTestStatLLH(DataVal, MCPred, w2);
  }
  return negLogL;
}

GetSampleTitle()

std::string SampleHandlerBase::GetSampleTitle ( const int  Sample ) const

inlinefinalvirtual

Get fancy title for specified samples.

Implements SampleHandlerInterface.

Definition at line 38 of file SampleHandlerBase.h.

{return SampleDetails[Sample].SampleTitle;}

GetSplineBins()

std::vector< std::vector< int > > SampleHandlerBase::GetSplineBins ( int  Event, BinnedSplineHandler *  BinnedSpline, bool &  ThrowCrititcal  ) const

protected

Retrieve the spline bin indices associated with a given event.

Warning

ThrowCrititcal argument will be eventually removed

Definition at line 1194 of file SampleHandlerBase.cpp.

                                                                                                                                     {
// ************************************************
  const int SampleIndex = MCEvents[Event].NominalSample;
  const auto SampleTitle = GetSampleTitle(SampleIndex);
  bool NoOscChannels = false;
  if(Oscillator == nullptr && GetNOscChannels(SampleIndex) == 1) {
    MACH3LOG_DEBUG("Assuming there are no osc channels in {}", __func__);
    NoOscChannels = true;
  }
  const int OscIndex = NoOscChannels ? 0 : GetOscChannel(SampleDetails[SampleIndex].OscChannels,
                                                         MCEvents[Event].nupdgUnosc, MCEvents[Event].nupdg);
  const int Mode = static_cast<int>(std::round(ReturnKinematicParameter("Mode", Event)));
  const double Etrue = MCEvents[Event].enu_true;
  std::vector< std::vector<int> > EventSplines;
  switch(GetNDim(SampleIndex)) {
    case 1:
      EventSplines = BinnedSpline->GetEventSplines(SampleTitle, OscIndex, Mode, Etrue, *(MCEvents[Event].KinVar[0]), 0.);
      break;
    case 2:
      EventSplines = BinnedSpline->GetEventSplines(SampleTitle, OscIndex, Mode, Etrue, *(MCEvents[Event].KinVar[0]), *(MCEvents[Event].KinVar[1]));
      break;
    default:
      if(ThrowCrititcal) {
        MACH3LOG_CRITICAL("MaCh3 doesn't support binned splines for more than 2D while you use {}", GetNDim(SampleIndex));
        MACH3LOG_CRITICAL("Will use 2D like approach");
        ThrowCrititcal = false;
      }
      EventSplines = BinnedSpline->GetEventSplines(SampleTitle, OscIndex, Mode, Etrue, *(MCEvents[Event].KinVar[0]), *(MCEvents[Event].KinVar[1]));
      break;
  }
  return EventSplines;
}

GetW2Array() [1/2]

auto SampleHandlerBase::GetW2Array ( ) const

inline

Return array storing W2 entries for every bin.

Definition at line 161 of file SampleHandlerBase.h.

                          {
    return SampleHandler_array_w2;
  }

GetW2Array() [2/2]

std::vector<double> SampleHandlerBase::GetW2Array ( const int  Sample ) const

inline

Return array storing W2 entries for single sample.

Definition at line 176 of file SampleHandlerBase.h.

                                                       {
    return GetArrayForSample(Sample, SampleHandler_array_w2);
  }

GetW2Hist() [1/2]

const TH1 * SampleHandlerBase::GetW2Hist ( const int  Sample )

finalvirtual

Get W2 histogram.

Implements SampleHandlerInterface.

Definition at line 895 of file SampleHandlerBase.cpp.

                                                        {
// ************************************************
  FillHist(Sample, SampleDetails[Sample].W2Hist, SampleHandler_array_w2);
  if(SampleDetails[Sample].W2Hist == nullptr) {
    MACH3LOG_ERROR("Can't access {} for {}Dimensions", __func__, GetNDim(Sample));
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  return SampleDetails[Sample].W2Hist;
}

GetW2Hist() [2/2]

const TH1 * SampleHandlerBase::GetW2Hist

(

const std::string &

Sample

)

Definition at line 906 of file SampleHandlerBase.cpp.

                                                               {
// ************************************************
  const int Index = GetSampleIndex(Sample);
  return GetW2Hist(Index);
}

Init()

virtual void SampleHandlerBase::Init ( )

protectedpure virtual

Initialise any variables that your experiment specific SampleHandler needs.

Implemented in PySampleHandlerBase.

Initialise()

void SampleHandlerBase::Initialise ( )

protected

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

Definition at line 170 of file SampleHandlerBase.cpp.

                                   {
// ************************************************
  TStopwatch clock;
  clock.Start();
 
  //First grab all the information from your sample config via your manager
  ReadConfig();
 
  //Now initialise all the variables you will need
  Init();
 
  nEvents = SetupExperimentMC();
  MCEvents.resize(nEvents);
  SetupMC();
 
  MACH3LOG_INFO("=============================================");
  MACH3LOG_INFO("Total number of events is: {}", GetNEvents());
  SetupOscParameters();
  MACH3LOG_INFO("Setting up Sample Binning..");
  SetBinning();
  MACH3LOG_INFO("Setting up Splines..");
  SetupSplines();
  MACH3LOG_INFO("Setting up Normalisation Pointers..");
  SetupNormParameters();
  MACH3LOG_INFO("Setting up Functional Pointers..");
  SetupFunctionalParameters();
  MACH3LOG_INFO("Setting up Additional Weight Pointers..");
  AddAdditionalWeightPointers();
  MACH3LOG_INFO("Setting up Kinematic Map..");
  SetupKinematicMap();
  clock.Stop();
  MACH3LOG_INFO("Finished loading MC for {}, it took {:.2f}s to finish", GetName(), clock.RealTime());
  MACH3LOG_INFO("Initialising Data");
  InititialiseData();
  MACH3LOG_INFO("=======================================================");
}

InitialiseNuOscillatorObjects()

void SampleHandlerBase::InitialiseNuOscillatorObjects ( )

protected

including Dan's magic NuOscillator

Definition at line 1040 of file SampleHandlerBase.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 M3::float_t*> OscParams = ParHandler->GetOscParsFromSampleName(SampleHandlerName);
  if (OscParams.empty()) {
    MACH3LOG_ERROR("OscParams is empty for sample '{}'.", GetName());
    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(0));
  // Add samples only if we don't use same binning
  if(!EqualBinningPerOscChannel) {
    // KS: Start from 1 because sample 0 already added
    for(int iSample = 1; iSample < GetNSamples(); iSample++) {
      Oscillator->AddSample(NuOscillatorConfigFile, GetNOscChannels(iSample));
    }
    for(int iSample = 0; iSample < GetNSamples(); iSample++) {
      for(int iChannel = 0; iChannel < GetNOscChannels(iSample); iChannel++) {
        std::vector<M3::float_t> EnergyArray;
        std::vector<M3::float_t> CosineZArray;
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wuseless-cast"
        for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
          if(MCEvents[iEvent].NominalSample != iSample) continue;
          // 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(SampleDetails[MCEvents[iEvent].NominalSample].OscChannels, MCEvents[iEvent].nupdgUnosc, MCEvents[iEvent].nupdg);
          //DB Remove NC events from the arrays which are handed to the NuOscillator objects
          if (!MCEvents[iEvent].isNC && Channel == iChannel) {
            EnergyArray.push_back(M3::float_t(MCEvents[iEvent].enu_true));
          }
        }
        std::sort(EnergyArray.begin(),EnergyArray.end());
 
        //============================================================================
        //DB Atmospheric only part, can only happen if truecz has been initialised within the experiment specific code
        if (MCEvents[0].coszenith_true != M3::_BAD_DOUBLE_) {
          for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
            if(MCEvents[iEvent].NominalSample != iSample) continue;
            // 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(SampleDetails[MCEvents[iEvent].NominalSample].OscChannels, MCEvents[iEvent].nupdgUnosc, MCEvents[iEvent].nupdg);
            //DB Remove NC events from the arrays which are handed to the NuOscillator objects
            if (!MCEvents[iEvent].isNC && Channel == iChannel) {
              CosineZArray.push_back(M3::float_t(MCEvents[iEvent].coszenith_true));
            }
          }
          std::sort(CosineZArray.begin(),CosineZArray.end());
        }
#pragma GCC diagnostic pop
        Oscillator->SetOscillatorBinning(iSample, iChannel, EnergyArray, CosineZArray);
      } // end loop over channels
    } // end loop over samples
  }
}

InitialiseSplineObject()

void SampleHandlerBase::InitialiseSplineObject ( )

protected

Definition at line 1356 of file SampleHandlerBase.cpp.

                                               {
// ************************************************
  if(auto BinnedSplines = dynamic_cast<BinnedSplineHandler*>(SplineHandler.get())) {
    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"],
                                                      (SampleHandlerName + "_SplineFile.root"), __FILE__, __LINE__);
    if(!LoadSplineFile) {
      for(int iSample = 0; iSample < GetNSamples(); iSample++) {
        std::vector<std::string> spline_filepaths = SampleDetails[iSample].spline_files;
 
        //Keep a track of the spline variables
        std::vector<std::string> SplineVarNames = {"TrueNeutrinoEnergy"};
        if (GetNDim(iSample) == 1) {
          SplineVarNames.push_back(GetKinVarName(iSample, 0));
        } else if (GetNDim(iSample) == 2) {
          SplineVarNames.push_back(GetKinVarName(iSample, 0));
          SplineVarNames.push_back(GetKinVarName(iSample, 1));
        } else {
          MACH3LOG_CRITICAL("{} Not implemented for dimension {}, will use 2D", __func__, GetNDim(iSample));
          SplineVarNames.push_back(GetKinVarName(iSample, 0));
          SplineVarNames.push_back(GetKinVarName(iSample, 1));
        }
        BinnedSplines->AddSample(SampleHandlerName, GetSampleTitle(iSample), spline_filepaths, SplineVarNames);
      }
      BinnedSplines->CountNumberOfLoadedSplines(false, 1);
      BinnedSplines->TransferToMonolith();
      if(PrepSplineFile) BinnedSplines->PrepareSplineFile(SplineFileName);
    } else {
      // KS: Skip default spline loading and use flattened spline format allowing to read stuff much faster
      BinnedSplines->LoadSplineFile(SplineFileName);
    }
    MACH3LOG_INFO("--------------------------------");
    MACH3LOG_INFO("Setup Far Detector splines");
 
    SetSplinePointers();
 
    BinnedSplines->cleanUpMemory();
  } else if (auto UnbinnedSpline = dynamic_cast<UnbinnedSplineHandler*>(SplineHandler.get())) {
    (void) UnbinnedSpline;
    SetSplinePointers();
  } else {
    MACH3LOG_ERROR("Unsupported spline type encountered.");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

InititialiseData()

virtual void SampleHandlerBase::InititialiseData ( )

protectedpure virtual

Function responsible for loading data from file or loading from file.

Implemented in PySampleHandlerBase.

IsEventSelected()

bool SampleHandlerBase::IsEventSelected ( const int  iSample, const int  iEvent  )

protected

DB Function which determines if an event is selected based on KinematicCut.

Definition at line 284 of file SampleHandlerBase.cpp.

                                                                                      {
// ************************************************
  const auto& SampleSelection = Selection[iSample];
  const int SelectionSize = static_cast<int>(SampleSelection.size());
  for (int iSelection = 0; iSelection < SelectionSize; ++iSelection) {
    const auto& Cut = SampleSelection[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 SampleHandlerBase::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 300 of file SampleHandlerBase.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 SampleHandlerBase::IsSubEventVarString

(

const std::string &

VarStr

)

const

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

Definition at line 1718 of file SampleHandlerBase.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;
}

LoadSingleSample()

void SampleHandlerBase::LoadSingleSample ( const int  iSample, const YAML::Node &  Settings  )

protected

Initialise single sample from config file.

Add new sample

Definition at line 108 of file SampleHandlerBase.cpp.

                                                                                          {
// ************************************************
  SampleInfo SingleSample;
  //SampleTitle has to be provided in the sample yaml otherwise this will throw an exception
  SingleSample.SampleTitle = Get<std::string>(SampleSettings["SampleTitle"], __FILE__ , __LINE__);
 
  Binning->SetupSampleBinning(SampleSettings["Binning"], SingleSample);
 
  auto MCFilePrefix  = Get<std::string>(SampleSettings["InputFiles"]["mtupleprefix"], __FILE__, __LINE__);
  auto MCFileSuffix  = Get<std::string>(SampleSettings["InputFiles"]["mtuplesuffix"], __FILE__, __LINE__);
  auto SplinePrefix  = Get<std::string>(SampleSettings["InputFiles"]["splineprefix"], __FILE__, __LINE__);
  auto SplineSuffix  = Get<std::string>(SampleSettings["InputFiles"]["splinesuffix"], __FILE__, __LINE__);
 
  int NChannels = static_cast<M3::int_t>(SampleSettings["OscChannels"].size());
  SingleSample.OscChannels.reserve(NChannels);
 
  int OscChannelCounter = 0;
  for (auto const &osc_channel : SampleSettings["OscChannels"]) {
    OscChannelInfo OscInfo;
    OscInfo.flavourName       = Get<std::string>(osc_channel["Name"], __FILE__ , __LINE__);
    OscInfo.flavourName_Latex = Get<std::string>(osc_channel["LatexName"], __FILE__ , __LINE__);
    OscInfo.InitPDG           = GetFromManager(osc_channel["nutype"], 0, __FILE__,__LINE__);
    OscInfo.FinalPDG          = GetFromManager(osc_channel["oscnutype"], 0, __FILE__,__LINE__);
    OscInfo.ChannelIndex      = OscChannelCounter;
 
    for (const auto& Existing : SingleSample.OscChannels) {
      if (Existing.InitPDG == OscInfo.InitPDG && Existing.FinalPDG == OscInfo.FinalPDG) {
        MACH3LOG_ERROR("Duplicate oscillation channel detected! InitPDG = {}, FinalPDG = {}"
                       "already defined in channel {} for sample {}",
                       OscInfo.InitPDG, OscInfo.FinalPDG, Existing.ChannelIndex, SingleSample.SampleTitle);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
    }
    auto MCFileNames = Get<std::vector<std::string>>(osc_channel["mtuplefile"], __FILE__ , __LINE__);
    for(size_t iFile = 0; iFile < MCFileNames.size(); iFile++){
      std::string FileName = MCFilePrefix + MCFileNames[iFile] + MCFileSuffix;
      MCFileNames[iFile] = FileName;
      FileToInitPDGMap[FileName] = NuPDG(OscInfo.InitPDG);
      FileToFinalPDGMap[FileName] = NuPDG(OscInfo.FinalPDG);
    }
 
    SingleSample.OscChannels.push_back(std::move(OscInfo));
    SingleSample.mc_files.push_back(MCFileNames);
    SingleSample.spline_files.push_back(SplinePrefix+osc_channel["splinefile"].as<std::string>()+SplineSuffix);
    OscChannelCounter++;
  }
  //Now grab the selection cuts from the manager
  for ( auto const &SelectionCuts : SampleSettings["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[iSample].emplace_back(CutObj);
  }
  Selection = StoredSelection;
  SampleDetails[iSample] = std::move(SingleSample);
}

PassesSelection()

template<typename ParT >

bool SampleHandlerBase::PassesSelection ( const ParT &  Par, std::size_t  iEvent  )

protected

Definition at line 2108 of file SampleHandlerBase.cpp.

                                                                         {
// ***************************************************************************
  bool IsSelected = true;
  if (Par.hasKinBounds) {
    const auto& kinVars = Par.KinematicVarStr;
    const auto& selection = Par.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], Par.name);
        IsSelected = false;
        break;
      }
    }
  }
  return IsSelected;
}

PrepFunctionalParameters()

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

{};

PrintIntegral()

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

Computes and prints the integral breakdown of all modes and oscillation channels for a given sample.

Definition at line 1802 of file SampleHandlerBase.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< std::unique_ptr<TH1> > > IntegralList;
  IntegralList.resize(Modes->GetNModes());
 
  std::vector<double> ChannelIntegral;
  ChannelIntegral.resize(GetNOscChannels(iSample));
  for (unsigned int i=0;i<ChannelIntegral.size();i++) {ChannelIntegral[i] = 0.;}
 
  for (int i=0;i<Modes->GetNModes();i++) {
    if (GetNDim(iSample) == 1) {
      IntegralList[i] = ReturnHistsBySelection1D(iSample, GetKinVarName(iSample, 0),1,i,WeightStyle);
    } else {
      IntegralList[i] = CastVector<TH2, TH1>(ReturnHistsBySelection2D(iSample, GetKinVarName(iSample, 0), GetKinVarName(iSample, 1),1,i,WeightStyle));
    }
  }
 
  MACH3LOG_INFO("-------------------------------------------------");
 
  if (printToFile) {
    outfile << "\\begin{table}[ht]" << std::endl;
    outfile << "\\begin{center}" << std::endl;
    outfile << "\\caption{Integral breakdown for sample: " << GetSampleTitle(iSample) << "}" << std::endl;
    outfile << "\\label{" << GetSampleTitle(iSample) << "-EventRate}" << std::endl;
 
    TString nColumns;
    for (int i=0;i<GetNOscChannels(iSample);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 :"<<GetSampleTitle(iSample)<<"\n";
  }
 
  MACH3LOG_INFO("Integral breakdown for sample: {}", GetSampleTitle(iSample));
  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(iSample); i++) {
    table_headings += fmt::format(" {:<17} |", GetFlavourName(iSample, i));
    table_footline += "--------------------";
    if (printToFile) {outfile << "&" << std::setw(space) << SampleDetails[iSample].OscChannels[i].flavourName_Latex << " ";}
    if (printToCSV)  {outcsv << GetFlavourName(iSample, 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();
  }
}

PrintRates()

void SampleHandlerBase::PrintRates ( const bool  DataOnly = false )

finalvirtual

Helper function to print rates for the samples with LLH.

Parameters

DataOnly

whether to print data only rates

Implements SampleHandlerInterface.

Definition at line 2036 of file SampleHandlerBase.cpp.

                                                      {
// ***************************************************************************
  if (!SampleHandler_data.size()) {
    MACH3LOG_ERROR("Data sample is empty!");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
  MACH3LOG_INFO("Printing for {}", GetName());
 
  if (!DataOnly) {
    const std::string sep_full(71, '-');
    MACH3LOG_INFO("{}", sep_full);
    MACH3LOG_INFO("{:<40}{:<10}{:<10}{:<10}|", "Sample", "Data", "MC", "-LLH");
  } else {
    const std::string sep_data(51, '-');
    MACH3LOG_INFO("{}", sep_data);
    MACH3LOG_INFO("{:<40}{:<10}|", "Sample", "Data");
  }
 
  double sumData = 0.0;
  double sumMC = 0.0;
  double likelihood = 0.0;
 
  for (int iSample = 0; iSample < GetNSamples(); ++iSample) {
    std::string name = GetSampleTitle(iSample);
    std::vector<double> DataArray = GetDataArray(iSample);
    double dataIntegral = std::accumulate(DataArray.begin(), DataArray.end(), 0.0);
    sumData += dataIntegral;
    if (!DataOnly) {
      std::vector<double> MCArray = GetMCArray(iSample);
      double mcIntegral = std::accumulate(MCArray.begin(), MCArray.end(), 0.0);
      sumMC += mcIntegral;
      likelihood = GetSampleLikelihood(iSample);
 
      MACH3LOG_INFO("{:<40}{:<10.2f}{:<10.2f}{:<10.2f}|", name, dataIntegral, mcIntegral, likelihood);
    } else {
      MACH3LOG_INFO("{:<40}{:<10.2f}|", name, dataIntegral);
    }
  }
  if (!DataOnly) {
    likelihood = GetLikelihood();
    MACH3LOG_INFO("{:<40}{:<10.2f}{:<10.2f}{:<10.2f}|", "Total", sumData, sumMC, likelihood);
    const std::string sep_full(71, '-');
    MACH3LOG_INFO("{}", sep_full);
  } else {
    MACH3LOG_INFO("{:<40}{:<10.2f}|", "Total", sumData);
    const std::string sep_data(51, '-');
    MACH3LOG_INFO("{}", sep_data);
  }
}

ReadConfig()

void SampleHandlerBase::ReadConfig ( )

protected

Load information about sample handler and corresponding samples from config file.

Definition at line 54 of file SampleHandlerBase.cpp.

{
  auto ModeName = Get<std::string>(SampleManager->raw()["MaCh3ModeConfig"], __FILE__ , __LINE__);
  Modes = std::make_unique<MaCh3Modes>(getenv("MACH3")+std::string("/") + ModeName);
  //SampleName has to be provided in the sample yaml otherwise this will throw an exception
  SampleHandlerName = Get<std::string>(SampleManager->raw()["SampleHandlerName"], __FILE__ , __LINE__);
 
  fTestStatistic = static_cast<TestStatistic>(SampleManager->GetMCStatLLH());
  if (CheckNodeExists(SampleManager->raw(), "LikelihoodOptions")) {
    UpdateW2 = GetFromManager<bool>(SampleManager->raw()["LikelihoodOptions"]["UpdateW2"], false);
  }
 
  if (!CheckNodeExists(SampleManager->raw(), "BinningFile")){
    MACH3LOG_ERROR("BinningFile not given in for sample handler {}, ReturnKinematicParameterBinning will not work", SampleHandlerName);
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  auto EnabledSasmples = Get<std::vector<std::string>>(SampleManager->raw()["Samples"], __FILE__ , __LINE__);
  // Get number of samples and resize relevant objects
  nSamples = static_cast<M3::int_t>(EnabledSasmples.size());
  SampleDetails.resize(nSamples);
  StoredSelection.resize(nSamples);
  for (int iSample = 0; iSample < nSamples; iSample++)
  {
    auto SampleSettings = SampleManager->raw()[EnabledSasmples[iSample]];
    LoadSingleSample(iSample, SampleSettings);
  } // end loop over enabling samples
 
  // 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 SampleHandlerBase::RegisterFunctionalParameters ( )

protectedpure virtual

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

Implemented in PySampleHandlerBase.

RegisterIndividualFunctionalParameter()

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

protected

HH - a helper function for RegisterFunctionalParameter.

Definition at line 468 of file SampleHandlerBase.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 SampleHandlerBase::ResetHistograms ( )

protected

Helper function to reset histograms.

Definition at line 457 of file SampleHandlerBase.cpp.

                                        {
// **************************************************
  // DB Reset values stored in PDF array to 0.
  // Don't openMP this; no significant gain
  const int nBins = Binning->GetNBins();
  std::fill_n(SampleHandler_array.begin(), nBins, 0.0);
  if (FirstTimeW2) {
    std::fill_n(SampleHandler_array_w2.begin(), nBins, 0.0);
  }
} // end function

ResetShifts()

virtual void SampleHandlerBase::ResetShifts ( const int  iEvent )

inlineprotectedvirtual

HH - reset the shifted values to the original values.

Definition at line 255 of file SampleHandlerBase.h.

{(void)iEvent;};

ReturnHistsBySelection1D()

std::vector< std::unique_ptr< TH1 > > SampleHandlerBase::ReturnHistsBySelection1D	(	const int	iSample,
		const std::string &	KinematicProjection,
		const int	Selection1,
		const int	Selection2 = -1,
		const int	WeightStyle = 0
	)

Definition at line 1944 of file SampleHandlerBase.cpp.

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

ReturnHistsBySelection2D()

std::vector< std::unique_ptr< TH2 > > SampleHandlerBase::ReturnHistsBySelection2D	(	const int	iSample,
		const std::string &	KinematicProjectionX,
		const std::string &	KinematicProjectionY,
		const int	Selection1,
		const int	Selection2 = -1,
		const int	WeightStyle = 0
	)

Definition at line 1983 of file SampleHandlerBase.cpp.

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

ReturnKinematicParameter() [1/2]

virtual double SampleHandlerBase::ReturnKinematicParameter ( const int  KinematicVariable, const int  iEvent  ) const

protectedpure virtual

Implemented in PySampleHandlerBase.

ReturnKinematicParameter() [2/2]

double SampleHandlerBase::ReturnKinematicParameter ( const std::string &  KinematicParameter, int  iEvent  ) const

inlineprotected

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

Definition at line 297 of file SampleHandlerBase.h.

                                                                                         {
    return ReturnKinematicParameter(ReturnKinematicParameterFromString(KinematicParameter), iEvent);
  }

ReturnKinematicParameterBinning()

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

finalprotectedvirtual

Return the binning used to draw a kinematic parameter.

Todo:

might be useful to allow overwriting this

Implements SampleHandlerInterface.

Definition at line 1665 of file SampleHandlerBase.cpp.

                                                                                                                                {
// ************************************************
  // If any of fit based variables return them
  if(Binning->IsUniform(Sample)) {
    for(int iDim = 0; iDim < GetNDim(Sample); iDim++) {
      if(KinematicParameter == GetKinVarName(Sample, iDim)) {
        return Binning->GetBinEdges(Sample, iDim);
      }
    } // end loop over Dimension
  } // end loop over IsUniform
 
  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(Sample));
  } else if (KinematicParameter == "Mode") {
    return MakeBins(Modes->GetNModes());
  }
 
  std::vector<double> BinningVect;
  // We first check if binning for a sample has been specified
  auto BinningConfig = M3OpenConfig(SampleManager->raw()["BinningFile"].as<std::string>());
  if(BinningConfig[GetSampleTitle(Sample)] && BinningConfig[GetSampleTitle(Sample)][KinematicParameter]){
    BinningVect = Get<std::vector<double>>(BinningConfig[GetSampleTitle(Sample)][KinematicParameter], __FILE__, __LINE__);
  } else {
    BinningVect = Get<std::vector<double>>(BinningConfig[KinematicParameter], __FILE__, __LINE__);
  }
 
  // Ensure binning is increasing
  auto IsIncreasing = [](const std::vector<double>& vec) {
    for (size_t i = 1; i < vec.size(); ++i) {
      if (vec[i] <= vec[i-1]) {
        return false;
      }
    }
    return true;
  };
 
  if (!IsIncreasing(BinningVect)) {
    MACH3LOG_ERROR("Binning for {} is not increasing [{}]", KinematicParameter, fmt::join(BinningVect, ", "));
    throw MaCh3Exception(__FILE__,__LINE__);
  }
 
  return BinningVect;
}

ReturnKinematicParameterFromString()

int SampleHandlerBase::ReturnKinematicParameterFromString

(

const std::string &

KinematicStr

)

const

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

Definition at line 1612 of file SampleHandlerBase.cpp.

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

ReturnKinematicVector() [1/2]

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

inlineprotectedvirtual

Definition at line 306 of file SampleHandlerBase.h.

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

ReturnKinematicVector() [2/2]

std::vector<double> SampleHandlerBase::ReturnKinematicVector ( const std::string &  KinematicParameter, const int  iEvent  ) const

inlineprotected

Definition at line 303 of file SampleHandlerBase.h.

                                                                                                       {
    return ReturnKinematicVector(ReturnKinematicVectorFromString(KinematicParameter), iEvent);
  }

ReturnKinematicVectorFromString()

int SampleHandlerBase::ReturnKinematicVectorFromString

(

const std::string &

KinematicStr

)

const

JM: Convert a kinematic vector name to its corresponding integer ID.

Definition at line 1639 of file SampleHandlerBase.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_;
}

ReturnStackedHistBySelection1D()

std::unique_ptr< THStack > SampleHandlerBase::ReturnStackedHistBySelection1D	(	const int	iSample,
		const std::string &	KinematicProjection,
		const int	Selection1,
		const int	Selection2 = -1,
		const int	WeightStyle = 0
	)

Definition at line 2014 of file SampleHandlerBase.cpp.

                                                                                                          {
// ************************************************
  auto HistList = ReturnHistsBySelection1D(iSample, KinematicProjection, Selection1, Selection2, WeightStyle);
  auto StackHist = std::make_unique<THStack>((GetSampleTitle(iSample)+"_"+KinematicProjection+"_Stack").c_str(),"");
  // Note: we use .release() to transfer ownership of each TH1 to THStack.
  for (unsigned int i=0;i<HistList.size();i++) {
    StackHist->Add(HistList[i].release());
  }
  return StackHist;
}

ReturnStackHistLegend()

const TLegend* SampleHandlerBase::ReturnStackHistLegend ( ) const

inline

Return the legend used for stacked histograms with sample info.

Definition at line 135 of file SampleHandlerBase.h.

{return THStackLeg;}

ReturnStringFromKinematicParameter()

std::string SampleHandlerBase::ReturnStringFromKinematicParameter

(

const int

KinematicVariable

)

const

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

Definition at line 1624 of file SampleHandlerBase.cpp.

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

ReturnStringFromKinematicVector()

std::string SampleHandlerBase::ReturnStringFromKinematicVector

(

const int

KinematicVariable

)

const

JM: Convert a kinematic vector integer ID to its corresponding name as a string.

Definition at line 1651 of file SampleHandlerBase.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 SampleHandlerBase::Reweight ( )

overridevirtual

main routine modifying MC prediction based on proposed parameter values

Implements SampleHandlerInterface.

Definition at line 319 of file SampleHandlerBase.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();
 
  // Update the functional parameter values to the latest proposed values
  PrepFunctionalParameters();
 
  //KS: If using CPU this does nothing, if on GPU need to make sure we finished copying memory from
  if(SplineHandler) SplineHandler->SynchroniseMemTransfer();
 
  #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 SampleHandlerBase::SaveAdditionalInfo ( TDirectory *  Dir )

finalvirtual

Store additional info in a chan.

Reimplemented from SampleHandlerInterface.

Definition at line 1310 of file SampleHandlerBase.cpp.

                                                          {
// ************************************************
  Dir->cd();
 
  YAML::Node Config = SampleManager->raw();
  TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + GetName()));
  ConfigSave.Write();
 
  for(int iSample = 0; iSample < GetNSamples(); iSample++)
  {
    std::unique_ptr<TH1> data_hist;
 
    if (GetNDim(iSample) == 1) {
      data_hist = M3::Clone<TH1D>(dynamic_cast<const TH1D*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      data_hist->GetXaxis()->SetTitle(GetKinVarName(iSample, 0).c_str());
      data_hist->GetYaxis()->SetTitle("Number of Events");
    } else if (GetNDim(iSample) == 2) {
      if(Binning->IsUniform(iSample)) {
        data_hist = M3::Clone<TH2D>(dynamic_cast<const TH2D*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      } else {
        data_hist = M3::Clone<TH2Poly>(dynamic_cast<const TH2Poly*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      }
      data_hist->GetXaxis()->SetTitle(GetKinVarName(iSample, 0).c_str());
      data_hist->GetYaxis()->SetTitle(GetKinVarName(iSample, 1).c_str());
      data_hist->GetZaxis()->SetTitle("Number of Events");
    } else {
      data_hist = M3::Clone<TH1D>(dynamic_cast<const TH1D*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      int nbins = Binning->GetNBins(iSample);
      for(int iBin = 0; iBin < nbins; iBin++) {
        auto BinName = Binning->GetBinName(iSample, iBin);
        data_hist->GetXaxis()->SetBinLabel(iBin+1, BinName.c_str());
      }
      data_hist->GetYaxis()->SetTitle("Number of Events");
    }
 
    if (!data_hist) {
      MACH3LOG_ERROR("nullptr data hist :(");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
 
    data_hist->SetTitle(("data_" + GetSampleTitle(iSample)).c_str());
    data_hist->Write();
  }
}

SetBinning()

void SampleHandlerBase::SetBinning ( )

protected

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

Definition at line 750 of file SampleHandlerBase.cpp.

                                   {
// ************************************************
  for(int iSample = 0; iSample < GetNSamples(); iSample++)
  {
    int Dimension = GetNDim(iSample);
    std::string HistTitle = GetSampleTitle(iSample);
 
    auto* SamDet = &SampleDetails[iSample];
    if(Dimension == 1) {
      auto XVec = Binning->GetBinEdges(iSample, 0);
      SamDet->DataHist = new TH1D(("d" + HistTitle).c_str(), HistTitle.c_str(), static_cast<int>(XVec.size()-1), XVec.data());
      SamDet->MCHist   = new TH1D(("h" + HistTitle).c_str(), HistTitle.c_str(), static_cast<int>(XVec.size()-1), XVec.data());
      SamDet->W2Hist   = new TH1D(("w" + HistTitle).c_str(), HistTitle.c_str(), static_cast<int>(XVec.size()-1), XVec.data());
 
      // Set all titles so most of projections don't have empty titles...
      SamDet->DataHist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->DataHist->GetYaxis()->SetTitle("Events");
      SamDet->MCHist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->MCHist->GetYaxis()->SetTitle("Events");
      SamDet->W2Hist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->W2Hist->GetYaxis()->SetTitle("Events");
    } else if (Dimension == 2){
      if(Binning->IsUniform(iSample)) {
        auto XVec = Binning->GetBinEdges(iSample, 0);
        auto YVec = Binning->GetBinEdges(iSample, 1);
        int nX = static_cast<int>(XVec.size() - 1);
        int nY = static_cast<int>(YVec.size() - 1);
 
        SamDet->DataHist = new TH2D(("d" + HistTitle).c_str(), HistTitle.c_str(), nX, XVec.data(), nY, YVec.data());
        SamDet->MCHist   = new TH2D(("h" + HistTitle).c_str(), HistTitle.c_str(), nX, XVec.data(), nY, YVec.data());
        SamDet->W2Hist   = new TH2D(("w" + HistTitle).c_str(), HistTitle.c_str(), nX, XVec.data(), nY, YVec.data());
      } else {
        auto AddBinsToTH2Poly = [](TH2Poly* hist, const std::vector<BinInfo>& bins) {
          for (const auto& bin : bins) {
            double xLow  = bin.Extent[0][0];
            double xHigh = bin.Extent[0][1];
            double yLow  = bin.Extent[1][0];
            double yHigh = bin.Extent[1][1];
 
            double x[4] = {xLow, xHigh, xHigh, xLow};
            double y[4] = {yLow, yLow, yHigh, yHigh};
 
            hist->AddBin(4, x, y);
          }
        };
        // Create all three histograms
        SamDet->DataHist = new TH2Poly();
        SamDet->DataHist->SetName(("d" + HistTitle).c_str());
        SamDet->DataHist->SetTitle(HistTitle.c_str());
 
        SamDet->MCHist = new TH2Poly();
        SamDet->MCHist->SetName(("h" + HistTitle).c_str());
        SamDet->MCHist->SetTitle(HistTitle.c_str());
 
        SamDet->W2Hist = new TH2Poly();
        SamDet->W2Hist->SetName(("w" + HistTitle).c_str());
        SamDet->W2Hist->SetTitle(HistTitle.c_str());
 
        // Add bins to each
        AddBinsToTH2Poly(static_cast<TH2Poly*>(SamDet->DataHist), Binning->GetNonUniformBins(iSample));
        AddBinsToTH2Poly(static_cast<TH2Poly*>(SamDet->MCHist),   Binning->GetNonUniformBins(iSample));
        AddBinsToTH2Poly(static_cast<TH2Poly*>(SamDet->W2Hist),   Binning->GetNonUniformBins(iSample));
      }
 
      // Set all titles so most of projections don't have empty titles...
      SamDet->DataHist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->DataHist->GetYaxis()->SetTitle(SamDet->VarStr[1].c_str());
      SamDet->MCHist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->MCHist->GetYaxis()->SetTitle(SamDet->VarStr[1].c_str());
      SamDet->W2Hist->GetXaxis()->SetTitle(SamDet->VarStr[0].c_str());
      SamDet->W2Hist->GetYaxis()->SetTitle(SamDet->VarStr[1].c_str());
    } else {
      int nbins = Binning->GetNBins(iSample);
      SamDet->DataHist = new TH1D(("d" + HistTitle).c_str(), HistTitle.c_str(), nbins, 0, nbins);
      SamDet->MCHist   = new TH1D(("h" + HistTitle).c_str(), HistTitle.c_str(), nbins, 0, nbins);
      SamDet->W2Hist   = new TH1D(("w" + HistTitle).c_str(), HistTitle.c_str(), nbins, 0, nbins);
 
      for(int iBin = 0; iBin < nbins; iBin++) {
        auto BinName = Binning->GetBinName(iSample, iBin);
        SamDet->DataHist->GetXaxis()->SetBinLabel(iBin+1, BinName.c_str());
        SamDet->MCHist->GetXaxis()->SetBinLabel(iBin+1, BinName.c_str());
        SamDet->W2Hist->GetXaxis()->SetBinLabel(iBin+1, BinName.c_str());
      }
 
      // Set all titles so most of projections don't have empty titles...
      SamDet->DataHist->GetYaxis()->SetTitle("Events");
      SamDet->MCHist->GetYaxis()->SetTitle("Events");
      SamDet->W2Hist->GetYaxis()->SetTitle("Events");
    }
 
    SamDet->DataHist->SetDirectory(nullptr);
    SamDet->MCHist->SetDirectory(nullptr);
    SamDet->W2Hist->SetDirectory(nullptr);
  }
 
  //Set the number of X and Y bins now
  SetupReweightArrays();
  FindNominalBinAndEdges();
}

SetSplinePointers()

void SampleHandlerBase::SetSplinePointers ( )

protected

Set pointers for each event to appropriate weights, for unbinned based on event number while for binned based on other kinematical properties.

Definition at line 1228 of file SampleHandlerBase.cpp.

                                          {
// ************************************************
  //Now loop over events and get the spline bin for each event
  if (auto BinnedSpline = dynamic_cast<BinnedSplineHandler*>(SplineHandler.get())) {
    bool ThrowCrititcal = true;
    auto SplineParsVec = ParHandler->GetSplineParsFromSampleName(SampleHandlerName);
    for (unsigned int j = 0; j < GetNEvents(); ++j) {
      auto EventSplines = GetSplineBins(j, BinnedSpline, ThrowCrititcal);
      const int NSplines = static_cast<int>(EventSplines.size());
      if(NSplines == 0) continue;
      auto& w_pointers = MCEvents[j].total_weight_pointers;
      w_pointers.reserve(w_pointers.size() + NSplines);
 
      for(int spline = 0; spline < NSplines; spline++) {
        int SystIndex = EventSplines[spline][2];
 
        bool IsSelected = PassesSelection(SplineParsVec[SystIndex], j);
        // Need to then break the event loop
        if(!IsSelected){
          MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", j, SplineParsVec[SystIndex].name);
          continue;
        }
        //Event Splines indexed as: sample name, oscillation channel, syst, mode, etrue, var1, var2 (var2 is a dummy 0 for 1D splines)
        w_pointers.push_back(BinnedSpline->RetPointer(EventSplines[spline][0], EventSplines[spline][1],
                                                      EventSplines[spline][2], EventSplines[spline][3],
                                                      EventSplines[spline][4], EventSplines[spline][5],
                                                      EventSplines[spline][6]));
      } // end loop over splines
      w_pointers.shrink_to_fit();
    } // end loop over events
  } else if (auto UnbinnedSpline = dynamic_cast<UnbinnedSplineHandler*>(SplineHandler.get())) {
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) {
      MCEvents[iEvent].total_weight_pointers.push_back(UnbinnedSpline->RetPointer(iEvent));
    }
  } else {
    MACH3LOG_ERROR("Not supported splines");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

SetupExperimentMC()

virtual int SampleHandlerBase::SetupExperimentMC ( )

protectedpure virtual

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

Implemented in PySampleHandlerBase.

SetupFunctionalParameters()

void SampleHandlerBase::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.

Todo:

KS: Instead of clearing it they should not be class members, so we must fix it in future

Definition at line 488 of file SampleHandlerBase.cpp.

                                                  {
// **************************************************
  funcParsGrid.resize(GetNEvents());
  if(ParHandler == nullptr) return;
  funcParsVec = ParHandler->GetFunctionalParametersFromSampleName(SampleHandlerName);
  // RegisterFunctionalParameters is implemented in experiment-specific code,
  // which calls RegisterIndividualFuncPar to populate funcParsNamesMap, funcParsNamesVec, and funcParsFuncMap
  RegisterFunctionalParameters();
  funcParsMap.resize(funcParsNamesMap.size());
 
  // For every functional parameter in XsecCov that matches the name in funcParsNames, add it to the map
  for (FunctionalParameter& fp : funcParsVec) {
    auto it = funcParsNamesMap.find(fp.name);
    // If we don't find a match, we need to throw an error
    if (it == funcParsNamesMap.end()) {
      MACH3LOG_ERROR("Functional parameter {} not found, did you define it in RegisterFunctionalParameters()?", fp.name);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    const std::size_t key = static_cast<std::size_t>(it->second);
    MACH3LOG_INFO("Adding functional parameter: {} to funcParsMap with key: {}",fp.name, key);
 
    const int ikey = it->second;
    fp.funcPtr = &funcParsFuncMap[ikey];
 
    funcParsMap[key].valuePtr = fp.valuePtr;
    funcParsMap[key].funcPtr  = fp.funcPtr;
  }
 
  // Mostly the same as CalcNormsBins
  // 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
      const int Mode = static_cast<int>(std::round(ReturnKinematicParameter("Mode", static_cast<int>(iEvent))));
      bool ModeMatch = MatchCondition(it->modes, Mode);
      if (!ModeMatch) {
        MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, Mode, it->name);
        continue;
      }
      // Now check whether within kinematic bounds
      bool IsSelected = PassesSelection((*it), iEvent);
      // Need to then break the event loop
      if(!IsSelected){
        MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", iEvent, it->name);
        continue;
      }
      const std::size_t key = static_cast<std::size_t>(funcParsNamesMap[it->name]);
      funcParsGrid[iEvent].push_back(&funcParsMap[key]);
    }
  }
  MACH3LOG_INFO("Finished setting up functional parameters");
 
  CleanVector(funcParsNamesVec);
 
  funcParsNamesMap.clear();
  funcParsNamesMap.rehash(0);
}

SetupKinematicMap()

void SampleHandlerBase::SetupKinematicMap ( )

protected

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

Definition at line 208 of file SampleHandlerBase.cpp.

                                          {
// ************************************************
  if(KinematicParameters == nullptr || ReversedKinematicParameters == nullptr) {
    MACH3LOG_ERROR("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__);
    }
  }
  // KS: Ensure some MaCh3 specific variables are defined
  std::vector<std::string> Vars = {"Mode", "OscillationChannel", "TargetNucleus"};
  for(size_t iVar = 0; iVar < Vars.size(); iVar++) {
    try {
      ReturnKinematicParameterFromString(Vars[iVar]);
    } catch (const MaCh3Exception&) {
      MACH3LOG_ERROR("MaCh3 expected variable: {} not found in KinematicParameters.", Vars[iVar]);
      MACH3LOG_ERROR("All keys in KinematicParameters:");
      for (const auto& pair : *KinematicParameters) {
        MACH3LOG_ERROR("Key: {}", pair.first);
      }
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
}

SetupMC()

virtual void SampleHandlerBase::SetupMC ( )

protectedpure virtual

Function which translates experiment struct into core struct.

Implemented in PySampleHandlerBase.

SetupNormParameters()

void SampleHandlerBase::SetupNormParameters ( )

protected

Setup the norm parameters by assigning each event with bin.

Definition at line 628 of file SampleHandlerBase.cpp.

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

SetupNuOscillatorPointers()

void SampleHandlerBase::SetupNuOscillatorPointers ( )

protected

Initialise pointer to oscillation weight to NuOscillator object.

Definition at line 1105 of file SampleHandlerBase.cpp.

                                                  {
// ************************************************
  auto AddOscPointer = GetFromManager<bool>(SampleManager->raw()["NuOsc"]["AddOscPointer"], true, __FILE__ , __LINE__);
  // Sometimes we don't want to add osc pointer to allow for some experiment specific handling of osc weight, like for example being able to shift osc weigh
  if(!AddOscPointer) {
    return;
  }
  for (unsigned int iEvent=0;iEvent<GetNEvents();iEvent++) {
   const M3::float_t* osc_w_pointer = GetNuOscillatorPointers(iEvent);
 
   // KS: Do not add unity
   if (osc_w_pointer != &M3::Unity) {
     MCEvents[iEvent].total_weight_pointers.push_back(osc_w_pointer);
   }
  }
}

SetupOscParameters()

void SampleHandlerBase::SetupOscParameters ( )

protected

Setup the osc parameters.

Definition at line 589 of file SampleHandlerBase.cpp.

                                           {
// ***************************************************************************
  // KS: Only make sense to setup osc if you have ParHandler
  if(ParHandler == nullptr ) return;
 
  auto OscParams = ParHandler->GetOscParsFromSampleName(SampleHandlerName);
  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("SampleHandler {} has {} osc params, while shared NuOsc has {} osc params", GetName(),
                        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");
    if(CheckNodeExists(SampleManager->raw(), "NuOsc")){
      MACH3LOG_ERROR("However config for SampleHandler {} has 'NuOsc' field", GetName());
      MACH3LOG_ERROR("This may indicate misconfiguration");
      MACH3LOG_ERROR("Either remove 'NuOsc' field from SampleHandler config or check your model.yaml and include oscillation for sample");
      throw MaCh3Exception(__FILE__, __LINE__);
    }
  }
}

SetupReweightArrays()

void SampleHandlerBase::SetupReweightArrays ( )

protected

Initialise data, MC and W2 histograms.

Definition at line 742 of file SampleHandlerBase.cpp.

                                            {
// ************************************************
  SampleHandler_array = std::vector<double>(Binning->GetNBins(),0);
  SampleHandler_array_w2 = std::vector<double>(Binning->GetNBins(),0);
  SampleHandler_data = std::vector<double>(Binning->GetNBins(),0);
}

SetupSplines()

virtual void SampleHandlerBase::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 PySampleHandlerBase.

Member Data Documentation

_modeNomWeightMap

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

protected

Definition at line 400 of file SampleHandlerBase.h.

Binning

std::unique_ptr<BinningHandler> SampleHandlerBase::Binning

protected

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

Definition at line 347 of file SampleHandlerBase.h.

FileToFinalPDGMap

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

private

Definition at line 424 of file SampleHandlerBase.h.

FileToInitPDGMap

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

private

Definition at line 423 of file SampleHandlerBase.h.

FirstTimeW2

bool SampleHandlerBase::FirstTimeW2

protected

KS:Super hacky to update W2 or not.

Definition at line 408 of file SampleHandlerBase.h.

funcParsFuncMap

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

protected

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

Definition at line 279 of file SampleHandlerBase.h.

funcParsGrid

std::vector<std::vector<FunctionalShifter*> > SampleHandlerBase::funcParsGrid

protected

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

Definition at line 262 of file SampleHandlerBase.h.

funcParsMap

std::vector<FunctionalShifter> SampleHandlerBase::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 267 of file SampleHandlerBase.h.

funcParsNamesMap

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

protected

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

Definition at line 276 of file SampleHandlerBase.h.

funcParsNamesVec

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

protected

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

Definition at line 281 of file SampleHandlerBase.h.

funcParsVec

std::vector<FunctionalParameter> SampleHandlerBase::funcParsVec

protected

HH - a vector that stores all the FuncPars struct.

Todo:

KS: Below functional variables are used only on setup, thus we should refactor them in such a way that they are removed as class members but this would be breaking change thus keep it for the time being.

Definition at line 273 of file SampleHandlerBase.h.

KinematicParameters

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

protected

Mapping between string and kinematic enum.

Definition at line 387 of file SampleHandlerBase.h.

KinematicVectors

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

protected

Definition at line 392 of file SampleHandlerBase.h.

MCEvents

std::vector<EventInfo> SampleHandlerBase::MCEvents

protected

Stores information about every MC event.

Definition at line 358 of file SampleHandlerBase.h.

Oscillator

std::shared_ptr<OscillationHandler> SampleHandlerBase::Oscillator

protected

Contains oscillator handling calculating oscillation probabilities.

Definition at line 219 of file SampleHandlerBase.h.

ParHandler

ParameterHandlerGeneric* SampleHandlerBase::ParHandler = nullptr

protected

ETA - All experiments will need an xsec, det and osc cov.

Definition at line 366 of file SampleHandlerBase.h.

ReversedKinematicParameters

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

protected

Mapping between kinematic enum and string.

Definition at line 389 of file SampleHandlerBase.h.

ReversedKinematicVectors

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

protected

Definition at line 393 of file SampleHandlerBase.h.

SampleDetails

std::vector<SampleInfo> SampleHandlerBase::SampleDetails

protected

Stores info about currently initialised sample.

Definition at line 360 of file SampleHandlerBase.h.

SampleHandler_array

std::vector<double> SampleHandlerBase::SampleHandler_array

protected

DB Array to be filled after reweighting.

Definition at line 349 of file SampleHandlerBase.h.

SampleHandler_array_w2

std::vector<double> SampleHandlerBase::SampleHandler_array_w2

protected

KS Array used for MC stat.

Definition at line 351 of file SampleHandlerBase.h.

SampleHandler_data

std::vector<double> SampleHandlerBase::SampleHandler_data

protected

DB Array to be filled in AddData.

Definition at line 353 of file SampleHandlerBase.h.

SampleHandlerName

std::string SampleHandlerBase::SampleHandlerName

protected

A unique ID for each sample based on which we can define what systematic should be applied.

Definition at line 370 of file SampleHandlerBase.h.

SampleManager

std::unique_ptr<Manager> SampleHandlerBase::SampleManager

protected

The manager object used to read the sample yaml file.

Definition at line 397 of file SampleHandlerBase.h.

Selection

std::vector< std::vector< KinematicCut > > SampleHandlerBase::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 383 of file SampleHandlerBase.h.

SplineHandler

std::unique_ptr<SplineBase> SampleHandlerBase::SplineHandler

protected

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

Definition at line 216 of file SampleHandlerBase.h.

StoredSelection

std::vector< std::vector< KinematicCut > > SampleHandlerBase::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 380 of file SampleHandlerBase.h.

THStackLeg

TLegend* SampleHandlerBase::THStackLeg = nullptr

protected

DB Miscellaneous Variables.

Definition at line 404 of file SampleHandlerBase.h.

UpdateW2

bool SampleHandlerBase::UpdateW2

protected

KS:Super hacky to update W2 or not.

Definition at line 410 of file SampleHandlerBase.h.

---

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

• Samples/SampleHandlerBase.h
• Samples/SampleHandlerBase.cpp