SampleHandlerFD Class Reference

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

#include <Samples/SampleHandlerFD.h>

Inheritance diagram for SampleHandlerFD:

Collaboration diagram for SampleHandlerFD:

Public Member Functions
	SampleHandlerFD (std::string ConfigFileName, ParameterHandlerGeneric *xsec_cov, const std::shared_ptr< OscillationHandler > &OscillatorObj_=nullptr)
	Constructor. More...
virtual	~SampleHandlerFD ()
	destructor More...
int	GetNDim (const int Sample) const override
	DB Function to differentiate 1D or 2D binning. More...
std::string	GetName () const override
std::string	GetSampleTitle (const int Sample) const override
std::string	GetKinVarName (const int iSample, const int Dimension) const override
	Return Kinematic Variable name for specified sample and dimension for example "Reconstructed_Neutrino_Energy". More...
std::string	GetXBinVarName (const int Sample) const
std::string	GetYBinVarName (const int Sample) const
const BinningHandler *	GetBinningHandler () const
	Get pointer to binning handler. More...
void	PrintIntegral (const int iSample, const TString &OutputName="/dev/null", const int WeightStyle=0, const TString &OutputCSVName="/dev/null")
void	AddData (const int Sample, TH1 *Data)
void	AddData (const int Sample, const std::vector< double > &Data_Array)
void	PrintRates (const bool DataOnly=false) override
	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...
TH1 *	GetDataHist (const int Sample) override
	Get Data histogram. More...
TH1 *	GetDataHist (const std::string &Sample)
TH1 *	GetMCHist (const int Sample) override
	Get MC histogram. More...
TH1 *	GetMCHist (const std::string &Sample)
TH1 *	GetW2Hist (const int Sample) override
	Get W2 histogram. More...
TH1 *	GetW2Hist (const std::string &Sample)
void	Reweight () override
M3::float_t	GetEventWeight (const int iEntry)
void	InitialiseNuOscillatorObjects ()
	including Dan's magic NuOscillator More...
void	SetupNuOscillatorPointers ()
const M3::float_t *	GetNuOscillatorPointers (const int iEvent) const
void	ReadConfig ()
void	LoadSingleSample (const int iSample, const YAML::Node &Settings)
int	GetNOscChannels (const int iSample) const override
std::string	GetFlavourName (const int iSample, const int iChannel) const override
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...
TH1 *	Get1DVarHist (const int iSample, const std::string &ProjectionVar, const std::vector< KinematicCut > &EventSelectionVec={}, int WeightStyle=0, TAxis *Axis=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec={}) override
TH2 *	Get2DVarHist (const int iSample, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &EventSelectionVec={}, int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec={}) override
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, TH2D *_h2DVar, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &SubEventSelectionVec={}, int WeightStyle=0)
TH1 *	Get1DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr) override
TH2 *	Get2DVarHistByModeAndChannel (const int iSample, const std::string &ProjectionVar_StrX, const std::string &ProjectionVar_StrY, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr) override
TH1 *	GetModeHist1D (const int iSample, int s, int m, int style=0)
TH2 *	GetModeHist2D (const int iSample, int s, int m, int style=0)
std::vector< TH1 * >	ReturnHistsBySelection1D (const int iSample, const std::string &KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=nullptr)
std::vector< TH2 * >	ReturnHistsBySelection2D (const int iSample, const std::string &KinematicProjectionX, const std::string &KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *XAxis=nullptr, TAxis *YAxis=nullptr)
THStack *	ReturnStackedHistBySelection1D (const int iSample, const std::string &KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=nullptr)
TLegend *	ReturnStackHistLegend ()
int	ReturnKinematicParameterFromString (const std::string &KinematicStr) const
	ETA function to generically convert a string from xsec cov to a kinematic type. More...
std::string	ReturnStringFromKinematicParameter (const int KinematicVariable) const
	ETA function to generically convert a kinematic type from xsec cov to a string. More...
void	SaveAdditionalInfo (TDirectory *Dir) override
	Store additional info in a chan. More...
int	ReturnKinematicVectorFromString (const std::string &KinematicStr) const
std::string	ReturnStringFromKinematicVector (const int KinematicVariable) const
bool	IsSubEventVarString (const std::string &VarStr)
	JM Check if a kinematic parameter string corresponds to a subevent-level variable. More...
std::vector< double >	GetDataArray () const
	Return array storing data entries for every bin. More...
std::vector< double >	GetMCArray () const
	Return array storing MC entries for every bin. More...
std::vector< double >	GetW2Array () const
	Return array storing W2 entries for every bin. More...
