CutFlow_MC.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Prod3Loaders.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/ReweightableSpectrum.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/Extrap/ModularExtrap.h"

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

#include "OscLib/OscCalcPMNSOpt.h"

#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"

using namespace ana;

#include <fstream>

struct Plot {
  std::string label;
  std::string fname;
  Binning bins;
  Var var;
};

void CutFlow_MC()
{
  // --- Set my loader
  Loaders *loader = new Prod3NomLoaders( ECAFType::kFullCAF, "full");
  loader->SetSpillCut(kStandardSpillCuts);
  
  // --- Set my oscillation calculator
  osc::OscCalcPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc);

  // --- Define my cuts and their names.
  std::vector<Cut> TieredCuts;
  std::vector<std::string> CutNames;
  CutNames.emplace_back("NoCut");       TieredCuts.emplace_back(kNoCut);
  CutNames.emplace_back("Quality");     TieredCuts.emplace_back(kNumuQuality);
  CutNames.emplace_back("Containment"); TieredCuts.emplace_back(kNumuQuality && kNumuContainFD2017);
  CutNames.emplace_back("Cosmic_Rej");  TieredCuts.emplace_back(kNumuQuality && kNumuContainFD2017 && kNumuCosmicRej2017);
  CutNames.emplace_back("Particle_ID"); TieredCuts.emplace_back(kNumuQuality && kNumuContainFD2017 && kNumuCosmicRej2017 && kNumuPID2017 );

  size_t NCuts = TieredCuts.size();

  // --- Define my objects.
  Plot plot = {"CC Energy [GeV]", "NuMuCCEnergy", kNumuCCEOptimisedBinning, kCCE};
  PredictionNoExtrap* ReconEnergy[NCuts];
  
  // --- Set my objects.
  for(size_t cut = 0; cut < NCuts; ++cut) {
    const Cut ThisCut = TieredCuts[cut];    
    ReconEnergy[cut] = new PredictionNoExtrap( *loader, plot.label, plot.bins, plot.var, ThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
  }

  // --- Set the loader off!
  loader->Go();

  // --- Define an outfile.
  std::string OutName = "MonteCarloCutFlow.root";
  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
  OutFile  -> cd();
  for(size_t cut = 0; cut < NCuts; ++cut){
    ReconEnergy[cut] -> Predict(&calc).SaveTo( OutFile,  TString(CutNames[cut])+TString("_FullMC"  ) ) ;
    ReconEnergy[cut] -> PredictComponent(&calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).SaveTo( OutFile,  TString(CutNames[cut])+TString("_NuMuCC"   ) ) ;
    ReconEnergy[cut] -> PredictComponent(&calc, Flavors::kAllNuE    , Current::kCC, Sign::kBoth).SaveTo( OutFile,  TString(CutNames[cut])+TString("_AllNue"   ) ) ;
    ReconEnergy[cut] -> PredictComponent(&calc, Flavors::kNuEToNuMu , Current::kCC, Sign::kBoth).SaveTo( OutFile,  TString(CutNames[cut])+TString("_NuEToMu"  ) ) ;
    ReconEnergy[cut] -> PredictComponent(&calc, Flavors::kAll       , Current::kNC, Sign::kBoth).SaveTo( OutFile,  TString(CutNames[cut])+TString("_AllNC"    ) ) ;
    ReconEnergy[cut] -> PredictComponent(&calc, Flavors::kAllNuTau  , Current::kCC, Sign::kBoth).SaveTo( OutFile,  TString(CutNames[cut])+TString("_AllTau"   ) ) ;
  }

} // End of function