make_nueFDprediction_kinematics.C

Go to the documentation of this file.

// Macro to make nominal and shifted predictions
// Takes names and datasets as input
// Create files for nominal and shifted samples
// set systematicSample to false for beam or genie systematics
// USAGE
// cafe -bq make_nueFDprediction_shifted.C bool "ShortName" "Long Name"        
//          "filename" "nd data" "nd mc dataset" "FDMC nonswap" "swap" "tau"
// See full examples in run_nueFDprediction*.sh 
//hastau,plus2,plus1,minu1,minu2,knoshift,compareto

// Analysis
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/SALoaders.h"
#include "CAFAna/Analysis/Style.h"

// Core
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"

// Cuts
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"


// Decomp
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/BENDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"

// Extrap
#include "3FlavorAna/Extrap/NueExtrap.h"

// Predict
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"

// Systs
#include "CAFAna/Systs/Systs.h"

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

#include "OscLib/OscCalcDumb.h"

#include "StandardRecord/Proxy/SRProxy.h" 


#include <iostream>
#include <iomanip>

using namespace ana;

void CheckFileOverwrite(TString);

osc::OscCalcDumb calc;

void make_nueFDprediction_kinematics( const std::string& outfilename = "FDprediction_kinematics.root",
                                  const bool hastau = false
)

{
  // ND loaders
  std::string loaderNDData = "prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
  std::string loaderNDMC = "prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
  
  // FD loaders
  std::string loaderFDnonswapMC = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
  std::string loaderFDfluxswapMC = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
  std::string loaderFDtauswapMC = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";

  //******* 1. Set up selection and variables ********
  //Only needed here. Keep updated! 
  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
  
  const Var kWOscDumb([](const caf::SRProxy* sr){return sr->mc.nu[0].woscdumb;});
    
  const int kNumSels = 1;
  const int kNumVars = 5;
  
  const HistDef defs[kNumVars] = {
      {"recoE", {"E_{reco} (GeV)", Binning::Simple(10, 0, 5), kNueEnergy2017}},
      {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
      {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
      {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
      {"CosTheta", {"CosTheta", Binning::Simple(500,0,1), kCosTheta}}
  };
  
  const Cut selFD[kNumSels] = {
      kNue2017FDAllSamples
      };
  const Cut selND[kNumSels] = {
      kNue2017NDPresel && kNueSecondAnaCVNeSsb
      };
  const std::string selNames[kNumSels] = {
      "AllSamples"
      };
  const Cut cutNDNumu = { 
      kNumuNDCvn
  };

  CheckFileOverwrite(outfilename.c_str());

  //********** 2. Set up loaders ********************* 
  SANueExtrapLoaders loaders;
  
  loaders.SetLoaderPath(loaderNDData, caf::kNEARDET, ana::Loaders::DataMC::kData);
  loaders.SetLoaderPath(loaderNDMC, caf::kNEARDET, ana::Loaders::DataMC::kMC);
  
  loaders.SetLoaderPath(loaderFDnonswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam, ana::Loaders::SwappingConfig::kNonSwap);
  loaders.SetLoaderPath(loaderFDfluxswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam,ana::Loaders::SwappingConfig::kFluxSwap);
  if(hastau)
    loaders.SetLoaderPath(loaderFDtauswapMC, caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam, ana::Loaders::SwappingConfig::kTauSwap);
  else
    loaders.DisableLoader(caf::kFARDET, ana::Loaders::DataMC::kMC, DataSource::kBeam,ana::Loaders::SwappingConfig::kTauSwap);

  loaders.SetSpillCut(kStandardSpillCuts);

  //******* 3. Set up the extrapolations and go *******
  SystShifts thisShift = kNoShift;

  NumuDecomp         * numuDecomp;
  ProportionalDecomp * propDecomp[kNumSels][kNumVars];
  ModularExtrap      * extrapProp[kNumSels][kNumVars];
  PredictionExtrap   * predicProp[kNumSels][kNumVars];
  PredictionNoExtrap * predicNoXP[kNumSels][kNumVars];  
  
  numuDecomp = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
  
  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
    for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){  
      
      const HistAxis& axisNue = defs[varIdx].axis;
      
      propDecomp[selIdx][varIdx] = new ProportionalDecomp  (loaders,
                                                                      axisNue, selND[selIdx],
                                                                      thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt );
      
      extrapProp[selIdx][varIdx] = new NueExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
      
      predicProp[selIdx][varIdx] = new PredictionExtrap    (&extrapProp[selIdx][varIdx]);
      
      predicNoXP[selIdx][varIdx] = new PredictionNoExtrap  (loaders,
                                                                      axisNue,
                                                                      selFD[selIdx], thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
      
    }//vars
  }//sels
  loaders.Go();
  
  //********** 4. Save everything ********************     
  TFile * file = new TFile(outfilename.c_str(),"RECREATE");
  TDirectory * dFDshi = file->mkdir("prediction");
  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
    for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
      TString varName = defs[varIdx].name.c_str();
      predicProp[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_"+selNames[selIdx]+"_"+varName) ;
      predicNoXP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextrap_"+selNames[selIdx]+"_"+varName) ;
    }  
  }     
  file->Close();

} // end

//****************************************
void CheckFileOverwrite(TString outfilename){
  if(!gSystem->AccessPathName(outfilename.Data())){
    std::cout << "\n\nThis will overwrite " << outfilename 
              << "\n\nAre you sure you want to continue? y/n ";
    string input = ""; getline(cin,input);
    if (input!="y") abort();
  }
}