std::vector< double >	GetArrayForSample (const int Sample, const double *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 SampleHandlerBase
	SampleHandlerBase ()
	The main constructor. More...
virtual	~SampleHandlerBase ()
	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
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	SetupFDMC ()=0
	Function which translates experiment struct into core struct. More...
void	Initialise ()
	Function which does a lot of the lifting regarding the workflow in creating different MC objects. More...
void	SetSplinePointers ()
	Set pointers for each event to appropriate weights, for unbinned based on event number while for binned based on other kinematical properties. 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, where Selection double looks like {{ND280KinematicTypes Var1, douuble LowBound}. More...
bool	IsSubEventSelected (const std::vector< KinematicCut > &SubEventCuts, const int iEvent, unsigned const int iSubEvent, size_t nsubevents)
	JM Function which determines if a subevent is selected. More...
virtual void	ResetShifts (const int iEvent)
	HH - reset the shifted values to the original values. 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
	Calculate the total weight weight for a given event. More...
virtual void	CalcWeightFunc (int iEvent)
	Calculate weights for function parameters. More...
virtual double	ReturnKinematicParameter (std::string KinematicParamter, int iEvent)=0
	Return the value of an associated kinematic parameter for an event. More...
virtual double	ReturnKinematicParameter (int KinematicVariable, int iEvent)=0
virtual std::vector< double >	ReturnKinematicVector (std::string KinematicParameter, int iEvent)
virtual std::vector< double >	ReturnKinematicVector (int KinematicVariable, int iEvent)
std::vector< double >	ReturnKinematicParameterBinning (const int Sample, const std::string &KinematicParameter) const override
	Return the binning used to draw a kinematic parameter. More...
virtual const double *	GetPointerToKinematicParameter (std::string KinematicParamter, int iEvent)=0
virtual const double *	GetPointerToKinematicParameter (double KinematicVariable, int iEvent)=0
const double *	GetPointerToOscChannel (const int iEvent) const
	Get pointer to oscillation channel associated with given event. Osc channel is const. More...
void	SetupNormParameters ()
	Setup the norm parameters by assigning each event with bin. More...
void	FillHist (const int Sample, TH1 *Hist, double *Array)
	Fill a histogram with the event-level information used in the fit. More...
void	FillArray ()
	DB Nice new multi-threaded function which calculates the event weights and fills the relevant bins of an array. More...
void	ResetHistograms ()
	Helper function to reset histograms. More...
void	InitialiseSplineObject ()
NuPDG	GetInitPDGFromFileName (const std::string &FileName) const
	Retrieve the initial neutrino PDG code associated with a given input file name. More...
NuPDG	GetFinalPDGFromFileName (const std::string &FileName) const
	Retrieve the final neutrino PDG code associated with a given input file name. More...
Protected Member Functions inherited from SampleHandlerBase
void	QuietPlease ()
	CW: Redirect std::cout to silence some experiment specific libraries. More...
void	NowTalk ()
	CW: Redirect std::cout to silence some experiment specific libraries. More...
template<typename T >
bool	MatchCondition (const std::vector< T > &allowedValues, const T &value)
	check if event is affected by following conditions, for example pdg, or modes etc More...

Protected Attributes
std::unique_ptr< 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...
double *	SampleHandlerFD_array
	DB Array to be filled after reweighting. More...
double *	SampleHandlerFD_array_w2
	KS Array used for MC stat. More...
double *	SampleHandlerFD_data
	DB Array to be filled in AddData. More...
std::vector< EventInfo >	MCSamples
	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 SampleHandlerBase
TestStatistic	fTestStatistic
	Test statistic tells what kind of likelihood sample is using. More...
std::streambuf *	buf
	Keep the cout buffer. More...
std::streambuf *	errbuf
	Keep the cerr buffer. More...
M3::int_t	nSamples
	Contains how many samples we've got. More...
unsigned int	nEvents
	Number of MC events are there. More...
std::unique_ptr< MaCh3Modes >	Modes
	Holds information about used Generator and MaCh3 modes. More...

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

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

Detailed Description

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

Author

Dan Barrow

Ed Atkin

Definition at line 21 of file SampleHandlerFD.h.

Member Enumeration Documentation

FDPlotType

enum SampleHandlerFD::FDPlotType

private

Enumerator
kModePlot
kOscChannelPlot

Definition at line 403 of file SampleHandlerFD.h.

                  {
    kModePlot = 0,
    kOscChannelPlot = 1
  };

Constructor & Destructor Documentation

SampleHandlerFD()

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

Constructor.

Parameters

ConfigFileName

Name of config to initialise the sample object

Definition at line 12 of file SampleHandlerFD.cpp.

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

~SampleHandlerFD()

SampleHandlerFD::~SampleHandlerFD ( )

virtual

destructor

Definition at line 48 of file SampleHandlerFD.cpp.

                                  {
  MACH3LOG_DEBUG("I'm deleting SampleHandlerFD");
  
  if (SampleHandlerFD_array != nullptr) delete[] SampleHandlerFD_array;
  if (SampleHandlerFD_array_w2 != nullptr) delete[] SampleHandlerFD_array_w2;
  //ETA - there is a chance that you haven't added any data...
  if (SampleHandlerFD_data != nullptr) delete[] SampleHandlerFD_data;
 
  if(THStackLeg != nullptr) delete THStackLeg;
}

Member Function Documentation

AddAdditionalWeightPointers()

virtual void SampleHandlerFD::AddAdditionalWeightPointers ( )

protectedpure virtual

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

Implemented in PySampleHandlerFD.

AddData() [1/2]

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

Definition at line 998 of file SampleHandlerFD.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__);
  }
 
  for (int idx = Start; idx < End; ++idx) {
    SampleHandlerFD_data[idx] = Data_Array[idx - Start];
  }
 
  FillHist(Sample, SampleDetails[Sample].DataHist, SampleHandlerFD_data);
}

AddData() [2/2]

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

Definition at line 935 of file SampleHandlerFD.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(SampleHandlerFD_data == nullptr) {
    MACH3LOG_ERROR("SampleHandlerFD_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__);
    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
      SampleHandlerFD_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(xBin + 1);
    }
    SampleDetails[Sample].DataHist->GetXaxis()->SetTitle(GetXBinVarName(Sample).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__);
      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
          SampleHandlerFD_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(xBin + 1, yBin + 1);
        }
      }
    } else{
      ChecHistType("TH2Poly", Dim, SampleDetails[Sample].DataHist, __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
        SampleHandlerFD_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(iBin + 1);
      }
    }
 
    SampleDetails[Sample].DataHist->GetXaxis()->SetTitle(GetXBinVarName(Sample).c_str());
    SampleDetails[Sample].DataHist->GetYaxis()->SetTitle(GetYBinVarName(Sample).c_str());
    SampleDetails[Sample].DataHist->GetZaxis()->SetTitle("Number of Events");
  } else {
    ChecHistType("TH1D", Dim, SampleDetails[Sample].DataHist, __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
      SampleHandlerFD_data[idx] = SampleDetails[Sample].DataHist->GetBinContent(iBin + 1);
    }
  }
}

ApplyShifts()

void SampleHandlerFD::ApplyShifts ( const int  iEvent )

protectedvirtual

ETA - generic function applying shifts.

Definition at line 562 of file SampleHandlerFD.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);
  }
}

ApplyTemporarySelection()

std::vector< std::vector< KinematicCut > > SampleHandlerFD::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.

Definition at line 1362 of file SampleHandlerFD.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 > SampleHandlerFD::BuildModeChannelSelection	(	const int	iSample,
		const int	kModeToFill,
		const int	kChannelToFill
	)		const

