make_nueFDprediction_kinematics_RHC_REW.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/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"
#include "3FlavorAna/Extrap/NueRHCExtrap.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 "CAFAna/Vars/HistAxes.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_RHC_REW(const std::string sample, const std::string wfile, const std::string& outfilename = "FDprediction_kinematics_REW.root")

{
  //********** 0. Basic Setup *********************
    
  // ND loaders
  std::string loaderNDData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns";
  std::string loaderNDMC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
  
  // FD loaders
  std::string loaderFDnonswapMC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1";
  std::string loaderFDfluxswapMC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_fluxswap_rhc_nova_v08_full_v1";
  std::string loaderFDtauswapMC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_tau_rhc_nova_v08_full_v1";
  
  //******* 1. Set up selection and variables ********
  const HistAxis axisNumu = kNumuCCAxis;
  
  const int kNumSels = 1;
  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}}
  };
  
  Cut nueFDsel = kNue2018FDAllSamples;
  
  /*if (sample == "CORE") nueFDsel = kNue2018FD;
  else if (sample == "PERI") nueFDsel = kNue2018FDPeripheral;*/
  
  const Cut selFD[kNumSels] = {
      nueFDsel
      };
  const Cut selND[kNumSels] = {
      kNue2018NDPresel && kNue2018CVNCut
      };
  const std::string selNames[kNumSels] = {
      sample
      };
      
  const Cut cutNDNumu = { 
      kNumuCutND2018
  };

  // open the weights file
    TFile *weight_file = new TFile(wfile.c_str(),"READ");
  
    // construction of trueQ2 weight
    TH1 *h;
    gDirectory->GetObject(("trueQ2_weight_"+sample).c_str(),h);
  
    std::cout << h->GetNbinsX() << endl;
    int nbins = h->GetNbinsX();
    double bins[nbins], binslow[nbins], binshigh[nbins];
    for(int i = 1; i <= nbins; ++i){
        bins[i-1] = h->GetBinContent(i);
        binslow[i-1] = h->GetBinLowEdge(i);
        binshigh[i-1] = h->GetBinLowEdge(i) + h->GetBinWidth(i);
    }
    
    const Var kTrueQ2Weight([nbins, &binslow, &binshigh, &bins](const caf::SRProxy *sr)
      {
        for(int i = 0; i < nbins; ++i){
            if((kRecoQ2(sr) > binslow[i]) && (kRecoQ2(sr) < binshigh[i])){return bins[i];}
        }
        return 1.0;
      }
    );
  
    //construction of PtP weight
    TH1 *hptp;
    gDirectory->GetObject(("PtP_weight_"+sample).c_str(),hptp);
  
    std::cout << hptp->GetNbinsX() << endl;
    int nbins_ptp = hptp->GetNbinsX();
    double bins_ptp[nbins_ptp], binslow_ptp[nbins_ptp], binshigh_ptp[nbins_ptp]; 
    for(int i = 1; i <= nbins_ptp; ++i){
        bins_ptp[i-1] = hptp->GetBinContent(i);
        binslow_ptp[i-1] = hptp->GetBinLowEdge(i);
        binshigh_ptp[i-1] = hptp->GetBinLowEdge(i) + hptp->GetBinWidth(i);
    }
  
    const Var kPtPWeight([nbins_ptp, &binslow_ptp, &binshigh_ptp, &bins_ptp](const caf::SRProxy *sr)
      {
        for(int i = 0; i < nbins_ptp; ++i){
            if((kPtP(sr) > binslow_ptp[i]) && (kPtP(sr) < binshigh_ptp[i])){return bins_ptp[i];}
        }
        return 1.0;
      }
    );
  
    //construction of Cos weight
    TH1 *hcos;
    gDirectory->GetObject(("CosNumi_weight_"+sample).c_str(),hcos);
  
    std::cout << hcos->GetNbinsX() << endl;
    int nbins_cos = hptp->GetNbinsX();
    double bins_cos[nbins_cos], binslow_cos[nbins_cos], binshigh_cos[nbins_cos]; 
    for(int i = 1; i <= nbins_cos; ++i){
        bins_cos[i-1] = hcos->GetBinContent(i);
        binslow_cos[i-1] = hcos->GetBinLowEdge(i);
        binshigh_cos[i-1] = hcos->GetBinLowEdge(i) + hcos->GetBinWidth(i);
    }
  
    const Var kCosWeight([nbins_cos, &binslow_cos, &binshigh_cos, &bins_cos](const caf::SRProxy *sr)
      {
        for(int i = 0; i < nbins_cos; ++i){
            if((kCosNumi(sr) > binslow_cos[i]) && (kCosNumi(sr) < binshigh_cos[i])){return bins_cos[i];}
        }
        return 1.0;
      }
    );
  
    // close weights file
    weight_file->Close();
    
  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;
  
  // reweighted numudecomps
  NumuDecomp         * numuDecomp_trueQ2;
  NumuDecomp         * numuDecomp_PtP;
  NumuDecomp         * numuDecomp_Cos;
  
  // propdecomp
  ProportionalDecomp * propDecomp[kNumSels][kNumVars];
  
  // extraps from different numudecomps
  ModularExtrap      * extrapProp_trueQ2[kNumSels][kNumVars];
  ModularExtrap      * extrapProp_PtP[kNumSels][kNumVars];
  ModularExtrap      * extrapProp_Cos[kNumSels][kNumVars];
   
  PredictionExtrap   * predicProp_trueQ2[kNumSels][kNumVars];
  PredictionExtrap   * predicProp_PtP[kNumSels][kNumVars];
  PredictionExtrap   * predicProp_Cos[kNumSels][kNumVars];
  
  PredictionNoExtrap * predicNoXP[kNumSels][kNumVars];
  
  numuDecomp_trueQ2 = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kTrueQ2Weight);
  numuDecomp_PtP = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kPtPWeight);
  numuDecomp_Cos = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt*kCosWeight);
  
  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 );

      // trueQ2 reweighted extrap
      extrapProp_trueQ2[selIdx][varIdx] = NueRHCExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_trueQ2, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
                                                                      
      // PtP reweighted extrap
      extrapProp_PtP[selIdx][varIdx] = new NueRHCExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_PtP, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
      // Cos reweighted extrap
      extrapProp_Cos[selIdx][varIdx] = new NueRHCExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_Cos, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
      
      
      predicProp_trueQ2[selIdx][varIdx] = new PredictionExtrap    (extrapProp_trueQ2[selIdx][varIdx]);
      predicProp_PtP[selIdx][varIdx] = new PredictionExtrap    (extrapProp_PtP[selIdx][varIdx]);
      predicProp_Cos[selIdx][varIdx] = new PredictionExtrap    (extrapProp_Cos[selIdx][varIdx]);
      
      predicNoXP[selIdx][varIdx] = new PredictionNoExtrap  (loaders,
                                                                      axisNue,
                                                                      selFD[selIdx], thisShift, kXSecCVWgt2018*kPPFXFluxCVWgt);
      
    }//vars
  }//sels

  loaders.Go();
  
  PredictionExtendToPeripheral * predicExte_trueQ2[kNumSels][kNumVars];
  PredictionExtendToPeripheral * predicExte_PtP[kNumSels][kNumVars];
  PredictionExtendToPeripheral * predicExte_Cos[kNumSels][kNumVars];
  
  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){    
    for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
      predicExte_trueQ2[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_trueQ2[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
      predicExte_PtP[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_PtP[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
      predicExte_Cos[selIdx][varIdx] = new PredictionExtendToPeripheral (predicProp_Cos[selIdx][varIdx], predicNoXP[selIdx][varIdx], kNue2018Binning, false);
    }
  }
  
  //********** 4. Save everything ********************     
  TFile * file = new TFile(outfilename.c_str(),"RECREATE");
  for(unsigned int selIdx = 0; selIdx < kNumSels; ++selIdx){
    TDirectory* dFDshi = file->mkdir((selNames[selIdx]+"_"+sample).c_str());
    for(unsigned int varIdx = 0; varIdx < kNumVars; ++varIdx){
      TString varName = defs[varIdx].name.c_str();
      predicExte_trueQ2[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_trueQ2_"+selNames[selIdx]+"_"+varName) ;
      predicExte_PtP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_PtP_"+selNames[selIdx]+"_"+varName) ;
      predicExte_Cos[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_Cos_"+selNames[selIdx]+"_"+varName) ;
      predicProp_trueQ2[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_trueQ2_"+selNames[selIdx]+"_"+varName) ;
      predicProp_PtP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_PtP_"+selNames[selIdx]+"_"+varName) ;
      predicProp_Cos[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_noextend_Cos_"+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();
  }
}