SingleSampleExperiment.cxx

Go to the documentation of this file.

#include "CAFAna/Experiment/SingleSampleExperiment.h"

#include "CAFAna/Core/HistCache.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/StanUtils.h"
#include "CAFAna/Core/Utilities.h"

#include "OscLib/func/IOscCalculator.h"

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

namespace ana
{
  const CosmicBkgScaleSyst kCosmicBkgScaleSyst;

  //----------------------------------------------------------------------
  SingleSampleExperiment::SingleSampleExperiment(const IPrediction* pred,
                                                 const Spectrum& data,
                                                 const Spectrum& cosmic,
                                                 double cosmicScaleError,
                                                 bool PoissonError)
    : fMC(pred), fData(data),
      fCosmic(cosmic.ToTH1(data.Livetime(), kLivetime)),
      fCosmicScaleError(cosmicScaleError),
      fPoissonError(PoissonError)
  {
    // Fill the vector of errors in case of PoissonError = true
    if(fPoissonError){
      TH1D* cosmicErrHist = cosmic.ToTH1(cosmic.Livetime(), kLivetime);

       TFeldmanCousins fc(0.6827);
       for(int i = 1; i<= cosmicErrHist->GetNbinsX(); i++){
           double y = cosmicErrHist->GetBinContent(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}); 
       }

       HistCache::Delete(cosmicErrHist);
    }

  }

  //----------------------------------------------------------------------
  SingleSampleExperiment::~SingleSampleExperiment()
  {
    HistCache::Delete(fCosmic);
  }

  //----------------------------------------------------------------------
  TH1D* SingleSampleExperiment::
  PredHistIncCosmics(osc::IOscCalculator* calc,
                     const SystShifts& syst) const
  {
    SystShifts systNoCosmic = syst;
    systNoCosmic.SetShift(&kCosmicBkgScaleSyst, 0);

    const Spectrum pred = fMC->PredictSyst(calc, systNoCosmic);
    TH1D* hpred = pred.ToTH1(fData.POT());
    if(fCosmic){
      if(fCosmicScaleError != 0 && !fPoissonError){
        const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
        hpred->Add(fCosmic, scale);
      }
      if(fPoissonError){
        TH1D* hcosm = HistCache::Copy(fCosmic);
        double sigma = syst.GetShift(&kCosmicBkgScaleSyst);
        for(int i = 1; i<= hcosm->GetNbinsX(); i++){
            if (sigma>0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].first));
            if (sigma<0) hcosm->SetBinContent(i, hcosm->GetBinContent(i)*(1+sigma*fBinErrors[i-1].second));
        }
        hpred->Add(hcosm);
        HistCache::Delete(hcosm);
      }
      else{
        hpred->Add(fCosmic);
      }
    }
    return hpred;
  }

  //----------------------------------------------------------------------
  double SingleSampleExperiment::ChiSq(osc::IOscCalculatorAdjustable* calc,
                                       const SystShifts& syst) const
  {
    TH1D* hpred = PredHistIncCosmics(calc, syst);
    TH1D* hdata = fData.ToTH1(fData.POT());

    // full namespace qualification to avoid degeneracy with method inherited from IExperiment
    const double ll = ana::LogLikelihood(hpred, hdata);

    HistCache::Delete(hpred);
    HistCache::Delete(hdata);

    return ll;
  }

  //----------------------------------------------------------------------
  void SingleSampleExperiment::
  Derivative(osc::IOscCalculator* calc,
             const SystShifts& shift,
             std::unordered_map<const ISyst*, double>& dchi) const
  {

    const double pot = fData.POT();

    std::unordered_map<const ISyst*, std::vector<double>> dp;
    for(auto it: dchi) dp[it.first] = {};
    fMC->Derivative(calc, shift, pot, dp);

    if(dp.empty()){ // prediction doesn't implement derivatives
      dchi.clear(); // pass on that info to our caller
      return;
    }

    TH1D* hpred = PredHistIncCosmics(calc, shift);
    TH1D* hdata = fData.ToTH1(pot);

    for(auto& it: dchi){
      if(it.first != &kCosmicBkgScaleSyst){
        it.second += LogLikelihoodDerivative(hpred, hdata, dp[it.first]);
      }
      else{
        const unsigned int N = fCosmic->GetNbinsX()+2;
        const double* ca = fCosmic->GetArray();
        std::vector<double> cosErr(N);
        if(fPoissonError){
          double sigma = shift.GetShift(&kCosmicBkgScaleSyst);
          for(unsigned int i = 1; i<= fBinErrors.size(); i++){
              double err = 0;
              if(sigma<0) err = fBinErrors[i-1].second;
              if(sigma>0) err = fBinErrors[i-1].first;
              cosErr[i] = ca[i]*err;
          }
        }
        else {
           for(unsigned int i = 1; i <= N; ++i) cosErr[i] = ca[i]*fCosmicScaleError;
        }
        it.second += LogLikelihoodDerivative(hpred, hdata, cosErr);
      }
    }
    exit(0);
    HistCache::Delete(hpred);
    HistCache::Delete(hdata);
  }

  //----------------------------------------------------------------------
  stan::math::var SingleSampleExperiment::LogLikelihood(osc::_IOscCalculatorAdjustable<stan::math::var> *osc,
                                                        const SystShifts &syst) const
  {
    auto  pred = syst.IsNominal()
                 ? fMC->Predict(osc).ToBins(fData.POT())
                 : fMC->PredictSyst(osc, syst).ToBins(fData.POT());
    TH1D* hdata = fData.ToTH1(fData.POT());

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

    HistCache::Delete(hdata);

    return ll;
  }
  //----------------------------------------------------------------------
  void SingleSampleExperiment::SaveTo(TDirectory* dir) const
  {
    TDirectory* tmp = dir;

    dir->cd();
    TObjString("SingleSampleExperiment").Write("type");

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

    if(fCosmic) fCosmic->Write("cosmic");

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<SingleSampleExperiment> SingleSampleExperiment::LoadFrom(TDirectory* dir)
  {
    TObjString* ptag = (TObjString*)dir->Get("type");
    assert(ptag);
    assert(ptag->GetString() == "SingleSampleExperiment");

    assert(dir->GetDirectory("mc"));
    assert(dir->GetDirectory("data"));
    

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

    TH1D* cosmic = 0;
    if(dir->Get("cosmic")) cosmic = (TH1D*)dir->Get("cosmic");

    auto ret = std::make_unique<SingleSampleExperiment>(mc, *data);
    if(cosmic) ret->fCosmic = cosmic;
    return ret;
  }
}