Definition at line 1691 of file SampleHandlerFD.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 SampleHandlerFD::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 642 of file SampleHandlerFD.cpp.

                                                                                                                        {
// ************************************************
  #ifdef DEBUG
  std::vector<int> VerboseCounter(norm_parameters.size(), 0);
  #endif
  for(unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent){
    std::vector< int > NormBins = {};
    if (ParHandler) {
      // Skip oscillated NC events
      // Not strictly needed, but these events don't get included in oscillated predictions, so
      // no need to waste our time calculating and storing information about xsec parameters
      // that will never be used.
      if (MCSamples[iEvent].isNC && (*MCSamples[iEvent].nupdg != *MCSamples[iEvent].nupdgUnosc) ) {
        MACH3LOG_TRACE("Event {}, missed NC/signal check", iEvent);
        continue;
      } //DB Abstract check on MaCh3Modes to determine which apply to neutral current
      for (std::vector<NormParameter>::iterator it = norm_parameters.begin(); it != norm_parameters.end(); ++it) {
        //Now check that the target of an interaction matches with the normalisation parameters
        bool TargetMatch = MatchCondition((*it).targets, *(MCSamples[iEvent].Target));
        if (!TargetMatch) {
          MACH3LOG_TRACE("Event {}, missed target check ({}) for dial {}", iEvent, *(MCSamples[iEvent].Target), (*it).name);
          continue;
        }
 
        //Now check that the neutrino flavour in an interaction matches with the normalisation parameters
        bool FlavourMatch = MatchCondition((*it).pdgs, *(MCSamples[iEvent].nupdg));
        if (!FlavourMatch) {
          MACH3LOG_TRACE("Event {}, missed PDG check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdg), (*it).name);
          continue;
        }
 
        //Now check that the unoscillated neutrino flavour in an interaction matches with the normalisation parameters
        bool FlavourUnoscMatch = MatchCondition((*it).preoscpdgs, *(MCSamples[iEvent].nupdgUnosc));
        if (!FlavourUnoscMatch){
          MACH3LOG_TRACE("Event {}, missed FlavourUnosc check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdgUnosc), (*it).name);
          continue;
        }
 
        //Now check that the mode of an interaction matches with the normalisation parameters
        bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
        if (!ModeMatch) {
          MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
          continue;
        }
 
        //Now check whether the norm has kinematic bounds
        //i.e. does it only apply to events in a particular kinematic region?
        // Now check whether within kinematic bounds
        bool IsSelected = 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 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 DEBUG
  MACH3LOG_DEBUG("┌──────────────────────────────────────────────────────────┐");
  for (std::size_t i = 0; i < norm_parameters.size(); ++i) {
    const auto& norm = norm_parameters[i];
    double eventRatio = static_cast<double>(VerboseCounter[i]) / static_cast<double>(GetNEvents());
 
    MACH3LOG_DEBUG("│ Param {:<15}, affects {:<8} events ({:>6.2f}%) │",
                  ParHandler->GetParFancyName(norm.index), VerboseCounter[i], eventRatio);
  }
  MACH3LOG_DEBUG("└──────────────────────────────────────────────────────────┘");
  #endif
}

CalcWeightFunc()

virtual void SampleHandlerFD::CalcWeightFunc ( int  iEvent )

inlineprotectedvirtual

Calculate weights for function parameters.

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

Definition at line 284 of file SampleHandlerFD.h.

{return; (void)iEvent;};

CalcWeightTotal()

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

protected

Calculate the total weight weight for a given event.

Definition at line 583 of file SampleHandlerFD.cpp.

                                                                                    {
// ***************************************************************************
  M3::float_t TotalWeight = 1.0;
  const int nNorms = static_cast<int>(MCEvent->norm_pointers.size());
  //Loop over stored normalisation and function pointers
  #ifdef MULTITHREAD
  #pragma omp simd reduction(*:TotalWeight)
  #endif
  for (int iParam = 0; iParam < nNorms; ++iParam)
  {
    #pragma GCC diagnostic push
    #pragma GCC diagnostic ignored "-Wuseless-cast"
    TotalWeight *= static_cast<M3::float_t>(*(MCEvent->norm_pointers[iParam]));
    #pragma GCC diagnostic pop
  }
 
  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 SampleHandlerFD::Fill1DSubEventHist	(	const int	iSample,
		TH1D *	_h1DVar,
		const std::string &	ProjectionVar,
		const std::vector< KinematicCut > &	SubEventSelectionVec = {},
		int	WeightStyle = 0
	)

Definition at line 1428 of file SampleHandlerFD.cpp.

                                                                                                                 {
// ************************************************
  int ProjectionVar_Int = ReturnKinematicVectorFromString(ProjectionVar_Str);
 
  //JM Loop over all events
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    const int EventSample = MCSamples[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 SampleHandlerFD::Fill2DSubEventHist	(	const int	iSample,
		TH2D *	_h2DVar,
		const std::string &	ProjectionVarX,
		const std::string &	ProjectionVarY,
		const std::vector< KinematicCut > &	SubEventSelectionVec = {},
		int	WeightStyle = 0
	)

Definition at line 1508 of file SampleHandlerFD.cpp.

                                                          {
// ************************************************
 
  bool IsSubEventVarX = IsSubEventVarString(ProjectionVar_StrX);
  bool IsSubEventVarY = IsSubEventVarString(ProjectionVar_StrY);   
 
  int ProjectionVar_IntX, ProjectionVar_IntY;
  if (IsSubEventVarX) ProjectionVar_IntX = ReturnKinematicVectorFromString(ProjectionVar_StrX);
  else ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
  if (IsSubEventVarY) ProjectionVar_IntY = ReturnKinematicVectorFromString(ProjectionVar_StrY);
  else ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY); 
 
  //JM Loop over all events
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    const int EventSample = MCSamples[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 SampleHandlerFD::FillArray ( )

protected

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

Function which does the core reweighting, fills the SampleHandlerFD::SampleHandlerFD_array vector with the weight calculated from reweighting

Definition at line 374 of file SampleHandlerFD.cpp.

                                {
//************************************************
  //DB Reset which cuts to apply
  Selection = StoredSelection;
  
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
    ApplyShifts(iEvent);
    const EventInfo* _restrict_ MCEvent = &MCSamples[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) {
      SampleHandlerFD_array[GlobalBin] += totalweight;
      if (FirstTimeW2) SampleHandlerFD_array_w2[GlobalBin] += totalweight*totalweight;
    }
  }
}

FillHist()

void SampleHandlerFD::FillHist ( const int  Sample, TH1 *  Hist, 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 SampleHandlerFD_array

Definition at line 264 of file SampleHandlerFD.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]);
    }
  }
}

FindNominalBinAndEdges()

void SampleHandlerFD::FindNominalBinAndEdges ( )

protected

Definition at line 839 of file SampleHandlerFD.cpp.

                                             {
// ************************************************
  for (unsigned int event_i = 0; event_i < GetNEvents(); event_i++) {
    int Sample = MCSamples[event_i].NominalSample;
    const int dim = GetNDim(Sample);
    MCSamples[event_i].KinVar.resize(dim);
    MCSamples[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++) {
      MCSamples[event_i].KinVar[iDim] = GetPointerToKinematicParameter(GetKinVarName(Sample, iDim), event_i);
      if (std::isnan(*MCSamples[event_i].KinVar[iDim]) || std::isinf(*MCSamples[event_i].KinVar[iDim])) {
        MACH3LOG_ERROR("Variable {} for sample {} and dimension {} is ill-defined and equal to {}",
                       GetKinVarName(Sample, iDim), GetSampleTitle(Sample), dim, *MCSamples[event_i].KinVar[iDim]);
        throw MaCh3Exception(__FILE__, __LINE__);
      }
      const int bin = Binning->FindNominalBin(Sample, iDim, *MCSamples[event_i].KinVar[iDim]);
      int NBins_i = static_cast<int>(Binning->GetBinEdges(Sample, iDim).size() - 1);
      SetNominalBin(bin, NBins_i, MCSamples[event_i].NomBin[iDim]);
    }
  }
}

Get1DVarHist()

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

overridevirtual

Implements SampleHandlerBase.

Definition at line 1384 of file SampleHandlerFD.cpp.

                                                                                                                        {
// ************************************************
  //DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
  //   Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
  auto tmp_Selection = ApplyTemporarySelection(iSample, EventSelectionVec);
 
  //DB Define the histogram which will be returned
  TH1D* _h1DVar = nullptr;;
  if (Axis) {
    _h1DVar = new TH1D("","",Axis->GetNbins(),Axis->GetXbins()->GetArray());
  } else {
    std::vector<double> xBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_Str);
    _h1DVar = new TH1D("", "", int(xBinEdges.size())-1, xBinEdges.data());
  }
  _h1DVar->GetXaxis()->SetTitle(ProjectionVar_Str.c_str());
 
  if (IsSubEventVarString(ProjectionVar_Str)) {
    Fill1DSubEventHist(iSample, _h1DVar, ProjectionVar_Str, SubEventSelectionVec, WeightStyle);
  } else {
    //DB Grab the associated enum with the argument string
    int ProjectionVar_Int = ReturnKinematicParameterFromString(ProjectionVar_Str);
 
    //DB Loop over all events
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
      const int EventSample = MCSamples[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()

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

overridevirtual

Implements SampleHandlerBase.

Definition at line 1742 of file SampleHandlerFD.cpp.

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

Get2DVarHist()

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

overridevirtual

Implements SampleHandlerBase.

Definition at line 1456 of file SampleHandlerFD.cpp.

                                                                                          {
// ************************************************
  //DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
  //   Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
  auto tmp_Selection = ApplyTemporarySelection(iSample, EventSelectionVec);
 
  //DB Define the histogram which will be returned
  TH2D* _h2DVar = nullptr;
  if (AxisX && AxisY) {
    _h2DVar = new TH2D("","",AxisX->GetNbins(),AxisX->GetXbins()->GetArray(),AxisY->GetNbins(),AxisY->GetXbins()->GetArray());
  } else {
    std::vector<double> xBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_StrX);
    std::vector<double> yBinEdges = ReturnKinematicParameterBinning(iSample, ProjectionVar_StrY);
    _h2DVar = new TH2D("", "", int(xBinEdges.size())-1, xBinEdges.data(), int(yBinEdges.size())-1, yBinEdges.data());
  }
  _h2DVar->GetXaxis()->SetTitle(ProjectionVar_StrX.c_str());
  _h2DVar->GetYaxis()->SetTitle(ProjectionVar_StrY.c_str());
 
  bool IsSubEventHist = IsSubEventVarString(ProjectionVar_StrX) || IsSubEventVarString(ProjectionVar_StrY);
  if (IsSubEventHist) Fill2DSubEventHist(iSample, _h2DVar, ProjectionVar_StrX, ProjectionVar_StrY, SubEventSelectionVec, WeightStyle);
  else {
    //DB Grab the associated enum with the argument string
    int ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
    int ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY);
 
    //DB Loop over all events
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
      const int EventSample = MCSamples[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()

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

overridevirtual

Implements SampleHandlerBase.

Definition at line 1750 of file SampleHandlerFD.cpp.

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

GetArrayForSample()

std::vector< double > SampleHandlerFD::GetArrayForSample	(	const int	Sample,
		const double *	array
	)		const

Return a sub-array for a given sample.

Definition at line 2055 of file SampleHandlerFD.cpp.

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

GetBinningHandler()

const BinningHandler* SampleHandlerFD::GetBinningHandler ( ) const

inline

Get pointer to binning handler.

Definition at line 45 of file SampleHandlerFD.h.

{return Binning.get();}

GetDataArray() [1/2]

std::vector<double> SampleHandlerFD::GetDataArray ( ) const

inline

Return array storing data entries for every bin.

Definition at line 162 of file SampleHandlerFD.h.

                                         {
    return std::vector<double>(SampleHandlerFD_data, SampleHandlerFD_data + Binning->GetNBins());
  }

GetDataArray() [2/2]

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

inline

Return array storing data entries for every bin.

Definition at line 177 of file SampleHandlerFD.h.

                                                         {
    return GetArrayForSample(Sample, SampleHandlerFD_data);
  }

GetDataHist() [1/2]

TH1 * SampleHandlerFD::GetDataHist ( const int  Sample )

overridevirtual

Get Data histogram.

Implements SampleHandlerBase.

Definition at line 919 of file SampleHandlerFD.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]

TH1 * SampleHandlerFD::GetDataHist

(

const std::string &

Sample

)

Definition at line 929 of file SampleHandlerFD.cpp.

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

GetEventWeight()

M3::float_t SampleHandlerFD::GetEventWeight

(

const int

iEntry

)

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

                                                          {
// ************************************************
 
  // Virtual by default does nothing, has to happen before CalcWeightTotal
  CalcWeightFunc(iEntry);
 
  const EventInfo* _restrict_ MCEvent = &MCSamples[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 SampleHandlerFD::GetFinalPDGFromFileName ( const std::string &  FileName ) const

inlineprotected

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

Definition at line 397 of file SampleHandlerFD.h.

{return FileToFinalPDGMap.at(FileName);}

GetFlavourName()

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

inlineoverridevirtual

Implements SampleHandlerBase.

Definition at line 95 of file SampleHandlerFD.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 SampleHandlerFD::GetInitPDGFromFileName ( const std::string &  FileName ) const

inlineprotected

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

Definition at line 395 of file SampleHandlerFD.h.

{return FileToInitPDGMap.at(FileName);}

GetKinVarName()

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

overridevirtual

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

Parameters

iSample	Sample index
Dimension	Dimension index

Implements SampleHandlerBase.

Definition at line 2045 of file SampleHandlerFD.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 SampleHandlerFD::GetLikelihood ( ) const

overridevirtual

DB Multi-threaded GetLikelihood.

Implements SampleHandlerBase.

Definition at line 1250 of file SampleHandlerFD.cpp.

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

GetMCArray() [1/2]

std::vector<double> SampleHandlerFD::GetMCArray ( ) const

inline

Return array storing MC entries for every bin.

Definition at line 166 of file SampleHandlerFD.h.

                                       {
    return std::vector<double>(SampleHandlerFD_array, SampleHandlerFD_array + Binning->GetNBins());
  }

GetMCArray() [2/2]

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

inline

Return array storing MC entries for every bin.

Definition at line 181 of file SampleHandlerFD.h.

                                                       {
    return GetArrayForSample(Sample, SampleHandlerFD_array);
  }

GetMCHist() [1/2]

TH1 * SampleHandlerFD::GetMCHist ( const int  Sample )

overridevirtual

Get MC histogram.

Implements SampleHandlerBase.

Definition at line 901 of file SampleHandlerFD.cpp.

                                                {
// ************************************************
  FillHist(Sample, SampleDetails[Sample].MCHist, SampleHandlerFD_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]

TH1 * SampleHandlerFD::GetMCHist

(

const std::string &

Sample

)

Definition at line 912 of file SampleHandlerFD.cpp.

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

GetModeHist1D()

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

inline

Definition at line 128 of file SampleHandlerFD.h.

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

GetModeHist2D()

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

inline

Definition at line 131 of file SampleHandlerFD.h.

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

GetName()

std::string SampleHandlerFD::GetName ( ) const

overridevirtual

Implements SampleHandlerBase.

Definition at line 1138 of file SampleHandlerFD.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 SampleHandlerFD::GetNDim ( const int  Sample ) const

inlineoverridevirtual

DB Function to differentiate 1D or 2D binning.

Implements SampleHandlerBase.

Definition at line 33 of file SampleHandlerFD.h.

{ return SampleDetails[Sample].nDimensions; }

GetNOscChannels()

int SampleHandlerFD::GetNOscChannels ( const int  iSample ) const

inlineoverridevirtual

Implements SampleHandlerBase.

Definition at line 93 of file SampleHandlerFD.h.

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

GetNuOscillatorPointers()

const M3::float_t * SampleHandlerFD::GetNuOscillatorPointers

(

const int

iEvent

)

const

Definition at line 1098 of file SampleHandlerFD.cpp.

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

GetPointerToKinematicParameter() [1/2]

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

protectedpure virtual

Implemented in PySampleHandlerFD.

GetPointerToKinematicParameter() [2/2]

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

protectedpure virtual

Implemented in PySampleHandlerFD.

GetPointerToOscChannel()

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

protected

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

Definition at line 1983 of file SampleHandlerFD.cpp.

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

GetSampleIndex()

int SampleHandlerFD::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 871 of file SampleHandlerFD.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 SampleHandlerFD::GetSampleLikelihood ( const int  isample ) const

overridevirtual

Get likelihood for single sample.

Implements SampleHandlerBase.

Definition at line 1228 of file SampleHandlerFD.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)
  {
    const double DataVal = SampleHandlerFD_data[idx];
    const double MCPred = SampleHandlerFD_array[idx];
    const double w2 = SampleHandlerFD_array_w2[idx];
 
    //KS: Calculate likelihood using Barlow-Beeston Poisson or even IceCube
    negLogL += GetTestStatLLH(DataVal, MCPred, w2);
  }
  return negLogL;
}

GetSampleTitle()

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

inlineoverridevirtual

Implements SampleHandlerBase.

Definition at line 35 of file SampleHandlerFD.h.

{return SampleDetails[Sample].SampleTitle;}

GetW2Array() [1/2]

std::vector<double> SampleHandlerFD::GetW2Array ( ) const

inline

Return array storing W2 entries for every bin.

Definition at line 170 of file SampleHandlerFD.h.

                                       {
    return std::vector<double>(SampleHandlerFD_array_w2, SampleHandlerFD_array_w2 + Binning->GetNBins());
  }

GetW2Array() [2/2]

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

inline

Return array storing W2 entries for single sample.

Definition at line 185 of file SampleHandlerFD.h.

                                                       {
    return GetArrayForSample(Sample, SampleHandlerFD_array_w2);
  }

GetW2Hist() [1/2]

TH1 * SampleHandlerFD::GetW2Hist ( const int  Sample )

overridevirtual

Get W2 histogram.

Implements SampleHandlerBase.

Definition at line 883 of file SampleHandlerFD.cpp.

                                                {
// ************************************************
  FillHist(Sample, SampleDetails[Sample].W2Hist, SampleHandlerFD_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]

TH1 * SampleHandlerFD::GetW2Hist

(

const std::string &

Sample

)

Definition at line 894 of file SampleHandlerFD.cpp.

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

GetXBinVarName()

std::string SampleHandlerFD::GetXBinVarName ( const int  Sample ) const

inline

Definition at line 42 of file SampleHandlerFD.h.

{return GetKinVarName(Sample, 0);}

GetYBinVarName()

std::string SampleHandlerFD::GetYBinVarName ( const int  Sample ) const

inline

Definition at line 43 of file SampleHandlerFD.h.

{return GetKinVarName(Sample, 1);}

Init()

virtual void SampleHandlerFD::Init ( )

protectedpure virtual

Initialise any variables that your experiment specific SampleHandler needs.

Implemented in PySampleHandlerFD.

Initialise()

void SampleHandlerFD::Initialise ( )

protected

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

Definition at line 163 of file SampleHandlerFD.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();
  MCSamples.resize(nEvents);
  SetupFDMC();
 
  MACH3LOG_INFO("=============================================");
  MACH3LOG_INFO("Total number of events is: {}", GetNEvents());
 
  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__);
    }
  }
  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("=======================================================");
}

InitialiseNuOscillatorObjects()

void SampleHandlerFD::InitialiseNuOscillatorObjects

(

)

including Dan's magic NuOscillator

Definition at line 1020 of file SampleHandlerFD.cpp.

                                                    {
// ************************************************
  auto NuOscillatorConfigFile = Get<std::string>(SampleManager->raw()["NuOsc"]["NuOscConfigFile"], __FILE__ , __LINE__);
  auto EqualBinningPerOscChannel = Get<bool>(SampleManager->raw()["NuOsc"]["EqualBinningPerOscChannel"], __FILE__ , __LINE__);
 
  // TN override the sample setting if not using binned oscillation
  if (EqualBinningPerOscChannel) {
    if (YAML::LoadFile(NuOscillatorConfigFile)["General"]["CalculationType"].as<std::string>() == "Unbinned") {
      MACH3LOG_WARN("Tried using EqualBinningPerOscChannel while using Unbinned oscillation calculation, changing EqualBinningPerOscChannel to false");
      EqualBinningPerOscChannel = false;
    }
  }
  std::vector<const double*> OscParams = ParHandler->GetOscParsFromSampleName(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;
 
        for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
          if(MCSamples[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[MCSamples[iEvent].NominalSample].OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
          //DB Remove NC events from the arrays which are handed to the NuOscillator objects
          if (!MCSamples[iEvent].isNC && Channel == iChannel) {
            EnergyArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_etru)));
          }
        }
        std::sort(EnergyArray.begin(),EnergyArray.end());
 
        //============================================================================
        //DB Atmospheric only part, can only happen if truecz has been initialised within the experiment specific code
        if (*(MCSamples[0].rw_truecz) != M3::_BAD_DOUBLE_) {
          for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
            if(MCSamples[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[MCSamples[iEvent].NominalSample].OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
            //DB Remove NC events from the arrays which are handed to the NuOscillator objects
            if (!MCSamples[iEvent].isNC && Channel == iChannel) {
              CosineZArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_truecz)));
            }
          }
          std::sort(CosineZArray.begin(),CosineZArray.end());
        }
        Oscillator->SetOscillatorBinning(iSample, iChannel, EnergyArray, CosineZArray);
      } // end loop over channels
    } // end loop over samples
  }
}

