SampleHandlerBase.h

Go to the documentation of this file.

#pragma once
 
//MaCh3 includes
#include "Splines/BinnedSplineHandler.h"
#include "Parameters/ParameterHandlerGeneric.h"
#include "Samples/SampleHandlerInterface.h"
#include "Samples/OscillationHandler.h"
#include "Samples/FarDetectorCoreInfoStruct.h"
#include "Samples/BinningHandler.h"
 
_MaCh3_Safe_Include_Start_ //{
#include "THStack.h"
#include "TLegend.h"
_MaCh3_Safe_Include_End_ //}
 
class SampleHandlerBase :  public SampleHandlerInterface
{
 public:
  //######################################### Functions #########################################
  SampleHandlerBase(std::string ConfigFileName, ParameterHandlerGeneric* xsec_cov,
                  const std::shared_ptr<OscillationHandler>& OscillatorObj_ = nullptr);
  virtual ~SampleHandlerBase();
 
  int GetNDim(const int Sample) const final { return SampleDetails[Sample].nDimensions; }
  std::string GetName() const final;
  std::string GetSampleTitle(const int Sample) const final {return SampleDetails[Sample].SampleTitle;}
 
  std::string GetKinVarName(const int iSample, const int Dimension) const final;
 
  void PrintIntegral(const int iSample, const TString& OutputName="/dev/null", const int WeightStyle=0, const TString& OutputCSVName="/dev/null");
 
  //===============================================================================
  // DB Reweighting and Likelihood functions
 
  //ETA - abstract these to SampleHandlerBase
  //DB Require these four functions to allow conversion from TH1(2)D to array for multi-threaded GetLikelihood
  void AddData(const int Sample, TH1* Data);
  void AddData(const int Sample, const std::vector<double>& Data_Array);
 
  void PrintRates(const bool DataOnly = false) final;
  double GetLikelihood() const override;
  double GetSampleLikelihood(const int isample) const override;
  //===============================================================================
 
  int GetSampleIndex(const std::string& SampleTitle) const;
 
  const TH1* GetDataHist(const int Sample) final;
  const TH1* GetDataHist(const std::string& Sample);
 
  const TH1* GetMCHist(const int Sample) final;
  const TH1* GetMCHist(const std::string& Sample);
 
  const TH1* GetW2Hist(const int Sample) final;
  const TH1* GetW2Hist(const std::string& Sample);
  void Reweight() override;
  M3::float_t GetEventWeight(const int iEntry);
 
  const M3::float_t* GetNuOscillatorPointers(const int iEvent) const;
 
  int GetNOscChannels(const int iSample) const final {return static_cast<int>(SampleDetails[iSample].OscChannels.size());};
 
  std::string GetFlavourName(const int iSample, const int iChannel) const final {
    if (iChannel < 0 || iChannel > GetNOscChannels(iSample)) {
      MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
      throw MaCh3Exception(__FILE__ , __LINE__);
    }
    return SampleDetails[iSample].OscChannels[iChannel].flavourName;
  }
  std::unique_ptr<TH1> Get1DVarHist(const int iSample, const std::string &ProjectionVar,
                                    const std::vector<KinematicCut> &EventSelectionVec = {}, int WeightStyle = 0,
                                    const std::vector<KinematicCut> &SubEventSelectionVec = {}) final;
  std::unique_ptr<TH2> Get2DVarHist(const int iSample, const std::string& ProjectionVarX, const std::string& ProjectionVarY,
                                    const std::vector< KinematicCut >& EventSelectionVec = {},
                                    int WeightStyle = 0, const std::vector< KinematicCut >& SubEventSelectionVec = {}) final;
  std::vector<KinematicCut> BuildModeChannelSelection(const int iSample, const int kModeToFill, const int kChannelToFill) const;
 
  void Fill1DSubEventHist(const int iSample, TH1D* _h1DVar, const std::string& ProjectionVar,
                          const std::vector< KinematicCut >& SubEventSelectionVec = {},
                          int WeightStyle=0);
  void Fill2DSubEventHist(const int iSample, TH2* _h2DVar, const std::string& ProjectionVarX, const std::string& ProjectionVarY,
                          const std::vector< KinematicCut >& SubEventSelectionVec = {}, int WeightStyle = 0);
 
  std::unique_ptr<TH1> Get1DVarHistByModeAndChannel(const int iSample, const std::string& ProjectionVar_Str,
                                                    const int kModeToFill = -1, const int kChannelToFill = -1,
                                                    const int WeightStyle = 0) final;
  std::unique_ptr<TH2> Get2DVarHistByModeAndChannel(const int iSample, const std::string& ProjectionVar_StrX,
                                                    const std::string& ProjectionVar_StrY, const int kModeToFill = -1,
                                                    const int kChannelToFill = -1, const int WeightStyle = 0) final;
 
