PredictionSystNumu2017.cxx

Go to the documentation of this file.

#include "3FlavorAna/Prediction/PredictionSystNumu2017.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"

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

#include "OscLib/IOscCalc.h"

#include "TObjString.h"

#include <iostream>

namespace ana
{

  REGISTER_LOADFROM("PredictionSystNumu2017", IPrediction, PredictionSystNumu2017);

  //----------------------------------------------------------------------

  const DummyNumu2017Syst kNumu2017MECq0ShapeSyst     ("numu2017MECq0Shape"           , "MECq0Shape"           , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017MaCCQE_reducedSyst ("numu2017MaCCQE_reduced"       , "MaCCQE_reduced"       , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017RPAShapeRESSyst    ("numu2017RPAShapeRES"          , "RPAShapeRES"          ,         {0, +1, +2}, false); 
  const DummyNumu2017Syst kNumu2017RPAShapeEnhSyst    ("numu2017RPAShapeenh"          , "RPAShapeEnh"          , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017MaCCRESSyst        ("numu2017MaCCRES"              , "MaCCRES"              , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017MaNCRESSyst        ("numu2017MaNCRES"              , "MaNCRES"              , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017MvCCRESSyst        ("numu2017MvCCRES"              , "MvCCRES"              , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017DISvnCC1piSyst     ("numu2017DISvnCC1pi"           , "DISvnCC1pi"           , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017AllSmallSysts      ("numu2017SumSmallXSecJoint2017", "SumSmallXSecJoint2017", {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017NumuSmallSysts     ("numu2017SumSmallXSecNumu2017" , "SumSmallXSecNumu2017" , {-2, -1, 0, +1, +2}, false);

  const DummyNumu2017Syst kNumu2017AllBeamSysts("numu2017beamTransportComb", "beamTransportComb", {-2, -1, 0, +1, +2}, false);

  const DummyNumu2017Syst kNumu2017RelHadEScale("numu2017RelHadEScale2017", "RelHadEScale", {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017AbsHadEScale("numu2017AbsHadEScale2017", "AbsHadEScale", {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017RelMuEScale ("numu2017RelMuEScale2017" , "RelMuEScale" , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017AbsMuEScale ("numu2017AbsMuEScale2017" , "AbsMuEScale" , {-2, -1, 0, +1, +2}, false);

  const DummyNumu2017Syst kNumu2017AbsCalibSyst  ("numu2017Calibration"  , "AbsCalib"     , {-1, 0, +1}, false);
  const DummyNumu2017Syst kNumu2017RelCalibSyst  ("numu2017RelativeCalib", "RelativeCalib", {-2, 0, +2}, false);
  const DummyNumu2017Syst kNumu2017CalibShapeSyst("numu2017CalibShape"   , "CalibShape"   ,     {0, +1}, false);
  const DummyNumu2017Syst kNumu2017LightLevelSyst("numu2017Lightlevel"   , "LightLevel"   , {-1, 0, +1}, false);
  const DummyNumu2017Syst kNumu2017CherenkovSyst ("numu2017Cherenkov"    , "Cherenkov"    ,     {0, +1}, false);

  const DummyNumu2017Syst kNumu2017PPFXPC00Syst ("numu2017ppfx_pc00", "PPFXPC00Sys", {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017PPFXPC01Syst ("numu2017ppfx_pc01", "PPFXPC01"   , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017PPFXPC02Syst ("numu2017ppfx_pc02", "PPFXPC02"   , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017PPFXPC03Syst ("numu2017ppfx_pc03", "PPFXPC03"   , {-2, -1, 0, +1, +2}, false);
  const DummyNumu2017Syst kNumu2017PPFXPC04Syst ("numu2017ppfx_pc04", "PPFXPC04"   , {-2, -1, 0, +1, +2}, false);

  //const DummyNumu2017Syst kNumu2017AbsCalibSyst2Sig("numu2017Calibration", "AbsCalib2Sig", {-2,-1, 0, +1, +2}, false);

  const DummyNumu2017Syst kNumu2017NormSystSig ("numu2017normSig" , "Sig. Norm."   , {0, +1}, false);
  const DummyNumu2017Syst kNumu2017NormSystBkg ("numu2017normBkg" , "Bkg. Norm."   , {0, +1}, false);
  const DummyNumu2017Syst kNumu2017NormSystBoth("numu2017NormBoth", "Sig+Bkg Norm.", {0, +1}, false);

  //----------------------------------------------------------------------
  PredictionSystNumu2017::PredictionSystNumu2017(ENumu2017ExtrapType extrap,
                                                 osc::IOscCalc* osc,
                                                 int WhichQuant,
                                                 const std::vector<const DummyNumu2017Syst*>& systs,
                                                 const std::string fname)
    : PredictionInterp()
  {
    fOscOrigin = osc->Copy();

    std::string extrapStr = "pred_";
    if (extrap == ENumu2017ExtrapType::kNuMu)         extrapStr += "xp_numu";
    if (extrap == ENumu2017ExtrapType::kNumuNoExtrap) extrapStr += "nxp";
    extrapStr += "_";

    std::string VarName = "numu_pred_Energy_Quant"+std::to_string(WhichQuant);

    std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;

    IPrediction *nom;
    std::string NomVar  = extrapStr+"Nominal/noShift/"+VarName;
    fPredNom = ana::LoadFromFile<IPrediction>(fname, NomVar);
    nom      = ana::LoadFromFile<IPrediction>(fname, NomVar).release();
    // --- The light level syst uses an alternative nominal.  Grab that now.
    std::string LNomVar = extrapStr+"Lightlevel/noShift/"+VarName;
    IPrediction* nomLightLevel = ana::LoadFromFile<IPrediction>(fname, LNomVar).release();
  
    for(const DummyNumu2017Syst* syst: systs){
      if(syst == &kNumu2017NormSystBoth ||
         syst == &kNumu2017NormSystSig ||
         syst == &kNumu2017NormSystBkg) continue;

      // Hack.  For sake of distinct names, every syst's fist 8 characters are numu2017
      std::string systName = syst->ShortName();
      systName.erase(0,8);

      std::cout << "Loading " << systName << "..." << std::endl;

      ShiftedPreds sp;
      sp.systName = systName;

      for(int sigma: syst->Shifts()){
        sp.shifts.push_back(sigma);

        if (sigma == 0){
          if (systName == "Lightlevel" || systName == "Cherenkov") {
            std::cout << "\tLooking at systName " << systName << ", so pushing back nomLightLevel " << std::endl;
            sp.preds.push_back(nomLightLevel);
          } else {
            sp.preds.push_back(nom);
          }
          continue;
        }

        std::string sigmaStr = "/";
        if(sigma == -2) sigmaStr += "minusTwo/";
        if(sigma == -1) sigmaStr += "minusOne/";
        if(sigma == +1) sigmaStr += "plusOne/";
        if(sigma == +2) sigmaStr += "plusTwo/"; 


        auto shCore = ana::LoadFromFile<PredictionExtrap>  (fname, extrapStr + systName + sigmaStr +VarName).release();
        sp.preds.push_back(shCore);
      } // end for sigma

      fPreds.emplace(syst, sp);
    } // end for syst

    // Calculate all the ratios etc
    InitFits();

    AddNormSysts(systs);
  }

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

  //----------------------------------------------------------------------
  void PredictionSystNumu2017::AddNormSysts(const std::vector<const DummyNumu2017Syst*>& systs)
  {
    TH1* h = Predict(fOscOrigin).ToTH1(1);
    const int nBins = h->GetNbinsX();
    delete h;


    // Got this from Alex. TODO add reference
    // pythag == sum in quadrature
    // leaving Sig and bkg hooks in place, but they are not used anywhere
    const double normBothSyst = 3.5/100;
    const double normBkgSyst = 3.5/100;
    const double normSigSyst = 3.5/100;

    // TODO find a way to write these three norm systs as a loop
    //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBoth)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "numu2017NormBoth";})){
      ShiftedPreds normsp;
      normsp.systName = "NormBoth";
      normsp.shifts = {0};
      normsp.preds = {fPredNom.get()};
      for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
        normsp.fits[coeff].resize(nBins+2);
        for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
          normsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, normBothSyst, 1));
        }
      }
      normsp.FillRemaps();
      fPreds.emplace(&kNumu2017NormSystBoth, normsp);
    }


    //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystBkg)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "numu2017normBkg";})){
      ShiftedPreds normbksp;
      normbksp.systName = "NumuNormBkg";
      normbksp.shifts = {0};
      normbksp.preds = {fPredNom.get()};
      for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
        normbksp.fits[coeff].resize(nBins+2);
        for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
          normbksp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff != kNumuSurv) ? normBkgSyst : 0, 1));
        }
      }
      normbksp.FillRemaps();
      fPreds.emplace(&kNumu2017NormSystBkg, normbksp);
    }


    //    if(std::count(systs.begin(), systs.end(), &kNumu2017NormSystSig)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "numu2017normSig";})){
      ShiftedPreds normsigsp;
      normsigsp.systName = "NumuNormSig";
      normsigsp.shifts ={0};
      normsigsp.preds = {fPredNom.get()};
      for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
        normsigsp.fits[coeff].resize(nBins+2);
        for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
          normsigsp.fits[coeff][binIdx].push_back(Coeffs(0, 0, (coeff == kNumuSurv) ? normSigSyst : 0, 1));
        }
      }
      normsigsp.FillRemaps();
      fPreds.emplace(&kNumu2017NormSystSig, normsigsp);
    }
  }

  //----------------------------------------------------------------------
  std::unique_ptr<PredictionSystNumu2017> PredictionSystNumu2017::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() == "PredictionSystNumu2017");
    delete tag;

    std::unique_ptr<PredictionSystNumu2017> ret(new PredictionSystNumu2017);
    LoadFromBody(dir, ret.get(), {&kNumu2017NormSystBoth, &kNumu2017NormSystSig, &kNumu2017NormSystBkg});

    // Just enable all norm systs unconditionally. Way too hard to do otherwise
    ret->AddNormSysts({&kNumu2017NormSystBoth, &kNumu2017NormSystSig, &kNumu2017NormSystBkg});

    delete dir;

    return ret;
  }

  //----------------------------------------------------------------------
  void PredictionSystNumu2017::SaveTo(TDirectory* dir, const std::string& name) const
  {
    PredictionInterp::SaveTo(dir, name);

    TDirectory* tmp = gDirectory;

    dir = (TDirectory*)dir->GetDirectory(name.c_str());
    dir->cd();
    // Overwrite "IPrediction"
    TObjString("PredictionSystNumu2017").Write("type");

    dir->Write();
    delete dir;

    tmp->cd();
  }

} // namespace