InitialiseSplineObject()

void SampleHandlerFD::InitialiseSplineObject ( )

protected

Definition at line 1314 of file SampleHandlerFD.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<SMonolith*>(SplineHandler.get())) {
    (void) UnbinnedSpline;
    SetSplinePointers();
  } else {
    MACH3LOG_ERROR("Unsupported spline type encountered.");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

IsEventSelected()

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

protected

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

Definition at line 303 of file SampleHandlerFD.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 SampleHandlerFD::IsSubEventSelected ( const std::vector< KinematicCut > &  SubEventCuts, const int  iEvent, unsigned const int  iSubEvent, size_t  nsubevents  )

protected

JM Function which determines if a subevent is selected.

Definition at line 319 of file SampleHandlerFD.cpp.

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

IsSubEventVarString()

bool SampleHandlerFD::IsSubEventVarString

(

const std::string &

VarStr

)

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

Definition at line 1676 of file SampleHandlerFD.cpp.

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

LoadSingleSample()

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

Add new sample

Definition at line 113 of file SampleHandlerFD.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 mtupleprefix  = Get<std::string>(SampleSettings["InputFiles"]["mtupleprefix"], __FILE__, __LINE__);
  auto mtuplesuffix  = 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["SubSamples"].size());
  SingleSample.OscChannels.reserve(NChannels);
 
  int OscChannelCounter = 0;
  for (auto const &osc_channel : SampleSettings["SubSamples"]) {
    std::string MTupleFileName = mtupleprefix+osc_channel["mtuplefile"].as<std::string>()+mtuplesuffix;
 
    OscChannelInfo OscInfo;
    OscInfo.flavourName       = osc_channel["Name"].as<std::string>();
    OscInfo.flavourName_Latex = osc_channel["LatexName"].as<std::string>();
    OscInfo.InitPDG           = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
    OscInfo.FinalPDG          = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
    OscInfo.ChannelIndex      = OscChannelCounter;
 
    SingleSample.OscChannels.push_back(std::move(OscInfo));
 
    FileToInitPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
    FileToFinalPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
 
    SingleSample.mc_files.push_back(MTupleFileName);
    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);
  }
  SampleDetails[iSample] = std::move(SingleSample);
}