  std::unique_ptr<TH1> GetModeHist1D(const int iSample, int s, int m, int style = 0) {
    return Get1DVarHistByModeAndChannel(iSample, GetKinVarName(iSample, 0), m, s, style);
  }
  std::unique_ptr<TH2> GetModeHist2D(const int iSample, int s, int m, int style = 0) {
    return Get2DVarHistByModeAndChannel(iSample, GetKinVarName(iSample, 0), GetKinVarName(iSample, 1), m, s, style);
  }
 
  std::vector<std::unique_ptr<TH1>> ReturnHistsBySelection1D(const int iSample, const std::string& KinematicProjection,
                                                             const int Selection1, const int Selection2 = -1,
                                                             const int WeightStyle = 0);
  std::vector<std::unique_ptr<TH2>> ReturnHistsBySelection2D(const int iSample, const std::string& KinematicProjectionX,
                                                             const std::string& KinematicProjectionY,
                                                             const int Selection1, const int Selection2=-1,
                                                             const int WeightStyle=0);
  std::unique_ptr<THStack> ReturnStackedHistBySelection1D(const int iSample, const std::string& KinematicProjection,
                                          const int Selection1, const int Selection2 = -1, const int WeightStyle = 0);
  const TLegend* ReturnStackHistLegend() const {return THStackLeg;}
 
  int ReturnKinematicParameterFromString(const std::string& KinematicStr) const;
  std::string ReturnStringFromKinematicParameter(const int KinematicVariable) const;
 
  void SaveAdditionalInfo(TDirectory* Dir) final;
 
  int ReturnKinematicVectorFromString(const std::string& KinematicStr) const;
  std::string ReturnStringFromKinematicVector(const int KinematicVariable) const;
  bool IsSubEventVarString(const std::string& VarStr) const;
 
  auto GetDataArray() const {
    return SampleHandler_data;
  }
  auto GetMCArray() const {
    return SampleHandler_array;
  }
  auto GetW2Array() const {
    return SampleHandler_array_w2;
  }
  std::vector<double> GetArrayForSample(const int Sample, std::vector<double> const & array) const;
 
  std::vector<double> GetDataArray(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandler_data);
  }
  std::vector<double> GetMCArray(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandler_array);
  }
  std::vector<double> GetW2Array(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandler_array_w2);
  }
 
 protected:
  void InitialiseNuOscillatorObjects();
  void SetupNuOscillatorPointers();
  void ReadConfig();
  void LoadSingleSample(const int iSample, const YAML::Node& Settings);
 
  virtual void AddAdditionalWeightPointers() = 0;
 
  void SetupKinematicMap();
 
  virtual void SetupSplines() = 0;
 
  //DB Require all objects to have a function which reads in the MC
  virtual void Init() = 0;
 
  virtual int SetupExperimentMC() = 0;
 
  virtual void SetupMC() = 0;
  virtual void InititialiseData() = 0;
 
  void Initialise();
 
  std::unique_ptr<SplineBase> SplineHandler;
 
  std::shared_ptr<OscillationHandler> Oscillator;
  //===============================================================================
  void SetSplinePointers();
  std::vector< std::vector<int> > GetSplineBins(int Event, BinnedSplineHandler* BinnedSpline, bool& ThrowCrititcal) const;
 
  //Functions which find the nominal bin and bin edges
  void FindNominalBinAndEdges();
 
  void SetBinning();
 
  void SetupReweightArrays();
  //===============================================================================
 
  // ----- Functional Parameters -----
  virtual void SetupFunctionalParameters();
  void RegisterIndividualFunctionalParameter(const std::string& fpName, int fpEnum, FuncParFuncType fpFunc);
  virtual void RegisterFunctionalParameters() = 0;
  virtual void PrepFunctionalParameters(){};
  virtual void ApplyShifts(const int iEvent);
 
  bool IsEventSelected(const int iSample, const int iEvent) _noexcept_;
  bool IsSubEventSelected(const std::vector<KinematicCut> &SubEventCuts, const int iEvent, unsigned const int iSubEvent, size_t nsubevents);
  virtual void ResetShifts(const int iEvent) {(void)iEvent;};
  virtual void FinaliseShifts(const int iEvent) {(void)iEvent;};
  std::vector<std::vector<FunctionalShifter*>> funcParsGrid;
  std::vector<FunctionalShifter> funcParsMap;
 
 
  std::vector<FunctionalParameter> funcParsVec;
  std::unordered_map<std::string, int> funcParsNamesMap;
  std::unordered_map<int, FuncParFuncType> funcParsFuncMap;
  std::vector<std::string> funcParsNamesVec = {};
 
  void CalcNormsBins(std::vector<NormParameter>& norm_parameters, std::vector< std::vector< int > >& norms_bins);
  template <typename ParT> bool PassesSelection(const ParT& Par, std::size_t iEvent);
  M3::float_t CalcWeightTotal(const EventInfo* _restrict_ MCEvent) const _noexcept_;
 
  virtual void CalcWeightFunc(const int iEvent) {return; (void)iEvent;};
 
  double ReturnKinematicParameter(const std::string& KinematicParameter, int iEvent) const {
    return ReturnKinematicParameter(ReturnKinematicParameterFromString(KinematicParameter), iEvent);
  }
  virtual double ReturnKinematicParameter(const int KinematicVariable, const int iEvent) const = 0;
 
  // === JM declare the same functions for kinematic vectors ===
  std::vector<double> ReturnKinematicVector(const std::string& KinematicParameter, const int iEvent) const {
    return ReturnKinematicVector(ReturnKinematicVectorFromString(KinematicParameter), iEvent);
  }
  virtual std::vector<double> ReturnKinematicVector(const int KinematicVariable, const int iEvent) const {
    return {}; (void)KinematicVariable; (void)iEvent;};
  // ===========================================================
 
  std::vector<double> ReturnKinematicParameterBinning(const int Sample, const std::string &KinematicParameter) const final;
 
  const double* GetPointerToKinematicParameter(const std::string& KinematicParameter, int iEvent) const {
    return GetPointerToKinematicParameter(ReturnKinematicParameterFromString(KinematicParameter), iEvent);
  }
  virtual const double* GetPointerToKinematicParameter(const int KinematicVariable, const int iEvent) const = 0;
 
  const double* GetPointerToOscChannel(const int iEvent) const;
  void SetupNormParameters();
  void SetupOscParameters();
  //===============================================================================
  void FillHist(const int Sample, TH1* Hist, std::vector<double> &Array);
 
