SingleSampleExperiment.cxx

Go to the documentation of this file.

#include "CAFAna/Experiment/SingleSampleExperiment.h"

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/LoadFromRegistry.h"
#include "CAFAna/Core/StanUtils.h"
#include "CAFAna/Core/EigenUtils.h"

#include "OscLib/IOscCalc.h"

#include "Utilities/func/Stan.h"

#include "TDirectory.h"
#include "TObjString.h"
#include "TH1.h"
#include "TFeldmanCousins.h"

namespace ana
{
  REGISTER_LOADFROM("SingleSampleExperiment", IExperiment, SingleSampleExperiment);

  const CosmicBkgScaleSyst kCosmicBkgScaleSyst;

  //----------------------------------------------------------------------
  SingleSampleExperiment::SingleSampleExperiment(const IPrediction* pred,
                                                 const Spectrum& data,
                                                 const Spectrum& cosmic,
                                                 double cosmicScaleError,
                                                 bool PoissonError)
    : fMC(pred), fData(data),
      fCosmic(cosmic),
      fCosmicScaleError(cosmicScaleError),
      fPoissonError(PoissonError)
  {
    // Fill the vector of errors in case of PoissonError = true
    if(fPoissonError){
      // Note this is scaled to the cosmic's intrinsic stats level
      Eigen::ArrayXd cosmicErrHist = cosmic.GetEigen(cosmic.Livetime(), kLivetime);

      TFeldmanCousins fc(0.6827);
      for(int i = 1; i<= cosmicErrHist.size()-2; i++){
        double y = cosmicErrHist[i];
        double errup = 0;
        double errlow = 0;
        if(y!=0){
          if ( y < 50 )   errlow  = (y - fc.CalculateLowerLimit(y,0))/y;
          else            errlow = (sqrt(y))/y;
          if ( y < 30 )   errup  = (fc.CalculateUpperLimit(y,0)-y)/y;
          else            errup = (sqrt(y))/y;
        }
        fBinErrors.push_back({errup, errlow}); 
      }
    }

  }

  //----------------------------------------------------------------------
  SingleSampleExperiment::~SingleSampleExperiment()
  {
  }

  //----------------------------------------------------------------------
  template <typename T>
  Eigen::Array<T, Eigen::Dynamic, 1>
  SingleSampleExperiment::PredHistIncCosmics(osc::_IOscCalc<T>* calc,
                                             const SystShifts& syst) const
  {
    using EigenArray = Eigen::Array<T, Eigen::Dynamic, 1>;

    SystShifts systNoCosmic = syst;
    systNoCosmic.RemoveShift(&kCosmicBkgScaleSyst);

    Spectrum specPred = fMC->PredictSyst(calc, systNoCosmic);
    EigenArray apred;
    if (specPred.HasStan())
    {
      if constexpr (std::is_same_v<T, stan::math::var>)
        apred = specPred.GetEigenStan(fData.POT());
    }
    else
      apred = specPred.GetEigen(fData.POT());

    if(fCosmic.has_value()){
      // note the conversion here from double to stan::math::var if T is stan::math::var
      EigenArray acosm = fCosmic->GetEigen(fData.Livetime(), kLivetime);
      if(fCosmicScaleError != 0 && !fPoissonError){
        const T scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
        acosm *= scale;
      }
      if(fPoissonError){
        T sigma = syst.GetShift(&kCosmicBkgScaleSyst);
        for(int i = 1; i<= acosm.size()-2; i++){
          if (sigma>0) acosm[i] *= (1+sigma*fBinErrors[i-1].first);
          if (sigma<0) acosm[i] *= (1+sigma*fBinErrors[i-1].second);
        }
      }

      apred += acosm;
    }

    return apred;
  }

  template Eigen::ArrayXd
  SingleSampleExperiment::PredHistIncCosmics(osc::IOscCalc* calc,
                                                     const SystShifts& syst) const;
  template Eigen::ArrayXstan
  SingleSampleExperiment::PredHistIncCosmics(osc::IOscCalcStan* calc,
                                                              const SystShifts& syst) const;

  //----------------------------------------------------------------------
  double SingleSampleExperiment::ChiSq(osc::IOscCalcAdjustable* calc,
                                       const SystShifts& syst) const
  {
    Eigen::ArrayXd pred = PredHistIncCosmics(calc, syst);
    Eigen::ArrayXd data = fData.GetEigen(fData.POT());

    // full namespace qualification to avoid degeneracy with method inherited from IExperiment
    return ana::LogLikelihood(pred, data);
  }

  //----------------------------------------------------------------------
  stan::math::var SingleSampleExperiment::LogLikelihood(osc::IOscCalcAdjustableStan* osc,
                                                        const SystShifts &syst) const
  {
    Eigen::ArrayXstan pred = PredHistIncCosmics(osc, syst);
    Eigen::ArrayXd data = fData.GetEigen(fData.POT());

    // fully-qualified so that we get the one in StanUtils.h
    return ana::LogLikelihood(pred, data) / -2.;  // LogLikelihood(), confusingly, returns chi2=-2*LL
  }

  //----------------------------------------------------------------------
  void SingleSampleExperiment::SaveTo(TDirectory* dir, const std::string& name) const
  {
    TDirectory* tmp = dir;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("SingleSampleExperiment").Write("type");

    fMC->SaveTo(dir, "mc");
    fData.SaveTo(dir, "data");

    if(fCosmic.has_value()) fCosmic->SaveTo(dir, "cosmic");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<SingleSampleExperiment> SingleSampleExperiment::LoadFrom(TDirectory* dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    TObjString* ptag = (TObjString*)dir->Get("type");
    assert(ptag);
    assert(ptag->GetString() == "SingleSampleExperiment");
    delete ptag;

    const IPrediction* mc = ana::LoadFrom<IPrediction>(dir, "mc").release();
    const std::unique_ptr<Spectrum> data = Spectrum::LoadFrom(dir, "data");

    auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);

    if(dir->Get("cosmic")){
      ret->fCosmic = *ana::LoadFrom<Spectrum>(dir, "cosmic");
    }

    delete dir;

    return ret;
  }
}