PredictionSystNue2017.cxx

Go to the documentation of this file.

#include "3FlavorAna/Prediction/PredictionSystNue2017.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("PredictionSystNue2017", IPrediction, PredictionSystNue2017);

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

  const DummyNue2017Syst kNue2017MECq0ShapeSyst      ("nue2017MECq0Shape",
                                                      "MECq0Shape",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017MaCCQE_reducedSyst  ("nue2017MaCCQE_reduced",
                                                      "MaCCQE_reduced",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017RPAShapeRESSyst     ("nue2017RPAShapeRES",
                                                      "RPAShapeRES",
                                                      {0, +1, +2}, false);
  const DummyNue2017Syst kNue2017RPAShapeEnhSyst     ("nue2017RPAShapeenh",
                                                      "RPAShapeEnh",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017MaCCRESSyst         ("nue2017MaCCRES",
                                                      "MaCCRES",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017MaNCRESSyst         ("nue2017MaNCRES",
                                                      "MaNCRES",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017MvCCRESSyst         ("nue2017MvCCRES",
                                                      "MvCCRES",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017RadCorrNueSyst      ("nue2017radcorrnue",
                                                      "RadCorrNue",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017RadCorrNueBarSyst   ("nue2017radcorrnuebar",
                                                      "RadCorrNueBar",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017SecondClassCurrSyst ("nue20172ndclasscurr",
                                                      "SecondClassCurr",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017DISvnCC1piSyst      ("nue2017DISvnCC1pi",
                                                      "DISvnCC1pi",
                                                      {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017AllSmallXSecSyst    ("nue2017SumSmallXSecJoint2017", 
                                                      "SumSmallXSecJoint2017", 
                                                      {-2, -1, 0, +1, +2}, false);


  const DummyNue2017Syst kNue2017AbsCalibSyst   ("nue2017Calibration",
                                                 "AbsCalib",
                                                 {-1, 0, +1}, false);
  const DummyNue2017Syst kNue2017RelCalibSyst   ("nue2017RelativeCalib",
                                                 "RelCalib",
                                                 {-2, 0, +2}, false);
  const DummyNue2017Syst kNue2017CalibShapeSyst ("nue2017CalibShape",
                                                 "CalibShape",
                                                 {0, +1}, false);
  const DummyNue2017Syst kNue2017LightLevelSyst ("nue2017Lightlevel",
                                                 "LightLevel",
                                                 {-1, 0, +1}, false);
  const DummyNue2017Syst kNue2017CherenkovSyst  ("nue2017Cherenkov",
                                                 "Cherenkov",
                                                 {0, +1}, false);

  const DummyNue2017Syst kNue2017PPFXPC00Syst ("nue2017ppfx_pc00",
                                               "PPFXPC00Sys",
                                               {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017PPFXPC01Syst ("nue2017ppfx_pc01",
                                               "PPFXPC01",
                                               {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017PPFXPC02Syst ("nue2017ppfx_pc02",
                                               "PPFXPC02",
                                               {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017PPFXPC03Syst ("nue2017ppfx_pc03",
                                               "PPFXPC03",
                                               {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017PPFXPC04Syst ("nue2017ppfx_pc04",
                                               "PPFXPC04",
                                               {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017BeamTransportSyst ("nue2017beamTransportComb",
                                                    "Combined Beam Transport",
                                                    {-2, -1, 0, +1, +2}, false);

  const DummyNue2017Syst kNue2017ExtrapSigSyst ("nue2017extrap_signalkin",
                                                "Nue Signal Extrapolation",
                                                {-2, -1, 0, +1, +2}, false);
  const DummyNue2017Syst kNue2017ExtrapBkgSyst ("nue2017extrap_bkg",
                                                "Nue Bkg Extrapolation",
                                                {-2, -1, 0, +1, +2}, false);

  const DummyNue2017Syst kNue2017NormSystSig ("nue2017NueNormSig",
                                              "Sig. Norm.",
                                              {0, +1}, false);
  const DummyNue2017Syst kNue2017NormSystBkg ("nue2017NueNormBkg",
                                              "Bkg. Norm.",
                                              {0, +1}, false);
  const DummyNue2017Syst kNue2017NormSystBoth("nue2017NormBoth",
                                              "Sig+Bkg Norm.",
                                              {0, +1}, false);


  //----------------------------------------------------------------------
  PredictionSystNue2017::PredictionSystNue2017(ENue2017ExtrapType extrap,
                                               osc::IOscCalc* osc,
                                               const std::vector<const DummyNue2017Syst*>& systs,
                                               const std::string fname)
    : PredictionInterp()
  {
    fOscOrigin = osc->Copy();

    std::string extrapStr = "pred_";
    if (extrap == ENue2017ExtrapType::kProportional) extrapStr += "xp_prop";
    if (extrap == ENue2017ExtrapType::kCombo)        extrapStr += "xp_combo";
    if (extrap == ENue2017ExtrapType::kNoExtrap)     extrapStr += "nxp";
    extrapStr += "_";

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

    IPrediction *nom;

    //std::string varname ="";
    std::string varname ="/nue_pred_EnergyCVN2D";
    if (extrap == ENue2017ExtrapType::kNoExtrap){
      fPredNom = LoadFromFile<IPrediction>
        (fname, extrapStr + "Nominal/noShift" + varname);
      nom = LoadFromFile<IPrediction>
        (fname, extrapStr + "Nominal/noShift" + varname).release();
    }
    else{
      // We need to pull out a PredictionExtrap and a PredictionNoExtrap
      // to make the final PredictionExtendToPeripheral
      PredictionExtrap*   nomCore     = LoadFromFile<PredictionExtrap>
        (fname, extrapStr + "Nominal/noShift" + varname).release();
      PredictionNoExtrap* nomNoExtrap = LoadFromFile<PredictionNoExtrap>
        (fname, "pred_nxp_Nominal/noShift" + varname).release();
      fPredNom = std::unique_ptr<PredictionExtendToPeripheral>
        (new PredictionExtendToPeripheral(nomCore,nomNoExtrap,true));
      
      nom = new PredictionExtendToPeripheral(nomCore,nomNoExtrap,true);

    }

    // We have exactly 1 syst for the 2017 that is non-standard is some way
    // The light level syst uses an alternative nominal.  Grab that now.
    IPrediction* nomLightLevel;

    if (extrap == ENue2017ExtrapType::kNoExtrap){
      nomLightLevel = LoadFromFile<IPrediction>(fname, extrapStr + "Lightlevel/noShift"+
                                                varname).release();
    }
    else{
      PredictionExtrap*   nomLLCore     = LoadFromFile<PredictionExtrap>
        (fname, extrapStr + "Lightlevel/noShift"+varname).release();
      PredictionNoExtrap* nomLLNoExtrap = LoadFromFile<PredictionNoExtrap>
        (fname,"pred_nxp_Nominal/noShift"+varname).release();
      nomLightLevel = new PredictionExtendToPeripheral(nomLLCore,nomLLNoExtrap,true);
    }

    for(const DummyNue2017Syst* syst: systs){
      if(syst == &kNue2017NormSystBoth ||
         syst == &kNue2017NormSystSig ||
         syst == &kNue2017NormSystBkg) continue;

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

      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";

        //std::string varname ="";
        std::string varname ="/nue_pred_EnergyCVN2D"; 
        auto shCore = LoadFromFile<PredictionExtrap>  
          (fname, extrapStr + systName + sigmaStr +varname).release();
        auto shNoExtrap = LoadFromFile<PredictionNoExtrap>
          (fname, "pred_nxp_" + systName + sigmaStr +varname).release();
        sp.preds.push_back(new PredictionExtendToPeripheral(shCore,shNoExtrap,true));
      } // end for sigma

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

    // Calculate all the ratios etc
    InitFits();

    AddNormSysts(systs);
  }

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

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

    // Got this 3.5% from Alex. TODO add reference 
    // leaving Sig and bkg hooks in place, but they are not used in the fit
    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(), &kNue2017NormSystBoth)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "nue2017NormBoth";})){
      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(&kNue2017NormSystBoth, normsp);
    }


    //    if(std::count(systs.begin(), systs.end(), &kNue2017NormSystBkg)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "nue2017NueNormBkg";})){
      ShiftedPreds normbksp;
      normbksp.systName = "NueNormBkg";
      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 != kNueApp) ? normBkgSyst : 0, 1));
        }
      }
      normbksp.FillRemaps();
      fPreds.emplace(&kNue2017NormSystBkg, normbksp);
    }


    //    if(std::count(systs.begin(), systs.end(), &kNue2017NormSystSig)){
    if(std::count_if(systs.begin(), systs.end(),
                     [](const ISyst* s){return s->ShortName() == "nue2017NueNormSig";})){
      ShiftedPreds normsigsp;
      normsigsp.systName = "NueNormSig";
      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 == kNueApp) ? normSigSyst : 0, 1));
        }
      }
      normsigsp.FillRemaps();
      fPreds.emplace(&kNue2017NormSystSig, normsigsp);
    }
  }

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

    std::unique_ptr<PredictionSystNue2017> ret(new PredictionSystNue2017);
    LoadFromBody(dir, ret.get(), {&kNue2017NormSystBoth, &kNue2017NormSystSig, &kNue2017NormSystBkg});

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

    delete dir;

    return ret;
  }

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

    TDirectory* tmp = gDirectory;

    dir = dir->GetDirectory(name.c_str());
    dir->cd();
    // Overwrite "IPrediction"
    TObjString("PredictionSystNue2017").Write((name+"/type").c_str());

    dir->Write();
    delete dir;

    tmp->cd();
  }

} // namespace