NuMu2020_BasicPIDPlots_FD.C

Go to the documentation of this file.

#include "CAFAna/Prediction/PredictionNoExtrap.h"

#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

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

#include "CAFAna/Vars/CVNProngVars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "3FlavorAna/Ana2020/NuMu_PIDOptimisation/NuMu2020_BasicPIDPlots_header.h"

#include "TH2.h"
#include "TFile.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include <iostream>

using namespace ana;


void NuMu2020_BasicPIDPlots_FD( bool isFHC, std::string period = "full" ) {

  std::string polarity = "fhc"; if (!isFHC) polarity = "rhc";
  std::string data     = "mc" ;
  std::string detector = "fd" ;

  std::cout << "\n================================= \n"
            << "\n\t isFHC  --> " << isFHC << " --> " << polarity
            << "\n\t Period --> " << period
            << "\n\t Far detector isFD --> " << detector
            << "\n================================= \n"
            << std::endl;
  
  // Make my vector of cuts.
  std::vector<std::pair<Cut, std::string> > MyCuts = BasicPIDPlots_Cuts( isFHC, true, false );
  const size_t NCuts = MyCuts.size();

  // What XSec and PPFX weights am I using?
  const Var kMyWeight = kXSecCVWgt2020 * kPPFXFluxCVWgt;
  
  // Figure out my definition, and then declare my spectrum loader.
  std::string MyNonSwapDef  = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_" +polarity+"_nova_v08_"+period+"_v1";
  std::string MyFluxSwapDef = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_"+polarity+"_nova_v08_"+period+"_v1";
  std::string MyTauSwapDef  = "prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_tau_"     +polarity+"_nova_v08_"+period+"_v1";
  Loaders *loader = new Loaders();
  loader->SetLoaderPath( MyNonSwapDef , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap  );
  loader->SetLoaderPath( MyFluxSwapDef, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap );
  loader->SetLoaderPath( MyTauSwapDef , caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap  );
  loader->SetSpillCut  ( kStandardSpillCuts );

  // Declare my spectra
  PredictionNoExtrap *sNuRecoE       [NCuts];
  PredictionNoExtrap *sReMIdScore    [NCuts];
  PredictionNoExtrap *sProngScoremuon[NCuts];
  PredictionNoExtrap *sCosRejScoreP4 [NCuts];
  PredictionNoExtrap *sCosRejScoreP5 [NCuts];
  PredictionNoExtrap *sOldPreselmuon [NCuts];
  PredictionNoExtrap *sOldPreselcosm [NCuts];
  PredictionNoExtrap *sLoosePTPmuon  [NCuts];
  PredictionNoExtrap *sLoosePTPcosm  [NCuts];

  for (size_t cc=0; cc < NCuts; ++cc) {
    Cut kThisCut = MyCuts[cc].first;

    sNuRecoE       [cc] = new PredictionNoExtrap( *loader, "Reco Nu E"            , kOptEnBins   , kCCE             , kThisCut, kNoShift, kMyWeight );
    /*
    sReMIdScore    [cc] = new PredictionNoExtrap( *loader, "ReMId Score"          , kRemidBinning, kRemID           , kThisCut, kNoShift, kMyWeight );
    sProngScoremuon[cc] = new PredictionNoExtrap( *loader, "Prong CVN Muon Score" , kSimpBins    , kCVNBestMuonScore, kThisCut, kNoShift, kMyWeight );
    sCosRejScoreP4 [cc] = new PredictionNoExtrap( *loader, "Prod4 CosRej Score"   , kSimpBins    , kNumuContPID2019 , kThisCut, kNoShift, kMyWeight );
    sCosRejScoreP5 [cc] = new PredictionNoExtrap( *loader, "Prod5 CosRej Score"   , kSimpBins    , kNumuContPID     , kThisCut, kNoShift, kMyWeight );
    sOldPreselmuon [cc] = new PredictionNoExtrap( *loader, "Old Presel Muon Score", kSimpBins    , kCVNm_oldpresel  , kThisCut, kNoShift, kMyWeight );
    sOldPreselcosm [cc] = new PredictionNoExtrap( *loader, "Old Presel Cosm Score", kSimpBins    , kCVNcos_oldpresel, kThisCut, kNoShift, kMyWeight );
    sLoosePTPmuon  [cc] = new PredictionNoExtrap( *loader, "Loose Pt/p Muon Score", kSimpBins    , kCVNm_looseptp   , kThisCut, kNoShift, kMyWeight );
    sLoosePTPcosm  [cc] = new PredictionNoExtrap( *loader, "Loose Pt/p Cosm Score", kSimpBins    , kCVNcos_looseptp , kThisCut, kNoShift, kMyWeight );
    */
  }

  // Set my loader going.
  loader->Go();

  SpectrumLoader SimpLoad(MyNonSwapDef.c_str());
  SimpLoad.SetSpillCut(kStandardSpillCuts);
  Spectrum* NoCutRecoNu = new Spectrum("NoCutRecoNu", kOptEnBins, SimpLoad, kCCE, kNoCut, kNoShift, kMyWeight );
  SimpLoad.Go();

  // Make my Output file
  std::string OutName = "BasicPIDPlots_2020_"+detector+"_"+data+"_"+polarity+"_"+period+".root";
  TFile *OutFile = TFile::Open(OutName.c_str(), "RECREATE");
  for (size_t cc=0; cc < NCuts; ++cc) {
    std::string ThisCutNa = MyCuts[cc].second;

    sNuRecoE       [cc] -> SaveTo( OutFile,  TString("RecoNuE"           )+TString(ThisCutNa) ) ;
    /*
    sReMIdScore    [cc] -> SaveTo( OutFile,  TString("ReMIdScore"        )+TString(ThisCutNa) ) ;
    sProngScoremuon[cc] -> SaveTo( OutFile,  TString("ProngCVNMuonScore" )+TString(ThisCutNa) ) ;
    sCosRejScoreP4 [cc] -> SaveTo( OutFile,  TString("Prod4CosRejScore"  )+TString(ThisCutNa) ) ;
    sCosRejScoreP5 [cc] -> SaveTo( OutFile,  TString("Prod5CosRejScore"  )+TString(ThisCutNa) ) ;
    sOldPreselmuon [cc] -> SaveTo( OutFile,  TString("OldPreselMuonScore")+TString(ThisCutNa) ) ;
    sOldPreselcosm [cc] -> SaveTo( OutFile,  TString("OldPreselCosmScore")+TString(ThisCutNa) ) ;
    sLoosePTPmuon  [cc] -> SaveTo( OutFile,  TString("LoosePtpMuonScore" )+TString(ThisCutNa) ) ;
    sLoosePTPcosm  [cc] -> SaveTo( OutFile,  TString("LoosePtpCosmScore" )+TString(ThisCutNa) ) ;
    */
  }
  NoCutRecoNu -> SaveTo( OutFile,  TString("NoCutRecoNu" ) ) ;
  
  OutFile -> Close();
}