PassesSelection()

template<typename ParT >

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

protected

Definition at line 2065 of file SampleHandlerFD.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 SampleHandlerFD::PrepFunctionalParameters ( )

inlineprotectedvirtual

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

Definition at line 241 of file SampleHandlerFD.h.

{};

PrintIntegral()

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

Definition at line 1757 of file SampleHandlerFD.cpp.

                                                                                                                                             {
  constexpr int space = 14;
 
  bool printToFile=false;
  if (OutputFileName.CompareTo("/dev/null")) {printToFile = true;}
 
  bool printToCSV=false;
  if(OutputCSVFileName.CompareTo("/dev/null")) printToCSV=true;
 
  std::ofstream outfile;
  if (printToFile) {
    outfile.open(OutputFileName.Data(), std::ios_base::app);
    outfile.precision(7);
  }
 
  std::ofstream outcsv;
  if(printToCSV){
    outcsv.open(OutputCSVFileName, std::ios_base::app); // Appened to CSV
    outcsv.precision(7);
  }
 
  double PDFIntegral = 0.;
 
  std::vector< std::vector< TH1* > > IntegralList;
  IntegralList.resize(Modes->GetNModes());
 
  std::vector<double> ChannelIntegral;
  ChannelIntegral.resize(GetNOscChannels(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, GetXBinVarName(iSample),1,i,WeightStyle);
    } else {
      IntegralList[i] = CastVector<TH2, TH1>(ReturnHistsBySelection2D(iSample, GetXBinVarName(iSample), GetYBinVarName(iSample),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();
  }
  // KS: Clean memory we could use smart pointers in future
  CleanContainer(IntegralList);
}

PrintRates()

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

overridevirtual

Helper function to print rates for the samples with LLH.

Parameters

DataOnly

whether to print data only rates

Implements SampleHandlerBase.

Definition at line 1993 of file SampleHandlerFD.cpp.

                                                    {
// ***************************************************************************
  if (SampleHandlerFD_data == nullptr) {
    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 SampleHandlerFD::ReadConfig

(

)

Definition at line 59 of file SampleHandlerFD.cpp.

{
  auto ModeName = Get<std::string>(SampleManager->raw()["MaCh3ModeConfig"], __FILE__ , __LINE__);
  Modes = std::make_unique<MaCh3Modes>(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 SampleHandlerFD::RegisterFunctionalParameters ( )

protectedpure virtual

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

Implemented in PySampleHandlerFD.

RegisterIndividualFunctionalParameter()

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

protected

HH - a helper function for RegisterFunctionalParameter.

Definition at line 483 of file SampleHandlerFD.cpp.

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

ResetHistograms()

void SampleHandlerFD::ResetHistograms ( )

protected

Helper function to reset histograms.

Definition at line 472 of file SampleHandlerFD.cpp.

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

ResetShifts()

virtual void SampleHandlerFD::ResetShifts ( const int  iEvent )

inlineprotectedvirtual

HH - reset the shifted values to the original values.

Definition at line 250 of file SampleHandlerFD.h.

{(void)iEvent;};

ReturnHistsBySelection1D()

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

Definition at line 1900 of file SampleHandlerFD.cpp.

                                                                                                                           {
// ************************************************
  std::vector<TH1*> hHistList;
  std::string legendEntry;
 
  if (THStackLeg != nullptr) {delete THStackLeg;}
  THStackLeg = new TLegend(0.1,0.1,0.9,0.9);
 
  int iMax = -1;
  if (Selection1 == FDPlotType::kModePlot) {
    iMax = Modes->GetNModes();
  }
  if (Selection1 == FDPlotType::kOscChannelPlot) {
    iMax = GetNOscChannels(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,XAxis));
      THStackLeg->AddEntry(hHistList[i],(Modes->GetMaCh3ModeName(i)+Form(" : (%4.2f)",hHistList[i]->Integral())).c_str(),"f");
 
      hHistList[i]->SetFillColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
      hHistList[i]->SetLineColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
    }
    if (Selection1 == FDPlotType::kOscChannelPlot) {
      hHistList.push_back(Get1DVarHistByModeAndChannel(iSample, KinematicProjection,Selection2,i,WeightStyle,XAxis));
      THStackLeg->AddEntry(hHistList[i],(GetFlavourName(iSample, i)+Form(" | %4.2f",hHistList[i]->Integral())).c_str(),"f");
    }
  }
 
  return hHistList;
}

ReturnHistsBySelection2D()

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

Definition at line 1938 of file SampleHandlerFD.cpp.

                                                                                        {
// ************************************************
  std::vector<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,XAxis,YAxis));
    }
    if (Selection1 == FDPlotType::kOscChannelPlot) {
      hHistList.push_back(Get2DVarHistByModeAndChannel(iSample, KinematicProjectionX,KinematicProjectionY,Selection2,i,WeightStyle,XAxis,YAxis));
    }
  }
 
  return hHistList;
}

