PredictionCombinePeriods.cxx

Go to the documentation of this file.

#include "CAFAna/Prediction/PredictionCombinePeriods.h"

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

#include "OscLib/IOscCalc.h"

#include "TDirectory.h"
#include "TH2.h"
#include "TObjString.h"
#include "TVectorD.h"

namespace ana
{
  REGISTER_LOADFROM("PredictionCombinePeriods", IPrediction, PredictionCombinePeriods);

  //----------------------------------------------------------------------
  PredictionCombinePeriods::PredictionCombinePeriods(const std::vector<std::pair<const IPrediction*, double>>& preds)
    : fPreds(preds)
  {
    assert(!fPreds.empty());
  }

  //----------------------------------------------------------------------
  PredictionCombinePeriods::~PredictionCombinePeriods()
  {
    // It's pretty normal for the elements of fPred to have come from the stack
    // for(auto it: fPreds) delete it.first;
  }

  //----------------------------------------------------------------------
  Spectrum PredictionCombinePeriods::Predict(osc::IOscCalc* calc) const
  {
    return PredictComponent(calc,
                            Flavors::kAll,
                            Current::kBoth,
                            Sign::kBoth);
  }

  //----------------------------------------------------------------------
  Spectrum PredictionCombinePeriods::PredictSyst(osc::IOscCalc* calc,
                                                 const SystShifts& syst) const
  {
    return PredictComponentSyst(calc, syst,
                                Flavors::kAll,
                                Current::kBoth,
                                Sign::kBoth);
  }

  //----------------------------------------------------------------------
  Spectrum PredictionCombinePeriods::
  PredictComponent(osc::IOscCalc* calc,
                   Flavors::Flavors_t flav,
                   Current::Current_t curr,
                   Sign::Sign_t sign) const
  {
    return PredictComponentSyst(calc, SystShifts::Nominal(),
                                flav, curr, sign);
  }

  //----------------------------------------------------------------------
  Spectrum PredictionCombinePeriods::
  PredictComponentSyst(osc::IOscCalc* calc,
                       const SystShifts& syst,
                       Flavors::Flavors_t flav,
                       Current::Current_t curr,
                       Sign::Sign_t sign) const
  {
    Eigen::ArrayXd asum;

    std::vector<std::string> labels;
    std::vector<Binning> bins;

    double potsum = 0;
    double livesum = 0;
    for(auto it: fPreds){
      const IPrediction* pred = it.first;
      const double targetPOT = it.second;
      const Spectrum s = pred->PredictComponentSyst(calc, syst,
                                                    flav, curr, sign);

      Eigen::ArrayXd h = s.GetEigen(targetPOT);
      potsum += targetPOT;
      // No. FD MC is going to have zero (unknown) livetime anyway. If someone
      // needs this livetime accounting, they can come and think hard about it
      // and then reinstate it.
      // livesum += s.Livetime()*targetPOT/s.POT();

      if(asum.size() > 0){
        asum += h;
      }
      else{
        asum = h;
        labels = s.GetLabels();
        bins = s.GetBinnings();
      }
    }

    return Spectrum(std::move(asum),
                    HistAxis(labels, bins),
                    potsum, livesum);
  }

  //----------------------------------------------------------------------
  OscillatableSpectrum PredictionCombinePeriods::ComponentCC(int from, int to) const
  {
    // TODO this is a total copy-paste job from the above.

    Eigen::MatrixXd asum;

    std::vector<std::string> labels;
    std::vector<Binning> bins;

    double potsum = 0;
    double livesum = 0;
    for(auto it: fPreds){
      const IPrediction* pred = it.first;
      const double targetPOT = it.second;
      const OscillatableSpectrum s = pred->ComponentCC(from, to);

      Eigen::MatrixXd h = s.GetEigen(targetPOT);
      potsum += targetPOT;
      // No. FD MC is going to have zero (unknown) livetime anyway. If someone
      // needs this livetime accounting, they can come and think hard about it
      // and then reinstate it.
      // livesum += s.Livetime()*targetPOT/s.POT();

      if(asum.size() > 0){
        asum += h;
      }
      else{
        asum = h;
        labels = s.GetLabels();
        bins = s.GetBinnings();
      }
    }

    return OscillatableSpectrum(std::move(asum),
                                HistAxis(labels, bins),
                                potsum, livesum);
  }

  //----------------------------------------------------------------------
  //nc
  Spectrum PredictionCombinePeriods::ComponentNCTotal() const
  {
    osc::NoOscillations noosc;
    return PredictComponent(&noosc, Flavors::kAll, Current::kNC, Sign::kBoth);
  }
  Spectrum PredictionCombinePeriods::ComponentNC() const
  {
    osc::NoOscillations noosc;
    return PredictComponent(&noosc, Flavors::kAll, Current::kNC, Sign::kNu);
  }
  Spectrum PredictionCombinePeriods::ComponentNCAnti() const
  {
    osc::NoOscillations noosc;
    return PredictComponent(&noosc, Flavors::kAll, Current::kNC, Sign::kAntiNu);
  }

  //end nc

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

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

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

    int idx = 0;
    for(auto it: fPreds)
      it.first->SaveTo(dir, TString::Format("pred%d", idx++).Data());

    TVectorD pots(fPreds.size());
    for(unsigned int i = 0; i < fPreds.size(); ++i)
      pots[i] = fPreds[i].second;

    pots.Write("pots");

    dir->Write();
    delete dir;

    tmp->cd();
  }

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

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "PredictionCombinePeriods");
    delete tag;

    TVectorD* pots = (TVectorD*)dir->Get("pots");
    assert(pots);

    std::vector<std::pair<const IPrediction*, double>> preds;

    for(int i = 0; i < pots->GetNrows(); ++i){
      const IPrediction* pred = ana::LoadFrom<IPrediction>(dir, TString::Format("pred%d", i).Data()).release();
      assert(pred);
      preds.emplace_back(pred, (*pots)[i]);
    }

    delete dir;

    return std::unique_ptr<PredictionCombinePeriods>(new PredictionCombinePeriods(preds));
  }
}