neutKEsyst.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SystShifts.h"

#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"

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

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Systs/EnergySysts.h"

#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TH2.h"
#include "TProfile.h"
#include "TSystem.h"

using namespace ana;

void neutKEsyst( bool IsFHC = false, bool IsFarDet = true)
{

  std::string sFHC_RHC = ( IsFHC    == true ? "fhc" : "rhc"  );
  std::string sFD_ND   = ( IsFarDet == true ? "FD"  : "ND"   );

  // --- Define my loader, and the file / dataset which it will load...
  std::string fname = "";
  if ( IsFHC ) {
    if (IsFarDet) {
      fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
//      fname = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
//      fname = "howard_neutronstudy_respin_birks_caf_fd_genie_nonswap_fhc_nova_v08_full_v1";
    } 
    else {
      fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
//      fname = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
//      fname = "howard_neutronstudy_respin_birks_caf_nd_genie_nonswap_fhc_nova_v08_full_v1";
    }
  } 
  else {
    if (IsFarDet) {
      fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
//      fname = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
//      fname = "howard_neutronstudy_respin_birks_caf_fd_genie_nonswap_rhc_nova_v08_full_v1";
    } 
    else {
      fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
//        fname = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
//      fname = "howard_neutronstudy_respin_birks_caf_nd_genie_nonswap_rhc_nova_v08_full_v1";
    }
  }

  SpectrumLoader loader(fname);
  loader.SetSpillCut(kStandardDQCuts);

  // weights, shifts
  const Var weight = kPPFXFluxCVWgt *  kXSecCVWgt2018;
  const ISyst* scale = &kNeutronVisEScalePrimariesSyst2018;

  // --- Define my output directory.
  std::string sOutFile = "Files_"+sFD_ND+"_"+sFHC_RHC+".root";

  // --- Define the number of GENIE interaction types I'm looking at....
  const std::vector<Cut>         GENIECuts = { kNoCut };
  const std::vector<std::string> GENIEStr  = { "All"  };
  //const std::vector<Cut>         GENIECuts = { kNoCut, kIsDytmanMEC, kIsDIS, kIsRes, kIsCoh, kIsQE };
  //const std::vector<std::string> GENIEStr  = { "All", "MEC", "DIS", "Resonance", "Coherent", "QuasiElactic" };
  const unsigned int NumGENIE = GENIECuts.size();

  // --- Define my broad cuts
  std::vector<std::string> CutNames; std::vector<Cut> kDetCuts;
  if (IsFarDet) {
    CutNames.push_back( "FD_2018full"   ); kDetCuts.push_back( kNumuQuality && kNumuContainFD2017 && kNumuPID2018 && kNumuCosmicRej2018 );
  } else {
    CutNames.push_back( "ND_2018full"   ); kDetCuts.push_back( kNumuQuality && kNumuContainND2017 && kNumuPID2018 );
  }
  unsigned int nDetCuts = kDetCuts.size();

  // --- Figure out where to load my quantile boundaries from.
  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__" + sFHC_RHC + "_full__numu2018.root";

  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile.c_str()).c_str() );
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
  // How many quantiles?
  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used 
  // Make my cuts.
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
  HadEFracQuantCuts.push_back( kNoCut );
  std::vector<std::string> QuantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
  const unsigned int nQuants = QuantNames.size();


  // Extra Vars
  // -----  Neutron multiplicity
  const Var kNNeutrons([](const caf::SRProxy* sr){
    if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
    return double(sr->mc.nu[0].nneutron);
  });

  // ----- Visible energy from neutron daughters
  const Var kVisEsum_Ndau([](const caf::SRProxy* sr){
    double visE = -5.0;
    if(sr->mc.nnu == 0) return visE;
    visE = 0.0;
    // loop over primary list
    for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
      // is it a neutron?
      if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
      visE+=sr->mc.nu[0].prim[PrimL].daughterVisEinslc;
    }
    return visE;
  });

  const MultiVar kVisE_Ndau([](const caf::SRProxy* sr){
    std::vector<double> visE;
    if(sr->mc.nnu == 0) return visE;
    for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL){
      if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
      visE.push_back(sr->mc.nu[0].prim[PrimL].daughterVisEinslc);
    }
    return visE; 
  });

  // ----- Number of 2D prongs
  const Var kNPng2D([](const caf::SRProxy* sr){
    double npng = -5.0;
    if(!sr->vtx.elastic.IsValid) return -5.0;
    npng = sr->vtx.elastic.fuzzyk.npng2d;
    return npng;
  });

  // ----- Number of 3D prongs
  const Var kNPng3D([](const caf::SRProxy* sr){
    double npng = -5.0;
    if(!sr->vtx.elastic.IsValid) return -5.0;
    npng = sr->vtx.elastic.fuzzyk.npng;
    return npng;
  });


  // ----- Cut -----
  const Cut kRemoveNoTruth([](const caf::SRProxy* sr) { 
      if (sr->mc.nnu == 0) return false; //if (sr->mc.nu.size() == 0) return false;
      return true;
    }); 


  // --- 1D Spectra
  // Energy plots
  Spectrum *NuE_nom  [NumGENIE][nDetCuts][nQuants];
  Spectrum *HadE_nom [NumGENIE][nDetCuts][nQuants]; 

  Spectrum *NuE_scale_p  [NumGENIE][nDetCuts][nQuants];
  Spectrum *HadE_scale_p [NumGENIE][nDetCuts][nQuants]; 

  Spectrum *NuE_scale_m  [NumGENIE][nDetCuts][nQuants];
  Spectrum *HadE_scale_m [NumGENIE][nDetCuts][nQuants]; 

  Spectrum *NNeu [NumGENIE][nDetCuts][nQuants];
  Spectrum *VisE_sum [NumGENIE][nDetCuts][nQuants];
  Spectrum *VisE [NumGENIE][nDetCuts][nQuants];

  // --- Define my Spectra which are to be filled from the loader
  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
    for (unsigned int det=0; det<nDetCuts; ++det) {
      for (unsigned int quant=0; quant<nQuants; ++quant) {
        // Define my new cut.
        Cut GENIECut = kRemoveNoTruth && GENIECuts[gen] && kDetCuts[det] && HadEFracQuantCuts[quant];
        std::string MyStr = " for GENIE Int. "+GENIEStr[gen]+", Cut "+CutNames[det]+", "+QuantNames[quant];
        // --- Now fill our 1D spectra!
        // Energy spectra
        NuE_nom [gen][det][quant] = new Spectrum("Reco #nu E (GeV), No Shift" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, kNoShift, weight);
        HadE_nom[gen][det][quant] = new Spectrum("Reco Had E (Gev), No Shift" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, kNoShift, weight);


        NuE_scale_p [gen][det][quant] = new Spectrum("Reco #nu E (GeV), Neutron Scale Up" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, SystShifts(scale, +1), weight);
        HadE_scale_p[gen][det][quant] = new Spectrum("Reco Had E (Gev), Neutron Scale Up" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, SystShifts(scale, +1), weight);

        NuE_scale_m [gen][det][quant] = new Spectrum("Reco #nu E (GeV), Neutron Scale Down" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, SystShifts(scale, -1), weight);
        HadE_scale_m[gen][det][quant] = new Spectrum("Reco Had E (Gev), Neutron Scale Down" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, SystShifts(scale, -1), weight);

        NNeu [gen][det][quant] = new Spectrum("Neutron Multiplicity" +MyStr, Binning::Simple(5,0,5), loader, kNNeutrons, GENIECut, kNoShift, weight); 
        VisE_sum [gen][det][quant] = new Spectrum("Summed visible energy from neutron daughters (GeV)" + MyStr, Binning::Simple(50, 0, 3), loader, kVisEsum_Ndau, GENIECut, kNoShift, weight);
        VisE [gen][det][quant] = new Spectrum("Visible energy from neutron daughters (GeV)" + MyStr, Binning::Simple(50, 0, 3), loader, kVisE_Ndau, GENIECut, kNoShift, weight);
  
      } // Loop through my Spectrum arrays.
    }
  }

  // --- Do it!
  loader.Go();

  // --- Open my outfile so that I save my histograms!
  TFile *OutFile = new TFile(sOutFile.c_str(),"RECREATE");

  // --- Save plots
  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
    for (unsigned int det=0; det<nDetCuts; ++det) {
      for (unsigned int quant=0; quant<nQuants; ++quant) {
        std::string MyStr = GENIEStr[gen]+"_"+CutNames[det]+"_"+QuantNames[quant];
        // --- Save my 1D histograms!
        // Energy spectra.
        NuE_nom   [gen][det][quant] -> SaveTo( OutFile,  TString("NuE_nom_")  + TString(MyStr) ) ;
        HadE_nom  [gen][det][quant] -> SaveTo( OutFile,  TString("HadE_nom_") + TString(MyStr) ) ;
        NuE_scale_p [gen][det][quant] -> SaveTo( OutFile,  TString("NuE_scale_p_")  + TString(MyStr) ) ;
        HadE_scale_p[gen][det][quant] -> SaveTo( OutFile,  TString("HadE_scale_p_") + TString(MyStr) ) ;
        NuE_scale_m [gen][det][quant] -> SaveTo( OutFile,  TString("NuE_scale_m_")  + TString(MyStr) ) ;
        HadE_scale_m[gen][det][quant] -> SaveTo( OutFile,  TString("HadE_scale_m_") + TString(MyStr) ) ;

        NNeu[gen][det][quant] -> SaveTo( OutFile,  TString("NNeu_") + TString(MyStr) ) ;
        VisE_sum[gen][det][quant] -> SaveTo( OutFile,  TString("VisE_sum_") + TString(MyStr) ) ;
        VisE[gen][det][quant] -> SaveTo( OutFile,  TString("VisE_") + TString(MyStr) ) ;
 
      }
    }
  }

}