PredictionSystJoint2018.cxx

Go to the documentation of this file.

#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/LoadFromRegistry.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Systs/DummySystStorage.h" // norm systs

#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Cuts/NueCuts2018.h" // kNue2018Binning

#include "OscLib/IOscCalc.h"

#include "TObjString.h"
#include "TROOT.h"

#include <iostream>
#include <unistd.h>

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

  //----------------------------------------------------------------------
  // Nue constructor
  PredictionSystJoint2018::PredictionSystJoint2018(ENu2018ExtrapType extrap,
                                                   osc::IOscCalc* osc,
                                                   const std::vector<const ISyst*>& systs,
                                                   std::string beam,
                                                   bool isFakeData,
                                                   bool mergePeripheral,
                                                   const std::string fname,
                                                   bool minimizeMemory)
    : fBeam(beam)
  {
    gROOT->SetMustClean(false);

    fOscOrigin = osc->Copy();

    fSplitBySign = IsFHC();

    std::cout << "mergePeripheral is "
              << (mergePeripheral ? "true" : "false") << std::endl;

    std::cout << "use FakeNDData predictions is "
              << (isFakeData ? "true" : "false") << std::endl;
    std::string extrapStr = "pred_";
    if (extrap == kProportional) extrapStr += "xp_prop";
    if (extrap == kCombo)        extrapStr += "xp_combo";
    if (extrap == kNoExtrap)     extrapStr += "nxp";

    if (extrap == kFake)         extrapStr += "xp_prop";

    extrapStr += "_";

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

    IPrediction *nom;

    std::string fnameSyst;
  
    if (fname != ""){
      std::cout << "Loading predictions from unofficial file..." << std::endl;
      std::cout << "\t--> " << fname << std::endl;
      fnameSyst = fname;
    }

    else if (IsFHC() && !isFakeData)
      fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if (!IsFHC() && !isFakeData)
      fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if (IsFHC() && isFakeData)
      fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnFakeNDData_2018-05-23/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_22_2018_pass1_merged.root";
    else if (!IsFHC() && isFakeData)
      fnameSyst = "/nova/ana/nu_e_ana/Ana2018/Predictions/provisional_singles/NueProd4SystematicsOnFakeNDData_2018-05-23/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_22_2018_pass1_merged.root";
    
    TFile* fSyst = TFile::Open(fnameSyst.c_str(), "read");
    if(fSyst->IsZombie()){
      std::cout << "Couldn't open " << fnameSyst << std::endl;
      abort();
    }

    std::string varname ="/nue_pred_Nue2018Axis";

    // make PredictionExtendOwning

    IPrediction* nomCore = ana::LoadFrom<IPrediction>(fSyst, extrapStr + "Nominal/noShift" + varname).release();
    PredictionNoExtrap* nomNoExtrap = ana::LoadFrom<PredictionNoExtrap>(fSyst, "pred_nxp_Nominal/noShift" + varname).release();

    fPredNom = std::make_unique<PredictExtendOwning>(nomCore, nomNoExtrap, kNue2018Binning, true);
    nom = new PredictExtendOwning(nomCore, nomNoExtrap, kNue2018Binning, true);
   

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

    IPrediction* nomLightLevel;
    
    IPrediction* nomLLCore = ana::LoadFrom<IPrediction>(
        fSyst, extrapStr + "Lightlevel/noShift"+varname).release();
    PredictionNoExtrap* nomLLNoExtrap = ana::LoadFrom<PredictionNoExtrap>(
        fSyst, "pred_nxp_Lightlevel/noShift"+varname).release();

    nomLightLevel = new PredictExtendOwning(
        nomLLCore, nomLLNoExtrap, true);
    
    for(const ISyst* syst: systs){

      if(syst == &kAna2018NormFHC ||
         syst == &kAna2018NormRHC) continue;
      
      std::string systName = syst->ShortName();
      
      std::cout << "Loading " << systName << "..." << std::endl;
      
      ShiftedPreds sp;
      sp.systName = systName;
      
      for(int sigma: GetDummyShifts(syst))
        {
        sp.shifts.push_back(sigma);
        
        if (sigma == 0)
          {
            if(syst == &kAnaLightlevelSyst ||
               syst == &kAnaCherenkovSyst)
              {
                std::cout << "\tLooking at systName " << systName  << ", so pushing back nomLightLevel " << std::endl;
                sp.preds.push_back(nomLightLevel);
              }
            else
              sp.preds.push_back(nom);
            continue;
          }
        
        const std::string sigmaStr = "/"+SigmaToString(sigma);
        
        IPrediction* shCore = ana::LoadFrom<IPrediction>(fSyst, extrapStr + systName + sigmaStr +varname).release();
        PredictionNoExtrap* shNoExtrap = ana::LoadFrom<PredictionNoExtrap>(fSyst, "pred_nxp_" + systName + sigmaStr +varname).release();

        sp.preds.push_back(new PredictExtendOwning(shCore, shNoExtrap, kNue2018Binning, true));
       
      } // end for sigma
     // std::cout<<" test find "<<systName<<" "<<syst->find(systName)<<std::endl;
      fPreds.emplace(syst, sp);
    } // end for syst

    // Calculate all the ratios etc
    InitFits();

    AddNormSysts(systs);

    // Clean up. Don't use if you need to make debug plots
    if(minimizeMemory){
      MinimizeMemory();
    }

    delete fSyst;
  }

  //----------------------------------------------------------------------
  // Numu constructor
  PredictionSystJoint2018::PredictionSystJoint2018(ENu2018ExtrapType extrap,
                                                   osc::IOscCalc* osc,
                                                   const std::string beam,
                                                   int WhichQuant,
                                                   const std::vector<const ISyst*>& systs,
                                                   const std::string fname,
                                                   bool minimizeMemory)
    : fBeam(beam)
  {
    gROOT->SetMustClean(false);

    fOscOrigin = osc->Copy();

    fSplitBySign = IsFHC();

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

    std::string VarName = "numu_pred_Energy_Quant"+std::to_string(WhichQuant);
    if (WhichQuant == -1) VarName = "numu_pred_NumuEnergy";
    
    std::cout << "\nLooking at VarName " << VarName << " --> Loading nominals now..." << std::endl;

    TFile* fin = TFile::Open(fname.c_str(), "read");

    IPrediction *nom;
    std::string NomVar  = extrapStr+"Nominal/noShift/"+VarName;
    fPredNom = ana::LoadFrom<IPrediction>(fin, NomVar);
    nom      = ana::LoadFrom<IPrediction>(fin, NomVar).release();
    // --- The light level syst uses an alternative nominal.  Grab that now.
    std::string LNomVar = extrapStr+"Lightlevel/noShift/"+VarName;
    IPrediction* nomLightLevel = ana::LoadFrom<IPrediction>(fin, LNomVar).release();

    for(const ISyst* syst: systs){

      if(syst == &kAna2018NormFHC ||
         syst == &kAna2018NormRHC ||
         syst->ShortName() == "RPAShapeRES2018Test" // || // As we aren't using it in Ana18, and so not in the Pred. files.
         // -- When using Ana18 Predictions, need to comment out the following XSec Systs.
         //syst->ShortName() == "RPAShapeenh2019" ||
         //syst->ShortName() == "RPAShapesupp2019" ||
         //syst->ShortName() == "RPAShapeRES2019"  ||
         //syst->ShortName() == "MECShape2018RPAFixNu" ||
         //syst->ShortName() == "MECShape2018RPAFixAntiNu"
         ) continue;

      std::string systName = syst->ShortName();

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

      ShiftedPreds sp;
      sp.systName = systName;

      for(int sigma: GetDummyShifts(syst)){
        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;
        }
        
        const std::string sigmaStr = "/"+SigmaToString(sigma)+"/";

        auto shCore = ana::LoadFrom<PredictionExtrap>  (fin, 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);

    if(minimizeMemory){
      MinimizeMemory();
    }

    delete fin;
  }

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

  //----------------------------------------------------------------------
  void PredictionSystJoint2018::
  AddNormSysts(const std::vector<const ISyst*>& systs)
  {
    // DocDb-27230-v7 and DocDb-27821 >v1
    // leave horn currents uncorrelated
    // e.g. FHC norm applies shift to fhc,
    // no shift to rhc, and vice versa
    const double kNormFHC = 0.0144;
    const double kNormRHC = 0.0064;

    // Include syst for wrong beam type to keep the fit happy, but set it to
    // zero magnitude
    if(IsFHC()){
      AddNormSyst(&kAna2018NormFHC, systs, kNormFHC);
      AddNormSyst(&kAna2018NormRHC, systs, 0);
    }
    else{
      AddNormSyst(&kAna2018NormFHC, systs, 0);
      AddNormSyst(&kAna2018NormRHC, systs, kNormRHC);
    }
  }

  //----------------------------------------------------------------------
  void PredictionSystJoint2018::
  AddNormSyst(const ISyst* syst,
              const std::vector<const ISyst*>& systs,
              double frac)
  {
    // Bail out if we're not actually requested to include this systematic
    if(std::find(systs.begin(), systs.end(), syst) == systs.end()) return;

    TH1* h = Predict(fOscOrigin).ToTH1(1);
    const int nBins = h->GetNbinsX();
    delete h;

    ShiftedPreds normsp;
    normsp.systName = syst->ShortName();
    normsp.shifts = {0};
    normsp.preds = {fPredNom.get()};
    normsp.fits.resize(kNCoeffTypes);
    if(fSplitBySign) normsp.fitsNubar.resize(kNCoeffTypes);

    for(int coeff = 0; coeff < kNCoeffTypes; ++coeff){
      normsp.fits[coeff].resize(nBins+2);
      if(fSplitBySign) normsp.fitsNubar[coeff].resize(nBins+2);
      for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
        normsp.fits[coeff][binIdx].emplace_back(0, 0, frac, 1);
        if(fSplitBySign) normsp.fitsNubar[coeff][binIdx].emplace_back(0, 0, frac, 1);
      }
    }

    normsp.FillRemaps();

    normsp.nCoeffs = 1;

    fPreds.emplace(syst, normsp);
  }

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

    std::unique_ptr<PredictionSystJoint2018> ret(new PredictionSystJoint2018);
    LoadFromBody(dir, ret.get(), {&kAna2018NormFHC, &kAna2018NormRHC});

    const TObjString* beam = (TObjString*)dir->Get("beam");
    assert(beam);
    ret->fBeam = beam->GetString();

    ret->fSplitBySign = ret->IsFHC();

    // Just enable all norm systs unconditionally. Way too hard to do otherwise
    // NB fBeam was set above. That's necessary for correctness.
    ret->AddNormSysts({&kAna2018NormFHC, &kAna2018NormRHC});

    delete dir;

    return ret;
  }

  //----------------------------------------------------------------------
  void PredictionSystJoint2018::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("PredictionSystJoint2018").Write("type");

    TObjString(fBeam.c_str()).Write("beam");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::string PredictionSystJoint2018::SigmaToString(int sigma) const
  {
    if(sigma == -2) return "minusTwo";
    if(sigma == -1) return "minusOne";
    if(sigma ==  0) return "";
    if(sigma == +1) return "plusOne";
    if(sigma == +2) return "plusTwo";

    std::cout << "PredictionSystJoint2018: unknown sigma " << sigma << std::endl;
    abort();
  }

  //----------------------------------------------------------------------
  bool PredictionSystJoint2018::IsFHC() const
  {
    if(fBeam == "fhc" ) return true;
    if(fBeam == "rhc") return false;

    std::cout << "PredictionSystJoint2018: Bad beam string: "
              << fBeam << std::endl;
    abort();
  }


} // namespace