makePred.C

Go to the documentation of this file.

// Makes all the predictions I need: 
//  SA-like, A17-like (per quartile) and OR

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"

#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/Analysis/Calcs.h"
#include "OscLib/OscCalcPMNSOpt.h"

#include "CAFAna/Prediction/PredictionExtrap.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "CAFAna/Prediction/PredictionInterp.h"


#include "TFile.h"


using namespace ana;

void setLoadersPathNumu2017(Loaders& loaders){

  std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu2017";
  std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017";
  std::string fd_nonswap  = "prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string fd_fluxswap = "prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
  std::string fd_tau      = "prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";

  //////////////////// set loaders ////////////////////
  loaders.SetLoaderPath( nd_data, caf::kNEARDET,        Loaders::kData, ana::kBeam, 
                         Loaders::kNonSwap );
  loaders.SetLoaderPath( nd_nonswap,  caf::kNEARDET,    Loaders::kMC,   ana::kBeam,
                         Loaders::kNonSwap );
  loaders.SetLoaderPath( fd_nonswap,  caf::kFARDET,     Loaders::kMC,   ana::kBeam,
                         Loaders::kNonSwap );
  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET,     Loaders::kMC,   ana::kBeam,
                         Loaders::kFluxSwap);
  loaders.SetLoaderPath( fd_tau,      caf::kFARDET,     Loaders::kMC,   ana::kBeam,
                         Loaders::kTauSwap );

}

void makePred(){                        

  // custom osc. parameters (non-max mixing) Get the numu result parameters!
  osc::OscCalcPMNSOpt calc_nonMax;
  ana::ResetOscCalcToDefault(&calc_nonMax);
  calc_nonMax.SetDmsq32(2.43071e-3);
  calc_nonMax.SetTh23(asin(sqrt(0.471974)));
  calc_nonMax.SetdCP(M_PI*0.915869);  // preliminary Ana2017 

  // Get FD files for HadEFrac vs. neutrino energy distribution
  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);

  Loaders loaders;
  loaders.SetSpillCut(kStandardSpillCuts);
  setLoadersPathNumu2017(loaders);

   const Cut kNumuCosmicRejSA(
         [](const caf::SRProxy* sr)
         { return ( sr->sel.cosrej.anglekal  > 0.5 &&
              sr->sel.cosrej.numuSAcontpid > 0.535 &&
              //sr->sel.cosrej.numucontpid > 0.535 &&
              sr->slc.nhit < 400);
         }
         );

  const Cut kNumuFDSA = kNumuQuality && kNumuContainFD && kNumuNCRej && kNumuCosmicRejSA;

  // make a vector with an entry for each quantiles
  std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
  std::vector<PredictionInterp*> predictionVec;

  NumuExtrapGenerator NumuCCExtrap_SA (kNumuCCAxis, kNumuFDSA, kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
  NumuExtrapGenerator NumuCCExtrap_OR (kNumuCCOptimisedAxis, kNumuCutFD2017 || kNumuFDSA, kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);

  PredictionInterp prediction_SA ({}, &calc_nonMax, NumuCCExtrap_SA, loaders);
  PredictionInterp prediction_OR ({}, &calc_nonMax, NumuCCExtrap_OR, loaders);

  for(auto thisCut : QuantCuts){
    NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuCutFD2017, thisCut && kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017));
    predictionVec.push_back(new PredictionInterp({}, &calc_nonMax, *(NumuCCExtrapVec.back()), loaders));
  } // loop over quantiles

  ////////////////////
  // GOOOOOOOO!!!!!
  loaders.Go();
  ////////////////////

  // Now we should save the predictions
  std::cerr << "Saving Predictions" << std::endl;
  std::string OutFileName = "pred_3A_prod3_AllQuants.root";

  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
  prediction_SA.SaveTo(f, "prediction_SA");
  prediction_OR.SaveTo(f, "prediction_OR"); 
  f->Close();                               
  delete f;


  for(size_t quant=1; quant <= predictionVec.size();quant++){
    auto prediction = predictionVec[quant-1];
    OutFileName = "pred_3A_prod3_" + (std::string)Form("Quant%d", quant) + ".root";
    std::cerr << "Saving Quantile " << quant << " to Prediction file: " << OutFileName << std::endl;
    TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
    prediction -> SaveTo(f, "prediction");
    f->Close();
    delete f;
  } 

}