ReturnKinematicParameter() [1/2]

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

protectedpure virtual

Implemented in PySampleHandlerFD.

ReturnKinematicParameter() [2/2]

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

protectedpure virtual

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

Implemented in PySampleHandlerFD.

ReturnKinematicParameterBinning()

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

overrideprotectedvirtual

Return the binning used to draw a kinematic parameter.

Todo:

might be useful to allow overwriting this

Implements SampleHandlerBase.

Definition at line 1622 of file SampleHandlerFD.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 SampleHandlerFD::ReturnKinematicParameterFromString

(

const std::string &

KinematicStr

)

const

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

Definition at line 1569 of file SampleHandlerFD.cpp.

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

ReturnKinematicVector() [1/2]

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

inlineprotectedvirtual

Definition at line 292 of file SampleHandlerFD.h.

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

ReturnKinematicVector() [2/2]

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

inlineprotectedvirtual

Definition at line 291 of file SampleHandlerFD.h.

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

ReturnKinematicVectorFromString()

int SampleHandlerFD::ReturnKinematicVectorFromString

(

const std::string &

KinematicStr

)

const

Definition at line 1596 of file SampleHandlerFD.cpp.

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

ReturnStackedHistBySelection1D()

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

Definition at line 1970 of file SampleHandlerFD.cpp.

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

