NuMu2019_BasicPIDPlots_Spectrum.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 "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 "3FlavorAna/Vars/NumuVars.h"

#include <iostream>

using namespace ana;

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

void NuMu2019_BasicPIDPlots_Spectrum( bool isFHC, bool isData, bool isFD, std::string period = "full" ) {

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

  std::cout << "\n================================= \n"
            << "\n\t isFHC  --> " << isFHC  << " --> " << polarity
            << "\n\t isData --> " << isData << " --> " << data
            << "\n\t Period --> " << period
            << "\n\t isFD   --> " << 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"       ) );
  if (!isData) {
    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...
  if (isFD) {
    SelCuts.push_back( std::make_pair( kNumuQuality&&kNumuContainFD2017, "_QualCont" ) );
    SelCuts.push_back( std::make_pair( kNumuCutFD2018                  , "_FullCut"  ) );
  } else {
    SelCuts.push_back( std::make_pair( kNumuQuality&&kNumuContainND2017, "_QualCont" ) );
    SelCuts.push_back( std::make_pair( kNumuCutND2018                  , "_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?
  Var         kMyWeight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
  if (isData) kMyWeight = kUnweighted;

  // Figure out my definition, and then declare my spectrum loader.
  std::string  MyDef = "";
  if      ( isFHC &&  isFD  ) MyDef ="prod_sumdecaf_R17-11-14-prod4reco.b_fd_cosmic_fhcTune_full_v1_goodruns_numu2018loose";
  else if ( isFHC && !isFD  ) MyDef ="prod_sumdecaf_R17-11-14-prod4reco.a_fd_cosmic_rhcTune_HighGain_v1_goodruns_snapshot_170116_numu2018loose";
  else if ( isFHC && !isData) MyDef ="prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1";
  else if ( isFHC &&  isData) MyDef ="prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
  else if (!isFHC && !isData) MyDef ="prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_rhc_nova_v08_full_v1";
  else if (!isFHC &&  isData) MyDef ="prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_numi_rhc_full_v1_goodruns";
  SpectrumLoader loader( MyDef );
  loader.SetSpillCut   ( kStandardSpillCuts );
  
  // Declare my spectra
  Spectrum *sNuRecoE       [NCuts];
  Spectrum *sReMIdScore    [NCuts];
  Spectrum *sProngScoremuon[NCuts];
  Spectrum *sCosRejScoreP4 [NCuts];
  Spectrum *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;
    if (isFD) kThisCut = kThisCut && kInCosmicTimingWindow;

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