#ifdef MULTITHREAD
  void FillArray_MP();
#endif
  void FillArray();
 
  void ResetHistograms();
 
  //===============================================================================
  //DB Variables required for GetLikelihood
  std::unique_ptr<BinningHandler> Binning;
  std::vector<double> SampleHandler_array;
  std::vector<double> SampleHandler_array_w2;
  std::vector<double> SampleHandler_data;
  //===============================================================================
 
  //===============================================================================
  std::vector<EventInfo> MCEvents;
  std::vector<SampleInfo> SampleDetails;
  //===============================================================================
 
  //===============================================================================
  //DB Covariance Objects
  ParameterHandlerGeneric *ParHandler = nullptr;
 
  //===============================================================================
  std::string SampleHandlerName;
 
  //===========================================================================
  //DB Vectors to store which kinematic cuts we apply
  //like in XsecNorms but for events in sample. Read in from sample yaml file
  //What gets used in IsEventSelected, which gets set equal to user input plus
  //all the vectors in StoreSelection
 
  std::vector< std::vector< KinematicCut > > StoredSelection;
  std::vector< std::vector< KinematicCut > > Selection;
   //===========================================================================
 
  const std::unordered_map<std::string, int>* KinematicParameters;
  const std::unordered_map<int, std::string>* ReversedKinematicParameters;
 
  // === JM mapping between string and kinematic vector enum ===
  const std::unordered_map<std::string, int>* KinematicVectors;
  const std::unordered_map<int, std::string>* ReversedKinematicVectors;
  // ===========================================================
 
  std::unique_ptr<Manager> SampleManager;
  void InitialiseSplineObject();
 
  std::unordered_map<std::string, double> _modeNomWeightMap;
 
  //===============================================================================
  TLegend* THStackLeg = nullptr;
  //===============================================================================
 
  bool FirstTimeW2;
  bool UpdateW2;
 
  NuPDG GetInitPDGFromFileName(const std::string& FileName) const {return FileToInitPDGMap.at(FileName);}
  NuPDG GetFinalPDGFromFileName(const std::string& FileName) const {return FileToFinalPDGMap.at(FileName);}
 
 private:
   std::vector<std::vector<KinematicCut>> ApplyTemporarySelection(const int iSample,
                                                                  const std::vector<KinematicCut>& ExtraCuts);
 
  std::unordered_map<std::string, NuPDG> FileToInitPDGMap;
  std::unordered_map<std::string, NuPDG> FileToFinalPDGMap;
 
  enum FDPlotType {
    kModePlot = 0,
    kOscChannelPlot = 1
  };
};