make_nueFDprediction_kinematics_FHC.C

Go to the documentation of this file.

// 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 "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"

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

// Extrap
#include "CAFAna/Extrap/ModularExtrap.h"

// Predict
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"

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

// Vars
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/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_FHC( const std::string& outfilename = "FDprediction_kinematics_FHC.root")

{
  // ND loaders
  std::string loaderNDData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
  std::string loaderNDMC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  
  // FD loaders
  std::string loaderFDnonswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string loaderFDfluxswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1";
  std::string loaderFDtauswapMC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1";

  //******* 1. Set up selection and variables ********
  //Only needed here. Keep updated! 
  const HistAxis axisNumu = kNumuCCAxis;
    
  const int kNumSels = 3;
  const int kNumVars = 5;
  
  const HistDef defs[kNumVars] = {
      {"cvn2d", kNue2018Axis},
      {"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] = {
      kNue2018FD,
      kNue2018FDPeripheral,
      kNue2018FDAllSamples
      };
  const Cut selND[kNumSels] = {
      kNue2018NDPresel && kNue2018CVNCut,
      kNue2018NDPresel && kNue2018CVNCut,
      kNue2018NDPresel && kNue2018CVNCut
      };
  const std::string selNames[kNumSels] = {
      "CORE",
      "PERI",
      "ALL"
      };
  const Cut cutNDNumu = { 
      kNumuCutND2018
  };

  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);
  loaders.SetLoaderPath(loaderFDtauswapMC, 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];  
  PredictionExtendToPeripheral * predicExte[kNumSels][kNumVars];
  
  numuDecomp = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*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, kXSecCVWgt2018*kPPFXFluxCVWgt );
      
      extrapProp[selIdx][varIdx] = new NueExtrap(NueExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt));
      
      predicProp[selIdx][varIdx] = new PredictionExtrap    (extrapProp[selIdx][varIdx]);
      
      predicNoXP[selIdx][varIdx] = new PredictionNoExtrap  (loaders,
                                                                      axisNue,
                                                                      selFD[selIdx], thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
      
    }//vars
  }//sels
  loaders.Go();
  
  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
    for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
      predicExte[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
    }
  }
  
  //********** 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();
      predicExte[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_"+selNames[selIdx]+"_"+varName) ;
      predicProp[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_"+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();
  }
}