ReturnStackHistLegend()

TLegend* SampleHandlerFD::ReturnStackHistLegend ( )

inline

Definition at line 144 of file SampleHandlerFD.h.

{return THStackLeg;}

ReturnStringFromKinematicParameter()

std::string SampleHandlerFD::ReturnStringFromKinematicParameter

(

const int

KinematicVariable

)

const

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

Definition at line 1581 of file SampleHandlerFD.cpp.

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

ReturnStringFromKinematicVector()

std::string SampleHandlerFD::ReturnStringFromKinematicVector

(

const int

KinematicVariable

)

const

Definition at line 1608 of file SampleHandlerFD.cpp.

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

Reweight()

void SampleHandlerFD::Reweight ( )

overridevirtual

Implements SampleHandlerBase.

Definition at line 338 of file SampleHandlerFD.cpp.

                               {
//************************************************
  //KS: Reset the histograms before reweight 
  ResetHistograms();
  
  //You only need to do these things if Oscillator has been initialised
  //if not then you're not considering oscillations
  if (Oscillator) Oscillator->Evaluate();
 
  // Calculate weight coming from all splines if we initialised handler
  if(SplineHandler) SplineHandler->Evaluate();
 
  // 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 SampleHandlerFD::SaveAdditionalInfo ( TDirectory *  Dir )

overridevirtual

Store additional info in a chan.

Reimplemented from SampleHandlerBase.

Definition at line 1269 of file SampleHandlerFD.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<TH1D*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      data_hist->GetXaxis()->SetTitle(GetXBinVarName(iSample).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<TH2D*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      } else {
        data_hist = M3::Clone<TH2Poly>(dynamic_cast<TH2Poly*>(GetDataHist(iSample)), "data_" + GetSampleTitle(iSample));
      }
      data_hist->GetXaxis()->SetTitle(GetXBinVarName(iSample).c_str());
      data_hist->GetYaxis()->SetTitle(GetYBinVarName(iSample).c_str());
      data_hist->GetZaxis()->SetTitle("Number of Events");
    } else {
      data_hist = M3::Clone<TH1D>(dynamic_cast<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 SampleHandlerFD::SetBinning ( )

protected

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

Definition at line 738 of file SampleHandlerFD.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 SampleHandlerFD::SetSplinePointers ( )

protected

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

Todo:

Fix this mess :(

Definition at line 1169 of file SampleHandlerFD.cpp.

                                        {
// ************************************************
  //Now loop over events and get the spline bin for each event
  if (auto BinnedSpline = dynamic_cast<BinnedSplineHandler*>(SplineHandler.get())) {
    bool ThrowCrititcal = true;
    for (unsigned int j = 0; j < GetNEvents(); ++j) {
      const int SampleIndex = MCSamples[j].NominalSample;
      const auto SampleTitle = GetSampleTitle(SampleIndex);
      const int OscIndex = GetOscChannel(SampleDetails[SampleIndex].OscChannels, (*MCSamples[j].nupdgUnosc), (*MCSamples[j].nupdg));
      const int Mode = int(*(MCSamples[j].mode));
      const double Etrue = *(MCSamples[j].rw_etru);
      std::vector< std::vector<int> > EventSplines;
      switch(GetNDim(SampleIndex)) {
        case 1:
          EventSplines = BinnedSpline->GetEventSplines(SampleTitle, OscIndex, Mode, Etrue, *(MCSamples[j].KinVar[0]), 0.);
          break;
        case 2:
          EventSplines = BinnedSpline->GetEventSplines(SampleTitle, OscIndex, Mode, Etrue, *(MCSamples[j].KinVar[0]), *(MCSamples[j].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, *(MCSamples[j].KinVar[0]), *(MCSamples[j].KinVar[1]));
          break;
      }
      const int NSplines = static_cast<int>(EventSplines.size());
      if(NSplines == 0) continue;
      const int PointersBefore = static_cast<int>(MCSamples[j].total_weight_pointers.size());
      MCSamples[j].total_weight_pointers.resize(PointersBefore + NSplines);
 
      for(int spline = 0; spline < NSplines; spline++) {
        //Event Splines indexed as: sample name, oscillation channel, syst, mode, etrue, var1, var2 (var2 is a dummy 0 for 1D splines)
        MCSamples[j].total_weight_pointers[PointersBefore+spline] = 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
    } // end loop over events
  } else if (auto UnbinnedSpline = dynamic_cast<SMonolith*>(SplineHandler.get())) {
    #ifdef _LOW_MEMORY_STRUCTS_
    for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) {
      MCSamples[iEvent].total_weight_pointers.push_back(UnbinnedSpline->retPointer(iEvent));
    }
    #else
    (void) UnbinnedSpline;
    MACH3LOG_ERROR("Unbinned splines only for float :(");
    throw MaCh3Exception(__FILE__, __LINE__);
    #endif
  } else {
    MACH3LOG_ERROR("Not supported splines");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
}

SetupExperimentMC()

virtual int SampleHandlerFD::SetupExperimentMC ( )

protectedpure virtual

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

Implemented in PySampleHandlerFD.

SetupFDMC()

virtual void SampleHandlerFD::SetupFDMC ( )

protectedpure virtual

Function which translates experiment struct into core struct.

Implemented in PySampleHandlerFD.

SetupFunctionalParameters()

void SampleHandlerFD::SetupFunctionalParameters ( )

protectedvirtual

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

Todo:

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

Definition at line 503 of file SampleHandlerFD.cpp.

                                                {
// **************************************************
  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());
  funcParsGrid.resize(GetNEvents());
 
  // 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
      bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
      if (!ModeMatch) {
        MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
        continue;
      }
      // Now check whether within kinematic bounds
      bool IsSelected = 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 SampleHandlerFD::SetupKinematicMap ( )

protected

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

Definition at line 227 of file SampleHandlerFD.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"};
  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__);
    }
  }
}

SetupNormParameters()

void SampleHandlerFD::SetupNormParameters ( )

protected

Setup the norm parameters by assigning each event with bin.

Definition at line 613 of file SampleHandlerFD.cpp.

                                          {
// ***************************************************************************
  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 MCSamples tree
  CalcNormsBins(norm_parameters, norms_bins);
 
  //DB Attempt at reducing impact of SystematicHandlerGeneric::calcReweight()
  int counter;
  for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) {
    counter = 0;
 
    MCSamples[iEvent].norm_pointers.resize(norms_bins[iEvent].size());
    for(auto const & norm_bin: norms_bins[iEvent]) {
      MCSamples[iEvent].norm_pointers[counter] = ParHandler->RetPointer(norm_bin);
      counter += 1;
    }
  }
}

