make_nueFDprediction_kinematics_REW.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 "CAFAna/Extrap/ModularExtrap.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_REW(const std::string& outfilename = "FDprediction_kinematics_REW.root",
                                  const bool hastau = false
                                  )

{
  //********** 0. Basic Setup *********************
    
  // 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 ********
  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
      
  const Var kWOscDumb([](const caf::SRProxy* sr){return sr->mc.nu[0].woscdumb;});
  
  // open the weights file
  TFile *weight_file = new TFile("/nova/app/users/nostom66/nue/extrapsysts/combined/Weights.root","READ");
  
  // construction of trueQ2 weight
  TH1 *h = new TH1();
  gDirectory->GetObject("trueQ2_weight_AllSamples",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 = new TH1();
  gDirectory->GetObject("PtP_weight_AllSamples",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 = new TH1();
  gDirectory->GetObject("CosNumi_weight_AllSamples",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();
  
  
  
  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 int kNumSels = 1;        
  
  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;

  
  // 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, kXSecCVWgt2017*kPPFXFluxCVWgt*kTrueQ2Weight);
  numuDecomp_PtP = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt*kPtPWeight);
  numuDecomp_Cos = new  NumuDecomp (loaders, axisNumu, cutNDNumu, thisShift, kNoShift, kXSecCVWgt2017*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, kXSecCVWgt2017*kPPFXFluxCVWgt );

      // trueQ2 reweighted extrap
      extrapProp_trueQ2[selIdx][varIdx] = NueExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_trueQ2, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
                                                                      
      // PtP reweighted extrap
      extrapProp_PtP[selIdx][varIdx] = NueExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_PtP, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
      // Cos reweighted extrap
      extrapProp_Cos[selIdx][varIdx] = NueExtrap      (loaders,
                                                                      *propDecomp[selIdx][varIdx], *numuDecomp_Cos, 
                                                                      axisNue, axisNumu, 
                                                                      selFD[selIdx], selND[selIdx], cutNDNumu,
                                                                      thisShift, kXSecCVWgt2017*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, 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_trueQ2[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_trueQ2_"+selNames[selIdx]+"_"+varName) ;
      predicProp_PtP[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_PtP_"+selNames[selIdx]+"_"+varName) ;
      predicProp_Cos[selIdx][varIdx]->SaveTo(dFDshi, "nue_pred_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();
  }
}