PredictionSyst3Flavor2020.cxx

Go to the documentation of this file.

#include "3FlavorAna/Prediction/PredictionSyst3Flavor2020.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionExtrapSum.h"
#include "3FlavorAna/Cuts/NueCuts2020.h" // kNue2020Binning
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/LoadFromRegistry.h"
#include "3FlavorAna/Systs/DummySystStorage.h" // norm systs
#include "3FlavorAna/Systs/NueAcceptSysts.h"  // accept syst to be execluded from numu
#include "3FlavorAna/Systs/MichelTaggingSyst.h" // michel tagging only for nue FHC
#include "3FlavorAna/Systs/AngleSysts.h" // lepton angle systs
#include "OscLib/IOscCalc.h"

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

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

namespace ana
{

  REGISTER_LOADFROM("PredictionSyst3Flavor2020", IPrediction, PredictionSyst3Flavor2020);

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

    fOscOrigin = osc->Copy();

    fSplitBySign = IsFHC();

    const bool isPtExtrap = extrap == ENu2020ExtrapType::kComboPtExtrap || extrap == ENu2020ExtrapType::kPropPtExtrap;

    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 == ENu2020ExtrapType::kProportional) extrapStr += "xp_prop";
    if (extrap == ENu2020ExtrapType::kCombo)        extrapStr += "xp_combo";
    if (extrap == ENu2020ExtrapType::kNoExtrap)     extrapStr += "nxp";

    if (extrap == ENu2020ExtrapType::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 {
      // --- Set some defaults. These aren't final, so they need to be updated, but that's better than nothing. 
      std::string FileDir = "/cvmfs/nova.osgstorage.org/analysis/3flavor/Predictions/";
      if(NERSC) {
        FileDir = std::string(std::getenv("FCHELPERANA2020_LIB_PATH")) + "/Predictions/";
      }
      if ( isPtExtrap ) {
        std::cout << "\t Loading the Pt Extrap files." << std::endl;
        if        ( IsFHC() && !isFakeData) { fnameSyst = FileDir + "pred_xp_nue_combo_Nue2020Axis_full_FHCAllSysts_PtQuants3_nueconcat_realNDData_20200427.root";
        } else if (!IsFHC() && !isFakeData) { fnameSyst = FileDir + "pred_xp_nue_prop_Nue2020Axis_full_RHCAllSysts_PtQuants3_nueconcat_realNDData_20200427.root";
        } else if ( IsFHC() &&  isFakeData) { fnameSyst = FileDir + "pred_xp_nue_combo_Nue2020Axis_full_FHCAllSysts_PtQuants3_nueconcat_fakeNDData_20200427.root";
        } else if (!IsFHC() &&  isFakeData) { fnameSyst = FileDir + "pred_xp_nue_prop_Nue2020Axis_full_RHCAllSysts_PtQuants3_nueconcat_fakeNDData_20200427.root";
        }
      } else {
        std::cout << "\t Loading the Non Pt Extrap files." << std::endl;
        if        ( IsFHC() && !isFakeData) { fnameSyst = FileDir + "pred_Nue2020Axis_full_FHCAllSysts_nueconcat_20200428.root";
        } else if (!IsFHC() && !isFakeData) { fnameSyst = FileDir + "pred_Nue2020Axis_full_RHCAllSysts_nueconcat_20200428.root";
        } else if ( IsFHC() &&  isFakeData) { fnameSyst = FileDir + "pred_xp_nue_combo_Nue2020Axis_full_FHCAllSysts_nueconcat_fakeNDData_20200427.root";
        } else if (!IsFHC() &&  isFakeData) { fnameSyst = FileDir + "pred_xp_nue_prop_Nue2020Axis_full_RHCAllSysts_nueconcat_fakeNDData_20200427.root";
        }
      }
      std::cout << "Loading the predictions from official file..." << std::endl;
      std::cout << "\t--> " << fnameSyst << std::endl;
    }
    
    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_Nue2020Axis";

    // make PredictionExtendOwning

    IPrediction* nomCore;
    if ( isPtExtrap ) {
      //std::cout << "In the isPtExtrap loading loop." << std::endl;
      nomCore = LoadNuePtExtrapPred( fSyst, extrap, "Nominal", "noShift" );
    } else {
      //std::cout << "In the non-PtExtrap loading loop." << std::endl;
      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, kNue2020Binning, true);
    nom = new PredictExtendOwning(nomCore, nomNoExtrap, kNue2020Binning, true);
    
    for(const ISyst* syst: systs)
    {
      if(syst == &kAna2020NormFHC ||
         syst == &kAna2020NormRHC ||
         syst == &kAna2020NormHornCorr) continue;
         
      // michel tagging only for FHC
      if(!IsFHC() && syst == &kMichelTaggingSyst2020) continue;

      // lepton angle syst only for PtExtrap
      if(!isPtExtrap &&
         (syst == &kLeptonAngleSystNDXZ2020 ||
          syst == &kLeptonAngleSystNDYZ2020 ||
          syst == &kLeptonAngleSystFDXZ2020 ||
          syst == &kLeptonAngleSystFDYZ2020)) 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) {sp.preds.push_back(nom); continue;}

            // Treat Cherenkov And CalibDrift as one sided by
            // fixing negative shifts to the nominal
            if((syst == &kAnaCherenkovSyst || syst == &kAnaCalibDriftSyst) && (sigma == -1 || sigma == -2))
            {
              sp.preds.push_back(nom);
              continue;
            }
        
        const std::string sigmaStr = "/"+SigmaToString(sigma);
        
        IPrediction* shCore;
        if ( isPtExtrap ) {
              //std::cout << "In the isPtExtrap loading loop." << std::endl;
              shCore = LoadNuePtExtrapPred( fSyst, extrap, systName, SigmaToString( sigma ) );
            } else {
              //std::cout << "In the non-PtExtrap loading loop." << std::endl; 
              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, kNue2020Binning, 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, extrap);

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

    delete fSyst;
  }

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

    fOscOrigin = osc->Copy();

    fSplitBySign = IsFHC();

    const bool isPtExtrap = extrap == ENu2020ExtrapType::kNuMuPtExtrap;

    std::cout << "use FakeNDData predictions is "
              << (isFakeData ? "true" : "false") << std::endl;

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

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

    std::string fnameSyst = fname;
    if (fname != "") {
      std::cout << "Loading predictions from unofficial file..." << std::endl;
      std::cout << "\t--> " << fnameSyst << std::endl;
    } else {
      // --- Set some defaults. These aren't final, so they need to be updated, but that's better than nothing. 
      std::string FileDir = "/cvmfs/nova.osgstorage.org/analysis/3flavor/Predictions/";
      if(NERSC) {
        FileDir = std::string(std::getenv("FCHELPERANA2020_LIB_PATH")) + "/Predictions/";
      }
      if ( isPtExtrap ) {
        std::cout << "\t Loading the Pt Extrap files." << std::endl;
        if        (  IsFHC() && !isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_FHCAllSysts_PtQuants3_numuconcat_realNDData_20200427.root";
        } else if ( !IsFHC() && !isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_RHCAllSysts_PtQuants3_numuconcat_realNDData_20200427.root";
        } else if (  IsFHC() &&  isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_FHCAllSysts_PtQuants3_numuconcat_fakeNDData_20200427.root";
        } else if ( !IsFHC() &&  isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_RHCAllSysts_PtQuants3_numuconcat_fakeNDData_20200427.root";
        }
      } else {
        std::cout << "\t Loading the Non-Pt Extrap files." << std::endl;
        if        (  IsFHC() && !isFakeData ) { fnameSyst = FileDir + "pred_NumuEnergy_AllQuantiles_full_FHCAllSysts_numuconcat_20200428.root";
        } else if ( !IsFHC() && !isFakeData ) { fnameSyst = FileDir + "pred_NumuEnergy_AllQuantiles_full_RHCAllSysts_numuconcat_20200428.root";
        } else if (  IsFHC() &&  isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_FHCAllSysts_numuconcat_fakeNDData_20200427.root";
        } else if ( !IsFHC() &&  isFakeData ) { fnameSyst = FileDir + "pred_xp_numu_NumuEnergy_AllQuantiles_full_RHCAllSysts_numuconcat_fakeNDData_20200427.root";
        }
      }
      std::cout << "Loading the predictions from official file..." << std::endl;
      std::cout << "\t--> " << fnameSyst << std::endl;
    }
    
    TFile* fin = TFile::Open(fnameSyst.c_str(), "read");

    IPrediction *nom;
    if ( isPtExtrap )
    {
      fPredNom = std::unique_ptr<IPrediction>( LoadNumuPtExtrapPred( fin, EhadFracQuant, "Nominal", "noShift" ) );
      nom      = LoadNumuPtExtrapPred( fin, EhadFracQuant, "Nominal", "noShift" );
    }
    else
    {
      std::string NomVar  = extrapStr+"Nominal/noShift/"+VarName;
      fPredNom = ana::LoadFrom<IPrediction>(fin, NomVar);
      nom      = ana::LoadFrom<IPrediction>(fin, NomVar).release();
    }

    for(const ISyst* syst: systs){

      if(syst == &kAna2020NormFHC ||
         syst == &kAna2020NormRHC ||
         syst == &kAna2020NormHornCorr ||
         // nue specific systs are not in numu files
         syst == &kNueAcceptSystSignalKinPtExtrap2020FHC ||
         syst == &kNueAcceptSystSignalKinPtExtrap2020RHC ||
         syst == &kNueAcceptSystSignalKin2020FHC ||
         syst == &kNueAcceptSystSignalKin2020RHC ||
         syst == &kMichelTaggingSyst2020) continue;

      // lepton angle syst only for PtExtrap
      if(!isPtExtrap &&
         (syst == &kLeptonAngleSystNDXZ2020 ||
          syst == &kLeptonAngleSystNDYZ2020 ||
          syst == &kLeptonAngleSystFDXZ2020 ||
          syst == &kLeptonAngleSystFDYZ2020)) 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) {sp.preds.push_back(nom);continue;}

        // Treat Cherenkov And CalibDrift as one sided by
        // fixing negative shifts to the nominal
        if((syst == &kAnaCherenkovSyst || syst == &kAnaCalibDriftSyst) &&
           (sigma == -1 || sigma == -2))
            {
              sp.preds.push_back(nom);
              continue;
            }
        
        const std::string sigmaStr = "/"+SigmaToString(sigma)+"/";

        IPrediction* shCore;
        if ( isPtExtrap ) shCore = LoadNumuPtExtrapPred( fin, EhadFracQuant, systName, SigmaToString( sigma ) );
        else              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, extrap);

    if(minimizeMemory){
      MinimizeMemory();
    }

    delete fin;
  }

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

  //----------------------------------------------------------------------
  void PredictionSyst3Flavor2020::
  AddNormSysts(const std::vector<const ISyst*>& systs, ENu2020ExtrapType extrap)
  {
    // Normalisation systematic consists of 6 pieces
    // 4 horn-currents uncorrelated
    // e.g. FHC norm applies shift to fhc,
    // no shift to rhc, and vice versa
    // POT  = 0.55% (docdb-42952)
    // Overlay (docdb-43838):
    // numu = 0.21%/0.48% (FHC/RHC)
    // nue  = 0.41%/0.53% (FHC/RHC)
    // overlay 'numu' applies to numu signal and nue signal
    // overlay 'nue' applies to nue bkg
    // FHC  
    const double kNormNumuFHC  = 0.0059;
    const double kNormNueFHC   = 0.0069;
    // RHC
    const double kNormNumuRHC  = 0.0073;
    const double kNormNueRHC   = 0.0076;
    
    // FD fully-uncorrelated mass uncertainty 
    // mass = 0.19% (docdb-20816)
    // gg mismatch
    // (Genie-Geant mismatch of the primary particle list) 
    // gg_mismatch = 0.8% (docdb-?????)
    // both horn-current correlated
    // Mass + gg_mismatch
    const double kNormHornCorr = 0.0082;
    
    // Remember: none of these normalization systematics will affect the nue rock component
    // final nue pred = rock prediction + fiducial prediction with systematics (this class)
    // That final nue pred is calculated at analysis time, not here 
    // (see GetNuePrediction2020 in joint_fit_2020_loader_tools.h)

    // Include horn-correlated norm 
    AddNormSyst(&kAna2020NormHornCorr, systs, kNormHornCorr, kNormHornCorr);
    
    // Include syst for wrong beam type to keep the fit happy, but set it to
    // zero magnitude
    if(IsFHC())
    {
      // Are we in numu? => keep only one kNormNumu
      if(extrap == ENu2020ExtrapType::kNuMu || extrap == ENu2020ExtrapType::kNuMuNoExtrap || extrap == ENu2020ExtrapType::kNuMuPtExtrap)
      {
        AddNormSyst(&kAna2020NormFHC, systs, kNormNumuFHC, kNormNumuFHC);
        AddNormSyst(&kAna2020NormRHC, systs, .0, .0);
      }
      // Are we in nue? => pass both kNormNumu and kNormNue
      else
      {
        AddNormSyst(&kAna2020NormFHC, systs, kNormNumuFHC, kNormNueFHC);
        AddNormSyst(&kAna2020NormRHC, systs, .0, .0);
      }
    }
    else
    {
      // Are we in numu? => keep only one kNormNumu
      if(extrap == ENu2020ExtrapType::kNuMu || extrap == ENu2020ExtrapType::kNuMuNoExtrap || extrap == ENu2020ExtrapType::kNuMuPtExtrap)
      {
        AddNormSyst(&kAna2020NormFHC, systs, .0, .0);
        AddNormSyst(&kAna2020NormRHC, systs, kNormNumuRHC, kNormNumuRHC);
      }
      // Are we in nue? => pass both kNormNumu and kNormNue
      else
      {
        AddNormSyst(&kAna2020NormFHC, systs, .0, .0);
        AddNormSyst(&kAna2020NormRHC, systs, kNormNumuRHC, kNormNueRHC);
      }
    }
  }
  //----------------------------------------------------------------------
  void PredictionSyst3Flavor2020::
  AddNormSyst(const ISyst* syst,
              const std::vector<const ISyst*>& systs,
              double fracNueApp, double fracOthers)
  {
    // Bail out if we're not actually requested to include this systematic
    if(std::find(systs.begin(), systs.end(), syst) == systs.end()) return;
    
    TH1D* 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);
    
    // The first fits coeff is for nue signal, the rest are for nue bkg and numu
    normsp.fits[CoeffsType::kNueApp].resize(nBins+2);
    if(fSplitBySign) normsp.fitsNubar[CoeffsType::kNueApp].resize(nBins+2);
    for(int binIdx = 0; binIdx < nBins+2; ++binIdx){
      normsp.fits[CoeffsType::kNueApp][binIdx].emplace_back(0, 0, fracNueApp, 1);
      if(fSplitBySign) normsp.fitsNubar[CoeffsType::kNueApp][binIdx].emplace_back(0, 0, fracNueApp, 1);
    }
    
    for(int coeff = 1; 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, fracOthers, 1);
        if(fSplitBySign) normsp.fitsNubar[coeff][binIdx].emplace_back(0, 0, fracOthers, 1);
      }
    }
    
    normsp.FillRemaps();
    normsp.nCoeffs = 1;
    fPreds.emplace(syst, normsp);
  }

  //----------------------------------------------------------------------
  std::unique_ptr<PredictionSyst3Flavor2020> PredictionSyst3Flavor2020::LoadFrom(TDirectory* dir, const std::string& name)
  {
    std::cout << "PredictionSyst3Flavor2020::LoadFrom() not implemented" << std::endl;
    abort();
    
    // TODO: implement with correct Norm syst and PtExtrap loading
    /*
    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);

    assert(tag->GetString() == "PredictionSyst3Flavor2020");

    std::unique_ptr<PredictionSyst3Flavor2020> ret(new PredictionSyst3Flavor2020);
    LoadFromBody(dir, ret.get(), {&kAna2020NormFHC, &kAna2020NormRHC});

    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({&kAna2020NormFHC, &kAna2020NormRHC});

    return ret;
    */
  }

  //----------------------------------------------------------------------
  void PredictionSyst3Flavor2020::SaveTo(TDirectory* dir, const std::string& name) const
  {
    std::cout << "PredictionSyst3Flavor2020::SaveTo() not implemented" << std::endl;
    abort();
    
    // TODO: implement with correct Norm syst (and PtExtrap) saving
    /*
    PredictionInterp::SaveTo(dir);

    TDirectory* tmp = gDirectory;

    dir->cd();
    // Overwrite "IPrediction"
    TObjString("PredictionSyst3Flavor2020").Write("type");

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

    tmp->cd();
    */
  }

  //----------------------------------------------------------------------
  std::string PredictionSyst3Flavor2020::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 << "PredictionSyst3Flavor2020: unknown sigma " << sigma << std::endl;
    abort();
  }

  //----------------------------------------------------------------------
  IPrediction* PredictionSyst3Flavor2020::LoadNuePtExtrapPred( TFile* file, ENu2020ExtrapType extrap, std::string syst_name, std::string sigma_str ) const
  {
    std::string signal_str = "pred_xp_nue_signal";

    std::string bkgd_str = "pred_xp_nue_bkgd_";
    if      ( extrap == ENu2020ExtrapType::kComboPtExtrap ) bkgd_str += "combo";
    else if ( extrap == ENu2020ExtrapType::kPropPtExtrap  ) bkgd_str += "prop";
    else
    {
      std::cout << "ENu2020ExtrapType must be kComboPtExtrap or kPropPtExtrap" << std::endl;
      abort();
    }

    std::string var = "nue_pred_Nue2020Axis";

    std::vector<std::string> pred_paths;

    // Path for background prediction
    pred_paths.push_back( Form( "%s_%s/%s/%s", bkgd_str.c_str(), syst_name.c_str(), sigma_str.c_str(), var.c_str() ) );

    // Paths for signal Pt quantile predictions
    int NPtQuants = 3;
    for ( int ptQuant = 1; ptQuant <= NPtQuants; ++ptQuant )
    {
      std::string signal_path = Form( "%s_%s/%s/%s_Signal_PtQuant%d", signal_str.c_str(), syst_name.c_str(), sigma_str.c_str(), var.c_str(), ptQuant );
      pred_paths.push_back( signal_path );
    }

    return new PredictionExtrapSum( file, pred_paths );
  }

  //----------------------------------------------------------------------
  IPrediction* PredictionSyst3Flavor2020::LoadNumuPtExtrapPred( TFile* file, int EhadFracQuant, std::string syst_name, std::string sigma_str ) const
  {
    if ( EhadFracQuant < 1 || EhadFracQuant > 4 )
    {
      std::cout << "Pt extrap requires 1 <= EhadFracQuant >= 4" << std::endl;
      abort();
    }

    std::string base_path = Form( "pred_xp_numu_%s/%s/numu_pred_NumuEnergy_EhadFracQuant%d", syst_name.c_str(), sigma_str.c_str(), EhadFracQuant ); 

    int NPtQuants = 3;
    std::vector<std::string> pred_paths;
    for ( int ptQuant = 1; ptQuant <= NPtQuants; ++ptQuant )
    {
      pred_paths.push_back( base_path + Form( "_PtQuant%dof%d", ptQuant, NPtQuants ) );
    }

    return new PredictionExtrapSum( file, pred_paths );
  }

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

    std::cout << "PredictionSyst3Flavor2020: Bad beam string: "
              << fBeam << std::endl;
    abort();
  }
} // namespace