make_prediction_rhc.C

Go to the documentation of this file.

// WARNING datasets correspond to old files (SA ideal condition RHC)
// Tested in release R17-02-22-legacy-cafs-v1 - UPDATED TO MC PROD3

#ifdef __CINT__
void make_prediction_rhc()
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/NumuSysts.h"
//#include "CAFAna/Systs/NueSystsSecondAna.h"
//#include "CAFAna/Systs/NueSysts.h"
#include "CAFAna/Systs/NueExtrapSysts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "OscLib/IOscCalc.h"
#include "TFile.h"

using namespace ana;


void make_prediction_rhc()
{

  //std::cerr << "\n\n WARNING datasets correspond to old files (SA ideal condition RHC)"
  //<<"\n Tested in release R17-02-22-legacy-cafs-v1 \n \n "; 
  
  TString fileName = "pred_tutorial_nue_numu_rhc-test.root";
  
  ////////////////////////////////////////////////////////////
  // Set up loaders
  ////////////////////////////////////////////////////////////
  Loaders loaders;
  
  //std::string ldrFDNonSwap = "prod_decaf_R16-03-03-prod2reco.e_fd_genie_nonswap_nogenierw_rhc_nova_v08_full_nue_or_numu_contain_v1";//6028 files
  std::string ldrFDNonSwap = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_nonswap_rhc_nova_v08_full_nue_or_numu_or_nus_contain_v1"; //6085 files
  //std::string ldrFDFlxSwap = "prod_decaf_R16-03-03-prod2reco.e_fd_genie_fluxswap_nogenierw_rhc_nova_v08_full_nue_or_numu_contain_v1"; //5994
  std::string ldrFDFlxSwap = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_rhc_nova_v08_full_nue_or_numu_or_nus_contain_v1"; //6029 files
  
  ldrFDNonSwap = "def_snapshot " + ldrFDNonSwap + " with stride 100";//New synonym for 'dataset_def_name_newest_snapshot'
  ldrFDFlxSwap = "defname: "+ldrFDFlxSwap + " with stride 100";//What if returns not snapshots?? // Can't find with stride 100?? //limit 100??
  
  loaders.SetLoaderPath(ldrFDNonSwap,caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(ldrFDFlxSwap,caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  
  loaders.SetSpillCut(kStandardSpillCuts);
  
  ////////////////////////////////////////////////////////////
  // Set up variables, cuts, etc
  ////////////////////////////////////////////////////////////
  
  // Note: Redundancies here for readability and future-proofing. 
  // You can use the kXXX directly in constructors 
  
  const Var weight = kTuftsWeightCC;
  const HistAxis numuAxis = kNumuCCAxis;
  const HistAxis nueAxis =  kNueAxisSecondAnaCVN;
  const Cut numuCutFD = kNumuFD;
  const Cut nueCutFD = kNueSAFDCVNSsb;
  
  //These are silly
  //std::vector< const ISyst * > numuSysts = {&kMECScaleSyst, &kSARelHadESyst, &kBeamAll}; //ORIGINAL-OLD
  std::vector< const ISyst * > numuSysts = {kMECScaleSystSA, &kFDHadEnergyScaleSyst, &kBeamAll}; //Note first one doesn't have ampersand '&'
  //std::vector< const ISyst * > numuSysts(1, &kNumuNormSyst); // Note syntaxis for single-element vector
  
  //std::vector< const ISyst * > nueSysts =  {&kMECScaleSyst, &kBeamAllTransport, &kNA49};// ORIGINAL-OLD
  std::vector< const ISyst * > nueSysts =  {kMECScaleSystSA, &kBeamAllTransport, &kNA49SA}; //Note first one doesn't have ampersand '&'
  
  ////////////////////////////////////////////////////////////
  //Set up the predictions
  ////////////////////////////////////////////////////////////

  NoExtrapGenerator numuGenerator (numuAxis, numuCutFD, weight);
  NoExtrapGenerator nueGenerator  (nueAxis,  nueCutFD,  weight);

  auto numuPrediction = numuGenerator.Generate(loaders);
  auto nuePrediction  = nueGenerator.Generate(loaders);

  auto calc = DefaultOscCalc();

  NoExtrapGenerator numuGeneratorCopy (numuAxis, numuCutFD, weight);
  NoExtrapGenerator nueGeneratorCopy  (nueAxis,  nueCutFD,  weight);
  auto numuPredSysts = new PredictionInterp (numuSysts,calc->Copy(),
                                             numuGeneratorCopy,loaders);
  auto nuePredSysts  = new PredictionInterp (nueSysts,calc->Copy(), 
                                             nueGeneratorCopy,loaders);
  ////////////////////////////////////////////////////////////
  // Go! And save everything
  ////////////////////////////////////////////////////////////
  
  loaders.Go();
  
  TFile * file = new TFile (fileName,"recreate");
  
  numuPrediction->SaveTo(file->mkdir("prediction_numu"));
  nuePrediction->SaveTo(file->mkdir("prediction_nue"));
  numuPredSysts->SaveTo(file->mkdir("prediction_systs_numu"));
  nuePredSysts->SaveTo(file->mkdir("prediction_systs_nue"));

}
#endif