NuMu2020_BasicPIDPlots_ND.C

Go to the documentation of this file.

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

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Var.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_ND( bool isFHC, bool isData, std::string period = "full" ) {

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

  std::cout << "\n================================= \n"
            << "\n\t isFHC  --> " << isFHC  << " --> " << polarity
            << "\n\t isData --> " << isData << " --> " << data
            << "\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, false, isData );
  const size_t NCuts = MyCuts.size();

  // What XSec and PPFX weights am I using?
  Var         kMyWeight = kXSecCVWgt2020 * kPPFXFluxCVWgt;
  if (isData) kMyWeight = kUnweighted;

  // Figure out my definition, and then declare my spectrum loader.
  std::string  MyDef = "prod_caf_R19-11-18-prod5reco.d_nd_numi_"+polarity+"_full_v1";
  if (!isData) MyDef = "prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_"+polarity+"_nova_v08_full_v1";
  SpectrumLoader loader( MyDef );
  loader.SetSpillCut   ( kStandardSpillCuts );
  
  // Declare my spectra
  Spectrum *sNuRecoE       [NCuts];
  Spectrum *sReMIdScore    [NCuts];
  Spectrum *sProngScoremuon[NCuts];
  Spectrum *sCosRejScoreP4 [NCuts];
  Spectrum *sCosRejScoreP5 [NCuts];
  Spectrum *sOldPreselmuon [NCuts];
  Spectrum *sOldPreselcosm [NCuts];
  Spectrum *sLoosePTPmuon  [NCuts];
  Spectrum *sLoosePTPcosm  [NCuts];

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

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

  // Set my loader going.
  loader.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) ) ;
    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) ) ;
    if (isFD) {
      sCosRejScoreP4[cc]-> SaveTo( OutFile,  TString("Prod4CosRejScore"  )+TString(ThisCutNa) ) ;
      sCosRejScoreP5[cc]-> SaveTo( OutFile,  TString("Prod5CosRejScore"  )+TString(ThisCutNa) ) ;
    }
  }
  OutFile -> Close();
}