caf_numu_reco_minus_true.C

Go to the documentation of this file.

// Uses truth information. Only run on MC

/*                                                                                                                                           
*/
#ifdef __CINT__
void caf_numu_reco_minus_true(std::string kInputFileName, std::string kOutputFileName, bool isFD)
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

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

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/NumuVars.h"

#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"

// #include "NumuEnergy/NumuEAlg.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TLine.h"
#include "TFile.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLine.h"
#include "TLatex.h"
#include "THStack.h"

using namespace ana;

// Put a "NOvA Preliminary" tag in the corner
void Preliminary(){
  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
  prelim->SetTextColor(kBlue);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
                       [](const caf::StandardRecord* sr)
                       {
                         if(sr->mc.nnu == 0) return false;
                         assert(sr->mc.nnu == 1);
                         return (sr->mc.nu[0].mode == caf::kMEC);
                       });

template<class T> class Tangible{

public:
  Tangible(const T& obj, const std::string&  shortName,
                                  const std::string&  blurb ):
  fObj(obj),
  fShortName(shortName),
  fBlurb(blurb)
  {};

  T           fObj;
  std::string fBlurb;
  std::string fShortName;

};


typedef Tangible<Cut> Selection;
typedef Tangible<HistAxis> TangibleAxis;


class UsefulHist{

public:

  // what's wrong with doing the following?
  // UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut):
  //   fSelA(selA), 
  //   fAxis(tanAxis),
  //   fHist(loader, fAxis.fObj, fSelA.fObj),
  //   fName(tanAxis.fShortName + "-" + selA.fShortName ){};

  UsefulHist(Selection selA, Selection selB, Selection selC, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut):
    fSelA(selA), 
    fSelB(selB), 
    fSelC(selC), 
    fAxis(tanAxis),
    fHist(loader, fAxis.fObj, fSelA.fObj && fSelB.fObj && fSelC.fObj && !kIsDytmanMEC),
    //fName(tanAxis.fShortName + "-" + selA.fShortName + "_" + selB.fShortName + "_" + selC.fShortName ){};
    fName(tanAxis.fShortName + "-" + selA.fShortName + "_" + selB.fShortName + "_" + selC.fShortName ){};
  
  Selection fSelA;
  Selection fSelB;
  Selection fSelC;
  TangibleAxis fAxis; 
  Spectrum fHist;
  std::string fName;

};


