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 "THStack.h"
#include "TLegend.h"
 
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 { return nDimensions; }
  std::string GetSampleName(int iSample = 0) const override;
  std::string GetTitle() const override {return SampleTitle;}
 
  std::string GetXBinVarName() {return XVarStr;}
  std::string GetYBinVarName() {return YVarStr;}
 
  void PrintIntegral(TString OutputName="/dev/null", int WeightStyle=0, 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(TH1D* Data);
  void AddData(TH2D* Data);
  void AddData(std::vector<double> &data);
  void AddData(std::vector< std::vector <double> > &data);
 
  double GetLikelihood() override;
  //===============================================================================
 
  TH1* GetMCHist(const int Dimension);
 
  TH1* GetW2Hist(const int Dimension);
 
  TH1* GetDataHist(const int Dimension);
 
  void Reweight() override;
  M3::float_t GetEventWeight(const int iEntry) const;
 
  void InitialiseNuOscillatorObjects();
  void SetupNuOscillatorPointers();
 
  void ReadSampleConfig();
 
  int GetNOscChannels() override {return static_cast<int>(OscChannels.size());}
 
  std::string GetFlavourName(const int iChannel) {
    if (iChannel < 0 || iChannel > static_cast<int>(OscChannels.size())) {
      MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
      throw MaCh3Exception(__FILE__ , __LINE__);      
    }
    return OscChannels[iChannel].flavourName;
  }
 
  TH1* Get1DVarHist(const std::string& ProjectionVar, const std::vector< KinematicCut >& EventSelectionVec = std::vector< KinematicCut >(),
                    int WeightStyle=0, TAxis* Axis=nullptr, const std::vector< KinematicCut >& SubEventSelectionVec = std::vector< KinematicCut >());
  TH2* Get2DVarHist(const std::string& ProjectionVarX, const std::string& ProjectionVarY,
                    const std::vector< KinematicCut >& EventSelectionVec = std::vector< KinematicCut >(),
                    int WeightStyle=0, TAxis* AxisX=nullptr, TAxis* AxisY=nullptr,
                    const std::vector< KinematicCut >& SubEventSelectionVec = std::vector< KinematicCut >());
  
  void Fill1DSubEventHist(TH1D* _h1DVar, const std::string& ProjectionVar, const std::vector< KinematicCut >& SubEventSelectionVec = std::vector< KinematicCut >(),
                    int WeightStyle=0);
  void Fill2DSubEventHist(TH2D* _h2DVar, const std::string& ProjectionVarX, const std::string& ProjectionVarY,
                    const std::vector< KinematicCut >& SubEventSelectionVec = std::vector< KinematicCut >(), int WeightStyle=0);
  
  TH1* Get1DVarHistByModeAndChannel(const std::string& ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis* Axis=nullptr);
  TH2* Get2DVarHistByModeAndChannel(const std::string& ProjectionVar_StrX, const std::string& ProjectionVar_StrY, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis* AxisX=nullptr, TAxis* AxisY=nullptr);
 
  TH1 *GetModeHist1D(int s, int m, int style = 0) {
    return Get1DVarHistByModeAndChannel(XVarStr,m,s,style);
  }
  TH2 *GetModeHist2D(int s, int m, int style = 0) {
    return Get2DVarHistByModeAndChannel(XVarStr,YVarStr,m,s,style);
  }
 
  std::vector<TH1*> ReturnHistsBySelection1D(std::string KinematicProjection, int Selection1,int Selection2=-1, int WeightStyle=0, TAxis* Axis=0);
  std::vector<TH2*> ReturnHistsBySelection2D(std::string KinematicProjectionX, std::string KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis* XAxis=0, TAxis* YAxis=0);
  THStack* ReturnStackedHistBySelection1D(std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis* Axis=0);
  TLegend* ReturnStackHistLegend() {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);
 
 protected:
  virtual void SetupWeightPointers() = 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<BinnedSplineHandler> SplineHandler;
 
  std::shared_ptr<OscillationHandler> Oscillator;
  //===============================================================================
  void FillSplineBins();
 
  //Functions which find the nominal bin and bin edges
  void FindNominalBinAndEdges1D();
  void FindNominalBinAndEdges2D();
 
  //DB Overrided functions of base class which calculate erec bin and boundaries for reweighting speedup in beam samples
  //ETA - this can be done using the core info stored in the new fdmc_struct
  void Set1DBinning(size_t nbins, double* boundaries);
  void Set2DBinning(size_t nbins1, double* boundaries1, size_t nbins2, double* boundaries2);
  void Set1DBinning(std::vector<double> &XVec) {Set1DBinning(XVec.size()-1, XVec.data());};
  void Set2DBinning(std::vector<double> &XVec, std::vector<double> &YVec) {Set2DBinning(XVec.size()-1, XVec.data(), YVec.size()-1, YVec.data());};
  void SetupSampleBinning();
  void SetupReweightArrays(const size_t numberXBins, const size_t numberYBins);
 
  std::string XVarStr, YVarStr;
  //===============================================================================
 
  // ----- 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(int iEvent);
 
  bool IsEventSelected(const int iEvent);
  bool IsSubEventSelected(const std::vector<KinematicCut> &SubEventCuts, const int iEvent, unsigned const int iSubEvent, size_t nsubevents);
  virtual void resetShifts(int iEvent) {(void)iEvent;};
  std::vector<FunctionalParameter> funcParsVec;
  std::unordered_map<std::string, int> funcParsNamesMap;
  std::vector<FunctionalParameter *> funcParsMap;
  std::unordered_map<int, FuncParFuncType> funcParsFuncMap;
  std::vector<std::vector<int>> funcParsGrid;
  std::vector<std::string> funcParsNamesVec = {};
 
  void CalcNormsBins(std::vector<NormParameter>& norm_parameters, std::vector< std::vector< int > >& xsec_norms_bins);
  M3::float_t CalcWeightSpline(const FarDetectorCoreInfo* MCEvent) const;
  M3::float_t CalcWeightNorm(const FarDetectorCoreInfo* 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 std::string& KinematicParameter);
  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();
 
  //===============================================================================
  //DB Functions required for reweighting functions
  //DB Replace previous implementation with reading bin contents from SampleHandlerFD_array
  void Fill1DHist();
  void Fill2DHist();
 
#ifdef MULTITHREAD
  void FillArray_MP();
#endif
  void FillArray();
 
  inline void ResetHistograms();
  
  //===============================================================================
  //DB Variables required for GetLikelihood
  SampleBinningInfo Binning;
 
  double** SampleHandlerFD_array;
  double** SampleHandlerFD_array_w2;
  double** SampleHandlerFD_data;
  //===============================================================================
 
  //===============================================================================
  //MC variables
  std::vector<FarDetectorCoreInfo> MCSamples;
  std::vector<OscChannelInfo> OscChannels;
  //===============================================================================
 
  //===============================================================================
  //DB Covariance Objects
  //ETA - All experiments will need an xsec, det and osc cov
  //these should be added to SampleHandlerBase to be honest
  ParameterHandlerGeneric *ParHandler = nullptr;
 
  //=============================================================================== 
 
  int nDimensions = M3::_BAD_INT_;
  std::string SampleName;
 
  std::string SampleTitle;
 
  //===========================================================================
  //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< KinematicCut > StoredSelection;
  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 InitialiseSingleFDMCObject();
  void InitialiseSplineObject();
 
  std::vector<std::string> mc_files;
  std::vector<std::string> spline_files;
 
  std::unordered_map<std::string, double> _modeNomWeightMap;
  
  //===============================================================================
  TLegend* THStackLeg = nullptr;
  //===============================================================================
 
  bool FirstTimeW2;
  bool UpdateW2;
 
  TH1D *dathist;
  TH2D *dathist2d;
 
  // binned PDFs
  TH1D* _hPDF1D;
  TH2D* _hPDF2D;
 
  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;
  
};