SampleHandlerFD.h

Go to the documentation of this file.

#pragma once
 
//MaCh3 includes
#include "Splines/BinnedSplineHandler.h"
#include "Parameters/ParameterHandlerGeneric.h"
#include "Samples/SampleHandlerBase.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 SampleHandlerFD :  public SampleHandlerBase
{
 public:
  //######################################### Functions #########################################
  SampleHandlerFD(std::string ConfigFileName, ParameterHandlerGeneric* xsec_cov,
                  const std::shared_ptr<OscillationHandler>& OscillatorObj_ = nullptr);
  virtual ~SampleHandlerFD();
 
  int GetNDim(const int Sample) const override { return SampleDetails[Sample].nDimensions; }
  std::string GetName() const override;
  std::string GetSampleTitle(const int Sample) const override {return SampleDetails[Sample].SampleTitle;}
 
  std::string GetKinVarName(const int iSample, const int Dimension) const override;
 
  std::string GetXBinVarName(const int Sample) const {return GetKinVarName(Sample, 0);}
  std::string GetYBinVarName(const int Sample) const {return GetKinVarName(Sample, 1);}
  const BinningHandler* GetBinningHandler() const {return Binning.get();}
 
  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 SampleHandlerFDBase
  //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) override;
  double GetLikelihood() const override;
  double GetSampleLikelihood(const int isample) const override;
  //===============================================================================
 
  int GetSampleIndex(const std::string& SampleTitle) const;
 
  TH1* GetDataHist(const int Sample) override;
  TH1* GetDataHist(const std::string& Sample);
 
  TH1* GetMCHist(const int Sample) override;
  TH1* GetMCHist(const std::string& Sample);
 
  TH1* GetW2Hist(const int Sample) override;
  TH1* GetW2Hist(const std::string& Sample);
 
  void Reweight() override;
  M3::float_t GetEventWeight(const int iEntry);
 
  void InitialiseNuOscillatorObjects();
  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 {return static_cast<int>(SampleDetails[iSample].OscChannels.size());};
 
  std::string GetFlavourName(const int iSample, const int iChannel) const override {
    if (iChannel < 0 || iChannel > GetNOscChannels(iSample)) {
      MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
      throw MaCh3Exception(__FILE__ , __LINE__);      
    }
    return SampleDetails[iSample].OscChannels[iChannel].flavourName;
  }
  std::vector<std::vector<KinematicCut>> ApplyTemporarySelection(const int iSample,
                                                                 const std::vector<KinematicCut>& ExtraCuts);
  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) {
    return Get1DVarHistByModeAndChannel(iSample, GetXBinVarName(iSample), m, s, style);
  }
  TH2 *GetModeHist2D(const int iSample, int s, int m, int style = 0) {
    return Get2DVarHistByModeAndChannel(iSample, GetXBinVarName(iSample),GetYBinVarName(iSample), m, s, style);
  }
 
  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() {return THStackLeg;}
  
  int ReturnKinematicParameterFromString(const std::string& KinematicStr) const;
  std::string ReturnStringFromKinematicParameter(const int KinematicVariable) const;
 
  void SaveAdditionalInfo(TDirectory* Dir) override;
 
  // === JM declare the same functions for kinematic vectors ===
  int ReturnKinematicVectorFromString(const std::string& KinematicStr) const;
  std::string ReturnStringFromKinematicVector(const int KinematicVariable) const;
  // ===========================================================
  bool IsSubEventVarString(const std::string& VarStr);
 
  std::vector<double> GetDataArray() const {
    return std::vector<double>(SampleHandlerFD_data, SampleHandlerFD_data + Binning->GetNBins());
  }
  std::vector<double> GetMCArray() const {
    return std::vector<double>(SampleHandlerFD_array, SampleHandlerFD_array + Binning->GetNBins());
  }
  std::vector<double> GetW2Array() const {
    return std::vector<double>(SampleHandlerFD_array_w2, SampleHandlerFD_array_w2 + Binning->GetNBins());
  }
  std::vector<double> GetArrayForSample(const int Sample, const double* array) const;
 
  std::vector<double> GetDataArray(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandlerFD_data);
  }
  std::vector<double> GetMCArray(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandlerFD_array);
  }
  std::vector<double> GetW2Array(const int Sample) const {
    return GetArrayForSample(Sample, SampleHandlerFD_array_w2);
  }
 
 protected:
  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 SetupFDMC() = 0;
 
  void Initialise();
  
  std::unique_ptr<SplineBase> SplineHandler;
 
  std::shared_ptr<OscillationHandler> Oscillator;
  //===============================================================================
  void SetSplinePointers();
 
  //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;};
  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;
 
  virtual void CalcWeightFunc(int iEvent) {return; (void)iEvent;};
 
  virtual double ReturnKinematicParameter(std::string KinematicParamter, int iEvent) = 0;
  virtual double ReturnKinematicParameter(int KinematicVariable, int iEvent) = 0;
  
  // === JM declare the same functions for kinematic vectors ===
  virtual std::vector<double> ReturnKinematicVector(std::string KinematicParameter, int iEvent) {return {}; (void)KinematicParameter; (void)iEvent;};
  virtual std::vector<double> ReturnKinematicVector(int KinematicVariable, int iEvent) {return {}; (void)KinematicVariable; (void)iEvent;};
  // ===========================================================
 
  std::vector<double> ReturnKinematicParameterBinning(const int Sample, const std::string &KinematicParameter) const override;
 
  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;
  void SetupNormParameters();
 
  //===============================================================================
  void FillHist(const int Sample, TH1* Hist, double* Array);
 
#ifdef MULTITHREAD
  void FillArray_MP();
#endif
  void FillArray();
 
  void ResetHistograms();
  
  //===============================================================================
  //DB Variables required for GetLikelihood
  std::unique_ptr<BinningHandler> Binning;
  double* SampleHandlerFD_array;
  double* SampleHandlerFD_array_w2;
  double* SampleHandlerFD_data;
  //===============================================================================
 
  //===============================================================================
  std::vector<EventInfo> MCSamples;
  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::unordered_map<std::string, NuPDG> FileToInitPDGMap;
  std::unordered_map<std::string, NuPDG> FileToFinalPDGMap;
  
  enum FDPlotType {
    kModePlot = 0,
    kOscChannelPlot = 1
  };
};