mrbrem_get_initial_spectra.C

Go to the documentation of this file.

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"

//#include "CAFAna/Cuts/BPFCuts.h"
//#include "CAFAna/Vars/BPFVars.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Vars/MRVars.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <cassert>

#include "TFile.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TH1.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TLeaf.h"
#include "stdio.h"
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <string>


using namespace ana;
using namespace std;

void mrbrem_get_initial_spectra(std::string mode = "fhc")
{

  std::string filenameCosmicsData = "prod_mrbremcaf_R19-11-18-prod5reco.i_fd_cosmic_fhc_full_v1_goodruns_underP9";
  std::string filenameCRYMC = "prod_mrbremcaf_R19-11-18-prod5reco.m_fd_cry_fhc_full_v3";
  std::string filenameGENIEFluxswap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1";
 
  if (mode == "rhc"){
      filenameCosmicsData = "prod_mrbremcaf_R19-11-18-prod5reco.i_fd_cosmic_rhc_full_v1_goodruns";
      filenameCRYMC = "prod_mrbremcaf_R19-11-18-prod5reco.m_fd_cry_rhc_full_v3";
      filenameGENIEFluxswap = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1";
  }

  SpectrumLoader loaderCosmicsData(filenameCosmicsData, kCosmic);
  SpectrumLoader loaderCRYMC(filenameCRYMC, kCosmic);
  SpectrumLoader loaderGENIEFluxswap(filenameGENIEFluxswap);
  loaderGENIEFluxswap.SetSpillCut(kStandardSpillCuts);
  loaderCosmicsData.SetSpillCut(kStandardSpillCuts_FD_Prod4MotivatedDQ);
  
  struct HistDef
  {
     std::string name;
     HistAxis axis;        
  };

  const int kNumVars = 1;
  const int kNumSels = 9;

  const HistDef defs[kNumVars] = {
        {"shwangle", {"Cos(#theta_{Z})", Binning::Simple(60,0,1.0), kShwAngle}},
  };

  struct CutDef
  {
     std::string name;
     Cut cut;
  };

  const CutDef sels[kNumSels] = {
      {"nocuts", kNoCut},
      {"core_presel", kNue2020CorePresel_MRBrem},
      {"periph_presel", kNue2020PeripheralPresel_MRBrem},
      {"periph_preselbdt", kNue2020PeripheralPresel_MRBrem && kCosPIDPeriBDT > 0.57},
      {"periph_preselcvn", kNue2020PeripheralPresel_MRBrem && kCVNe_looseptp >= 0.97},
      {"periph_preselen", kNue2020PeripheralPresel_MRBrem && kShwE > 1},
      {"core_full", kNue2020FD_MRBrem},
      {"periph_full", kNue2020FDPeripheral_MRBrem},
      {"periph_fullen", kNue2020FDPeripheral_MRBrem && kShwE > 1},
  };

  Spectrum* spec_mrbrem_data[kNumSels][kNumVars];
  Spectrum* spec_mrbrem_cry[kNumSels][kNumVars];
  Spectrum* spec_genie_fd[kNumSels][kNumVars];

  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
      for (int varIdx=0; varIdx<kNumVars;varIdx++){
           const HistAxis& axis = defs[varIdx].axis;
           spec_mrbrem_data[selIdx][varIdx] = new Spectrum(loaderCosmicsData, axis, sels[selIdx].cut && kInCosmicTimingWindow_FD_MR);
           spec_mrbrem_cry[selIdx][varIdx] = new Spectrum(loaderCRYMC, axis, sels[selIdx].cut && kInCosmicTimingWindow_FD_MR);
           spec_genie_fd[selIdx][varIdx] = new Spectrum(loaderGENIEFluxswap, axis, sels[selIdx].cut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
      }
  }

  loaderCosmicsData.Go();
  loaderCRYMC.Go();
  loaderGENIEFluxswap.Go();

  TFile *file = new TFile(("MRbrem_extract_weights_"+mode+".root").c_str(), "recreate");
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
      TDirectory* d = file->mkdir(sels[selIdx].name.c_str());
      for (int varIdx=0; varIdx<kNumVars;varIdx++){
           const char* name = defs[varIdx].name.c_str();
           spec_mrbrem_data[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_data_%s", name));
           spec_mrbrem_cry[selIdx][varIdx]->SaveTo(d, TString::Format("mrbrem_mc_%s", name));
           spec_genie_fd[selIdx][varIdx]->SaveTo(d, TString::Format("genie_fd_%s", name));
      }
  }

  file->Close();
  return;

}