NuMu2019_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/SpectrumLoader.h"

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

#include "TH2.h"
#include "TFile.h"

#include <iostream>

using namespace ana;

// ---------------------------------------------- //
// This has to be ran out of R19-04-17-2019ana.a  //
// ---------------------------------------------- //

void NuMu2019_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> > TrueCuts, SelCuts, MyCuts;
  TrueCuts.clear(); SelCuts.clear(); MyCuts.clear();
  // Create my true cuts...
  TrueCuts.push_back( std::make_pair( kNoCut                , "_NoTrue"       ) );
  TrueCuts.push_back( std::make_pair( kHasNeutrino          , "_TrueNu"       ) );
  TrueCuts.push_back( std::make_pair( kIsNumuCC && kIsNu    , "_TrueNuMu"     ) );
  TrueCuts.push_back( std::make_pair( kIsNumuCC && kIsAntiNu, "_TrueAntiNuMu" ) );
  // Create my selection cuts...
  SelCuts.push_back( std::make_pair( kNumuQuality&&kNumuContainFD2017, "_QualCont" ) );
  SelCuts.push_back( std::make_pair( kNumuCutFD2018                  , "_FullCut"  ) );
  // Put them together...
  for (size_t sel=0; sel<SelCuts.size(); ++sel) {
    for (size_t tr=0; tr<TrueCuts.size(); ++tr) {
      Cut kThisCut          = SelCuts[sel].first && TrueCuts[tr].first;
      std::string kThisName = SelCuts[sel].second + TrueCuts[tr].second;
      MyCuts.push_back( std::make_pair( kThisCut, kThisName ) );
    }
  }
  const size_t NCuts = MyCuts.size();

  // What XSec and PPFX weights am I using?
  const Var kMyWeight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
  
  // Figure out my definition, and then declare my spectrum loader.
  std::string sTag = ".d";
  if (!isFHC) sTag = ".e";
  std::string MyNonSwapDef  = "prod_caf_R17-11-14-prod4reco"+sTag+"_fd_genie_nonswap_" +polarity+"_nova_v08_full_v1";
  std::string MyFluxSwapDef = "prod_caf_R17-11-14-prod4reco"+sTag+"_fd_genie_fluxswap_"+polarity+"_nova_v08_full_v1";
  std::string MyTauSwapDef  = "prod_caf_R17-11-14-prod4reco"+sTag+"_fd_genie_tau_"     +polarity+"_nova_v08_full_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 *sCVNScoremuon  [NCuts];

  // Declare some binnings.
  const Binning kOptEnBins = kNumuCCEOptimisedBinning;
  const Binning kSimpBins  = Binning::Simple(102, -0.01, 1.01 );

  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    , kNumuContPID     , kThisCut, kNoShift, kMyWeight );
    sCVNScoremuon  [cc] = new PredictionNoExtrap( *loader, "CVN Muon Score"       , kSimpBins    , kCVNm            , kThisCut, kNoShift, kMyWeight );
  }

  // Set my loader going.
  loader->Go();
  
  // Make my Output file
  std::string OutName = "BasicPIDPlots_2019_"+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) ) ;
    sCVNScoremuon  [cc] -> SaveTo( OutFile,  TString("CVNMuonScore"      )+TString(ThisCutNa) ) ;
  }
  OutFile -> Close();
}