void caf_numu_reco_minus_true(std::string kInputFileName, std::string kOutputFileName, bool isFD){

  std::cout << "\nrun  : ---- Running CAF NuMU ND Validation.\n";
  std::cout << "run  :     input file name:       " << kInputFileName   << std::endl;
  std::cout << "run  :    output file name:       " << kOutputFileName  << std::endl;
  std::cout << "Are we using FD binning/cuts?     " << isFD  << std::endl;

  // make loader with user input dataset
  SpectrumLoader loader(kInputFileName);

  // make spill histogram tokeep record of spills. Going to use to count total number of events
  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});

  // set spill cuts
  loader.SetSpillCut(kStandardSpillCuts);
  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
  

  //////////////////// Cuts ////////////////////////
  // do cuts for contained and interaction type and true numu CC

  // kIsNumuCC
  // kIsQE
  // !kIsQE
  // kNumuFD and kNumuND. Use ifFD boolean to switch
  // kTruthContainedND. Only apply if ifFD==0
  // !kTruthContainedND. Only apply if ifFD==0  

  std::vector<Selection> selections_base;
  std::vector<Selection> selections_int;
  std::vector<Selection> selections_contain;

  //base
  selections_base.emplace_back(kIsNumuCC, "true_numucc",
                          "Selected events are truly numu cc");
  if(isFD) selections_base.emplace_back(kNumuFD, "numuFD",
                                        "Selected events pass the numu FD selection");
  else selections_base.emplace_back(kNumuND, "numuND",
                                    "Selected events pass the numu ND selection");
  //int.
  selections_int.emplace_back(kNoCut, "no_cut",
                          "No selection applied");
  selections_int.emplace_back(kIsQE, "true_qe",
                          "Selected events are truly quasi-elastic");
  selections_int.emplace_back(!kIsQE, "true_non_qe",
                          "Selected events are truly non quasi-elastic");
  //contain
  selections_contain.emplace_back(kNoCut, "no_cut",
                                  "No selection applied");

  selections_contain.emplace_back(kTruthContainedND, "true_nd_contained",
                                  "Selected events are truly contained in the ND");
  selections_contain.emplace_back(!kTruthContainedND, "true_nd_uncontained",
                                  "Selected events are truly un-contained in the ND");

  ////////////////////////////////////////////////////////////////////////////////


  //////////////////// Vars ////////////////////////
  const Var kTrueNHitSlc({"mc.nu", "mc.nnu", "mc.nu.nhitslc"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nnu == 0) ? 0 : sr->mc.nu[0].nhitslc;});
  const Var kTrueNHitTot({"mc.nu", "mc.nnu", "mc.nu.nhittot"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nnu == 0) ? 0 : sr->mc.nu[0].nhittot;});
  const Var kRecoNHitTot({"slc.nhit"},
                         [](const caf::StandardRecord* sr)
                         {return sr->slc.nhit;});
  const Var kTruehadE({"mc.nu", "mc.nnu", "mc.nu.*"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nnu == 0 || sr->mc.nu[0].prim.size() == 0) ? 0 : 
                             (sr->mc.nu[0].visE - sr->mc.nu[0].prim[0].visE);
                         });
  const Var kTruehadESlc({"mc.nu", "mc.nnu", "mc.nu.*"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nnu == 0 || sr->mc.nu[0].prim.size() == 0) ? 0 : 
                             (sr->mc.nu[0].visEinslc - sr->mc.nu[0].prim[0].visEinslc);
                         });
  const Var kRecohadE({"energy.numusimp.hadcalE"},
                         [](const caf::StandardRecord* sr)
                         {return (sr -> energy.numusimp.hadcalE);
                         });
  const Var kTrueMuE({"mc.nu", "mc.nnu", "mc.nu.*"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nu[0].prim.size() == 0) ? 0 : 
                             (sr->mc.nu[0].prim[0].p.E());
                         });
  const Var kTrueMuESlc({"mc.nu", "mc.nnu", "mc.nu.*"},
                         [](const caf::StandardRecord* sr)
                         {return (sr->mc.nu[0].prim.size() == 0) ? 0 : 
                             sr->mc.nu[0].prim[0].p.E();
                         });
  const Var kRecoMuE({"energy.mutrkE.E","sel.remid.bestidx"},
                         [](const caf::StandardRecord* sr)
                         {
                            if(sr->sel.remid.bestidx < 0.f) return 0.f;
                            if(sr->sel.remid.bestidx == 999) return 0.f;
                            return (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
                         });
  const Var kTrueNuEnergy({"mc.nu", "mc.nnu", "mc.nu.E"},
                          [](const caf::StandardRecord* sr)
                          {return (sr->mc.nnu == 0) ? 0 :
                              sr->mc.nu[0].E;
                          });
  const Var kRecoNuEnergy({"energy.numusimp.trkccE"},
                          [](const caf::StandardRecord* sr)
                          {return (sr->mc.nnu == 0) ? 0 :
                              sr->energy.numusimp.trkccE;
                          });


  // reco - true / true
  const Var kRMTOT_NHitSlc({"mc.nu", "mc.nnu", "mc.nu.nhitslc", "slc.nhit"},
                           [](const caf::StandardRecord* sr)
                           {
                             if(sr->mc.nnu == 0) return -10.;
                             
                             double nhitReco = (double)sr->slc.nhit;
                             double nhitTrue = (double)sr->mc.nu[0].nhitslc;
                             return ((nhitReco - nhitTrue) / nhitTrue);

                           });

  
  const Var kRMTOT_Energy({"mc.nu", "mc.nnu", "mc.nu.E", "energy.numusimp.trkccE"},
                          [](const caf::StandardRecord* sr)
                          {return (sr->mc.nnu == 0) ? -10 : 
                              ( (sr->energy.numusimp.trkccE - sr->mc.nu[0].E) / (sr->mc.nu[0].E) );
                          });
  //ehad
  const Var kRMTOT_hadE({"mc.nu", "mc.nnu", "mc.nu.*", "energy.numusimp.hadcalE"},
                          [](const caf::StandardRecord* sr)
                          {
                            if(sr->mc.nnu == 0) return -10.f;
                            if(sr->mc.nu[0].prim.size() == 0) return -10.f;
                            float trueHadE = (sr -> mc.nu[0].E) - (sr -> mc.nu[0].prim[0].p.E());
                            float muE=sr->energy.mutrkE[sr->sel.remid.bestidx].E;
                            float recoHadE=sr->energy.numusimp.trkccE-muE;

                            if(trueHadE == 0) return -10.f;
                            return ( (recoHadE - trueHadE) / trueHadE );
                          });
  // E_mu
  const Var kRMTOT_Emu({"mc.nu", "mc.nnu", "mc.nu.*", "energy.mutrkE.E", "sel.remid.bestidx"},
                          [](const caf::StandardRecord* sr)
                          {
                            if(sr->mc.nnu == 0) return -10.f;
                            if(sr->mc.nu[0].prim.size() == 0) return -10.f;
                            if(sr->sel.remid.bestidx < 0.f) return -10.f;
                            if(sr->sel.remid.bestidx == 999) return -10.f;
                            float trueMuE = (sr -> mc.nu[0].prim[0].p.E());
                            float recoMuE = (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
                            if(trueMuE == 0) return -10.f;
                            return ( (recoMuE - trueMuE) / trueMuE );
                          });

  ////////////////////////////////////////////////////////////////////////////////

  // Binning
  const Binning kXYBins = Binning::Simple(55,-2.20,2.20);
  const Binning kZBins = Binning::Simple(32,  0,16.00);
  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
  const Binning kFDZBins = Binning::Simple(55,0,60.00);
  const Binning kEnergyBinning = Binning::Simple(50,0,5);
  const Binning kHadronicEnergyBinning = Binning::Simple(50,0,5);

  //HistAxis
  const HistAxis axRMTOT_NHitSlc("nhit in slc (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_NHitSlc);
  const HistAxis axRMTOT_Energy("Neutrino energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_Energy);
  const HistAxis axRMTOT_HadE("Hadronic energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_hadE);
  const HistAxis axRMTOT_Emu("Muon energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_Emu);
  
  const HistAxis axTrueNHit("NHits", Binning::Simple(100,0,500), kTrueNHitTot);
  const HistAxis axTrueNHitSlc("NHits", Binning::Simple(100,0,500), kTrueNHitSlc);
  const HistAxis axRecoNHit("NHits", Binning::Simple(100,0,500), kRecoNHitTot);
  const HistAxis axTrueHadE("Had E. (GeV)", kHadronicEnergyBinning, kTruehadE);
  const HistAxis axTrueHadESlc("Had E. (GeV)", kHadronicEnergyBinning, kTruehadESlc);
  const HistAxis axRecoHadE("Had E. (GeV)", kHadronicEnergyBinning, kRecohadE);
  const HistAxis axTrueMuE("Muon E. (GeV)", kEnergyBinning, kTrueMuE);
  const HistAxis axTrueMuESlc("Muon E. (GeV)", kEnergyBinning, kTrueMuESlc);
  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
  const HistAxis axTrueNuEnergy("Neutrino E. (GeV)", kEnergyBinning, kTrueNuEnergy);
  const HistAxis axRecoNuEnergy("Neutrino E. (GeV)", kEnergyBinning, kRecoNuEnergy);
  
  ////////////////////////////////////////////////////////////////////////////////

  std::vector<TangibleAxis> variables;
  variables.emplace_back(axRMTOT_NHitSlc, "res-nhit_slc_reco_minus_true_over_true",
                         "Number of hits in slc. (reco-true)/true");
  variables.emplace_back(axRMTOT_Energy, "res-neutrino_energy_reco_minus_true_over_true",
                         "Neutrino energy (reco-true)/true");
  variables.emplace_back(axRMTOT_HadE, "res-hadronic_energy_reco_minus_true_over_true",
                         "Hadronic energy (reco-true)/true");
  variables.emplace_back(axRMTOT_Emu, "res-muon_energy_reco_minus_true_over_true",
                         "Muon energy (reco-true)/true");

  std::vector<TangibleAxis> dist_variables;
  // include these when the variables are corrected
  // dist_variables.emplace_back(axTrueNHit, "true-nhit_truth", "NHits");
  // dist_variables.emplace_back(axTrueNHitSlc, "true-nhit_slc_truth", "NHits");
  // dist_variables.emplace_back(axRecoNHit, "true-nhit_reco", "NHits");
  // dist_variables.emplace_back(axTrueHadE, "true-hade_truth", "HadE (GeV)");
  // dist_variables.emplace_back(axTrueHadESlc, "true-hade_slc_truth", "HadE (GeV)");
  // dist_variables.emplace_back(axRecoHadE, "true-hade_reco", "HadE (GeV)");
  // dist_variables.emplace_back(axTrueMuE, "true-mue_truth", "Mu E (GeV)");
  // dist_variables.emplace_back(axTrueMuESlc, "true-mue_slc_truth", "Mu E (GeV)");
  // dist_variables.emplace_back(axRecoMuE, "true-mue_reco", "Mu E (GeV)");
  // dist_variables.emplace_back(axTrueNuEnergy, "true-nu_energy_truth", "Nu E (GeV)");
  // dist_variables.emplace_back(axRecoNuEnergy, "true-nu_energy_reco", "Nu E (GeV)");

  ////////////////////////////////////////////////////////////////////////////////

  std::vector<UsefulHist> hists;
  hists.reserve(selections_base.size() * selections_int.size() * selections_contain.size() * variables.size() + selections_base.size() * selections_int.size() * selections_contain.size() * dist_variables.size());

  
  if(!isFD){
    for(const auto& sel_base:selections_base){
      for(const auto& sel_int:selections_int){
        for(const auto& sel_contain:selections_contain){
          for(const auto& variable:variables){
            hists.emplace_back(sel_base, sel_int, sel_contain, variable, loader);
          }
          // include these when the variables are corrected
          // for(const auto& dist:dist_variables){
          //   hists.emplace_back(sel_base, sel_int, sel_contain, dist, loader);
          // }
        }
      }
    }
  }

  if(isFD){
    for(const auto& sel_base:selections_base){
      for(const auto& sel_int:selections_int){
        for(const auto& variable:variables){
          hists.emplace_back(sel_base, sel_int, selections_contain[0], variable, loader);
        }
        // for(const auto& dist:dist_variables){
        //   hists.emplace_back(sel_base, sel_int, selections_contain[0], dist, loader);
        // }
      }
    }
  }

  std::cout << "\nrun  : --- run loader.\n";
  loader.Go();
  std::cout << "\nrun  : --- done.\n";
  double pot = hists[0].fHist.POT();

  std::cout << "\nrun  : --- save output.\n";
  // TFile  inputfile(kInputFileName.c_str(),     "READ");
  TFile outputfile(kOutputFileName.c_str(), "RECREATE");

  for(const auto& hist:hists){
    TH1 *temp = hist.fHist.ToTH1(pot);
    std::string tempName = hist.fName;
    temp->Write(tempName.c_str());
  }

  std::cout << "\nrun  : --- save PoT.\n";

  TH1F *pothist = new TH1F("meta-TotalPOT",    "TotalPOT", 1, 0, 1);
  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
  pothist->SetBinContent(1, pot);
  evthist->SetBinContent(1, spillHist->GetEntries());
  pothist->Write("meta-TotalPOT");  
  evthist->Write("meta-TotalEvents");

  std::cout << "\nrun  : --- close output files.\n";
  outputfile.Close();
  // inputfile.Close(); // will produce error if running over SAM dataset
  std::cout << "\nrun  : --- done.\n";

  // dictionary to translate histogram names 
  /*

  */

  return;
}

#endif