SetupNuOscillatorPointers()

void SampleHandlerFD::SetupNuOscillatorPointers

(

)

Definition at line 1081 of file SampleHandlerFD.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) {
     MCSamples[iEvent].total_weight_pointers.push_back(osc_w_pointer);
   }
  }
}

SetupReweightArrays()

void SampleHandlerFD::SetupReweightArrays ( )

protected

Initialise data, MC and W2 histograms.

Definition at line 724 of file SampleHandlerFD.cpp.

                                          {
// ************************************************
  SampleHandlerFD_array = new double[Binning->GetNBins()];
  SampleHandlerFD_array_w2 = new double[Binning->GetNBins()];
  SampleHandlerFD_data = new double[Binning->GetNBins()];
 
  for (int i = 0; i < Binning->GetNBins(); ++i) {
    SampleHandlerFD_array[i] = 0.0;
    SampleHandlerFD_array_w2[i] = 0.0;
    SampleHandlerFD_data[i] = 0.0;
  }
}

SetupSplines()

virtual void SampleHandlerFD::SetupSplines ( )

protectedpure virtual

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

Todo:

abstract the spline initialisation completely to core

Implemented in PySampleHandlerFD.

Member Data Documentation

_modeNomWeightMap

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

protected

Definition at line 382 of file SampleHandlerFD.h.

Binning

std::unique_ptr<BinningHandler> SampleHandlerFD::Binning

protected

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

Definition at line 329 of file SampleHandlerFD.h.

FileToFinalPDGMap

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

private

Definition at line 401 of file SampleHandlerFD.h.

FileToInitPDGMap

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

private

Definition at line 400 of file SampleHandlerFD.h.

FirstTimeW2

bool SampleHandlerFD::FirstTimeW2

protected

KS:Super hacky to update W2 or not.

Definition at line 390 of file SampleHandlerFD.h.

funcParsFuncMap

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

protected

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

Definition at line 269 of file SampleHandlerFD.h.

funcParsGrid

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

protected

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

Definition at line 252 of file SampleHandlerFD.h.

funcParsMap

std::vector<FunctionalShifter> SampleHandlerFD::funcParsMap

protected

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

Definition at line 257 of file SampleHandlerFD.h.

funcParsNamesMap

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

protected

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

Definition at line 266 of file SampleHandlerFD.h.

funcParsNamesVec

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

protected

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

Definition at line 271 of file SampleHandlerFD.h.

funcParsVec

std::vector<FunctionalParameter> SampleHandlerFD::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 263 of file SampleHandlerFD.h.

KinematicParameters

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

protected

Mapping between string and kinematic enum.

Definition at line 369 of file SampleHandlerFD.h.

KinematicVectors

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

protected

Definition at line 374 of file SampleHandlerFD.h.

MCSamples

std::vector<EventInfo> SampleHandlerFD::MCSamples

protected

Stores information about every MC event.

Definition at line 340 of file SampleHandlerFD.h.

Oscillator

std::shared_ptr<OscillationHandler> SampleHandlerFD::Oscillator

protected

Contains oscillator handling calculating oscillation probabilities.

Definition at line 217 of file SampleHandlerFD.h.

ParHandler

ParameterHandlerGeneric* SampleHandlerFD::ParHandler = nullptr

protected

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

Definition at line 348 of file SampleHandlerFD.h.

ReversedKinematicParameters

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

protected

Mapping between kinematic enum and string.

Definition at line 371 of file SampleHandlerFD.h.

ReversedKinematicVectors

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

protected

Definition at line 375 of file SampleHandlerFD.h.

SampleDetails

std::vector<SampleInfo> SampleHandlerFD::SampleDetails

protected

Stores info about currently initialised sample.

Definition at line 342 of file SampleHandlerFD.h.

SampleHandlerFD_array

double* SampleHandlerFD::SampleHandlerFD_array

protected

DB Array to be filled after reweighting.

Definition at line 331 of file SampleHandlerFD.h.

SampleHandlerFD_array_w2

double* SampleHandlerFD::SampleHandlerFD_array_w2

protected

KS Array used for MC stat.

Definition at line 333 of file SampleHandlerFD.h.

SampleHandlerFD_data

double* SampleHandlerFD::SampleHandlerFD_data

protected

DB Array to be filled in AddData.

Definition at line 335 of file SampleHandlerFD.h.

SampleHandlerName

std::string SampleHandlerFD::SampleHandlerName

protected

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

Definition at line 352 of file SampleHandlerFD.h.

SampleManager

std::unique_ptr<Manager> SampleHandlerFD::SampleManager

protected

The manager object used to read the sample yaml file.

Definition at line 379 of file SampleHandlerFD.h.

Selection

std::vector< std::vector< KinematicCut > > SampleHandlerFD::Selection

protected

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

Definition at line 365 of file SampleHandlerFD.h.

SplineHandler

std::unique_ptr<SplineBase> SampleHandlerFD::SplineHandler

protected

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

Definition at line 214 of file SampleHandlerFD.h.

StoredSelection

std::vector< std::vector< KinematicCut > > SampleHandlerFD::StoredSelection

protected

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

Definition at line 362 of file SampleHandlerFD.h.

THStackLeg

TLegend* SampleHandlerFD::THStackLeg = nullptr

protected

DB Miscellaneous Variables.

Definition at line 386 of file SampleHandlerFD.h.

UpdateW2

bool SampleHandlerFD::UpdateW2

protected

KS:Super hacky to update W2 or not.

Definition at line 392 of file SampleHandlerFD.h.

---

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

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