PredictionSystJointDemo.cxx

Go to the documentation of this file.

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

#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Cuts/NueCuts2018.h" // kNue2018Binning
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/LoadFromRegistry.h"

#include "3FlavorAna/Systs/DummySystStorage.h" // norm systs

#include "OscLib/IOscCalc.h"

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

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

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

  //----------------------------------------------------------------------
  // Nue constructor
  PredictionSystJointDemo::PredictionSystJointDemo(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 = "Nue_";
    if (extrap == kProportional) extrapStr += "PropExtrap";
    if (extrap == kCombo)        extrapStr += "ComboExtrap";
    if (extrap == kNoExtrap)     extrapStr += "NoExtrap";

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

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

    IPrediction *nom;
    std::string fnameSyst;
  
    fnameSyst = fname;
    TFile fSyst(fnameSyst.c_str());
    if(fSyst.IsZombie()){
      std::cout << "Couldn't open " << fnameSyst << std::endl;
      abort();
    }

    // make PredictionExtendOwning

    IPrediction* nomCore = ana::LoadFrom<IPrediction>(&fSyst, extrapStr + "/Nominal/NoShift").release();
    PredictionNoExtrap* nomNoExtrap = ana::LoadFrom<PredictionNoExtrap>(&fSyst, "Nue_NoExtrap/Nominal/NoShift").release();

    fPredNom = std::make_unique<PredictExtendOwning>(nomCore, nomNoExtrap, kNue2018Binning, true);
    nom = new PredictExtendOwning(nomCore, nomNoExtrap, kNue2018Binning, 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;
                // 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 + "/LightNom/NoShift").release();
                PredictionNoExtrap* nomLLNoExtrap = ana::LoadFrom<PredictionNoExtrap>(
                    &fSyst, "Nue_NoExtrap/LightNom/NoShift").release();

                nomLightLevel = new PredictExtendOwning(
                    nomLLCore, nomLLNoExtrap, true);
    
                sp.preds.push_back(nomLightLevel);
              }
            else
              sp.preds.push_back(nom);
            continue;
          }
        
        const std::string sigmaStr = "/"+SigmaToString(sigma);
        
        IPrediction* shCore = ana::LoadFrom<PredictionExtrap>  (&fSyst, extrapStr + "/"+systName + sigmaStr).release();
        PredictionNoExtrap* shNoExtrap = ana::LoadFrom<PredictionNoExtrap>(&fSyst, "Nue_NoExtrap/" + systName + sigmaStr).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();
    }
  }

  //----------------------------------------------------------------------
  // Numu constructor
  PredictionSystJointDemo::PredictionSystJointDemo(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 = "Numu_";
    if (extrap == kNuMu)             extrapStr += "Extrap";
    if (extrap == kNuMuNoExtrap)     extrapStr += "NoExtrap";
    extrapStr += "_";

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

    TFile fin(fname.c_str());

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

            sp.preds.push_back(nomLightLevel);
          } else {
            sp.preds.push_back(nom);
          }
          continue;
        }
        
        // const std::string sigmaStr = "/"+SigmaToString(sigma)+"/";
        const std::string sigmaStr = "/"+SigmaToString(sigma);

        // auto shCore = ana::LoadFrom<PredictionExtrap>  (&fin, extrapStr + systName + sigmaStr +VarName).release();
        auto shCore = ana::LoadFrom<PredictionExtrap>  (&fin, extrapStr + VarName+"/"+systName + sigmaStr).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();
    }
  }

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

  //----------------------------------------------------------------------
  void PredictionSystJointDemo::
  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 PredictionSystJointDemo::
  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<PredictionSystJointDemo> PredictionSystJointDemo::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() == "PredictionSystJointDemo");
    delete tag;

    std::unique_ptr<PredictionSystJointDemo> ret(new PredictionSystJointDemo);
    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 PredictionSystJointDemo::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("PredictionSystJointDemo").Write("type");

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

    dir->Write();
    delete dir;

    tmp->cd();
  }

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

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

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


} // namespace