numu_validation_numuvars.C

Go to the documentation of this file.

// Uses truth information. Only run on MC

/*                                                                                                                                           
*/

#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 "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Vars/GenieWeights.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();
}

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:

  UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut, const SystShifts &shift=kNoShift, const Var& wei = kTuftsWeightCC):
    fSelA(selA), 
    fAxis(tanAxis),
    fHist(loader, fAxis.fObj, fSelA.fObj, shift, wei),
    fName(tanAxis.fShortName + "-" + selA.fShortName ){};
  
  Selection fSelA;
  TangibleAxis fAxis; 
  Spectrum fHist;
  std::string fName;

};


void numu_validation_numuvars(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([](const caf::SRSpill* spill) {return spill->det;});

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

  //////////////////// Cuts ////////////////////////
  std::vector<Selection> selections_base;

  const Cut kNumuFD_cvn = kNumuQuality && kNumuContainFD && kCVNmSel && kNumuCosmicRej;
  const Cut kNumuND_cvn = kNumuQuality && kNumuContainND && kCVNmSel;

  if(isFD){
    selections_base.emplace_back(kNumuFD, "numuFD_remid",
                                 "Selected events pass the (remid) numu FD selection");
    selections_base.emplace_back(kNumuFD_cvn, "numuFD_cvn",
                                 "Selected events pass the (cvn) numu FD selection");
  }
  else if(!isFD){
    selections_base.emplace_back(kNumuND, "numuND_remid",
                                 "Selected events pass the (remid) numu ND selection");
    selections_base.emplace_back(kNumuND_cvn, "numuND_cvn",
                                 "Selected events pass the (cvn) numu ND selection");
  }

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

  // variables
  const Var kRecoMuE([](const caf::SRProxy* sr)
                     {
                       return (sr -> energy.numu.recomuonE);
                     });

  // Binning
  const Binning kXYBins = Binning::Simple(55,-2.20,2.20); //ND
  const Binning kZBins = Binning::Simple(32,  0,16.00); //ND
  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);

  std::vector<TangibleAxis> variables;
  variables.emplace_back(kNumuCCAxis, "numuE",
                         "CC energy estimator (GeV)  ");
  variables.emplace_back( HistAxis("Muon energy (GeV)", kEnergyBinning, kRecoMuE),
                          "muE", "Muon energy (GeV)");
  variables.emplace_back( HistAxis("Slice N_{Hit} (GeV)",  Binning::Simple(50, 0, 500), kNHit),
                         "slcNHit", "Number of hits in slice.  " );
  variables.emplace_back( HistAxis("Slice Calorimetric Energy (GeV)",  kEnergyBinning, kCaloE),
                         "calE", "Calorimetric energy of slice.  " );
  variables.emplace_back( HistAxis("Hadronic Energy (GeV)",  kHadronicEnergyBinning, kHadE),
                         "hadE",
                         "Hadronic enery, i.e. numu energy estimate minus muon track energy.  " );
  variables.emplace_back( HistAxis("Average Hadronic Energy Per Hit (GeV)",
                                   Binning::Simple(40, 0, 0.04), kHadEPerNHit),
                         "hadEPerNHit",
                         "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit.  " );
  variables.emplace_back( HistAxis("Track Energy Per Hit (GeV)",  Binning::Simple(40, 0, 0.04),
                                  kTrkEPerNHit),
                         "trkEPerNHit",
                         "Average energy per hit on primary track, i.e. trk_E/trk_n_hit.  " );
  variables.emplace_back( HistAxis("Hadronic N_{Hit} (GeV)",  Binning::Simple(50, 0, 100), kHadNHit),
                         "hadNHit", "Number of hits in hadronic cluster.  ");
  variables.emplace_back( HistAxis("Calorimetric Hadronic Energy (GeV)",
                                  kHadronicEnergyBinning, kNumuHadCalE),
                         "hadcalE", "Sum of calibrated energy deposits in hadronic cluster.  " );
  // variables.emplace_back( HistAxis("Slice Maximum Y (m)",  kXYBins, kMaxY), //not in decafs
  //                     "maxy", "Maximum Y position of slice.  " );
  variables.emplace_back( HistAxis("Number of Tracks in Slice",  Binning::Simple(14, 1, 15), kNKalman),
                         "nkal", "Number of tracks in slice.  " );
  variables.emplace_back( HistAxis("Number of Hits in Slice",  Binning::Simple(50, 0, 500), kNHit),
                         "nhit", "Number of hits in slice.  " );
  variables.emplace_back( HistAxis("Muon Track cos(#theta_{Z})",  Binning::Simple(50, 0,1), kDirZ),
                         "dirZ", "Z-direction of muon track.  " );
  variables.emplace_back( HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,3000), kSliceDuration),
                         "sliceDuration", "Slice duration.  " );
  variables.emplace_back( HistAxis("Number of Hits in Primary Track",  Binning::Simple(50,0,500), kTrkNhits),
                         "trkNhits", "Number of hits on primary track.  ");
  variables.emplace_back( HistAxis("Length of Primary Track (m)",  Binning::Simple(50,0,16), kTrkLength),
                         "trkLength", "Primary track length.  " );
  variables.emplace_back( HistAxis("Scattering Log-likelihood",  Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
                         "scatLL", "ReMId scattering log log-likelihood for primary track.  " );
  variables.emplace_back( HistAxis("dE/dx Log-likelihood",  Binning::Simple(50,-3,1), kReMIdDEDxLLH),
                         "dedxLL", "ReMId dE/dx log log-likelihood for primary track.  " );
  variables.emplace_back( HistAxis("Non-hadronic Plane Fraction",
                                  Binning::Simple(50,0,1), kReMIdMeasFrac),
                         "nonHadPlaneFrac", "ReMId Non-hadronic plane fraction.  " );
  variables.emplace_back( HistAxis("Muon ID",  kRemidBinning, kRemID),
                         "remid", "ReMId kNN score.  " );

  // detector specific variables:
  if(!isFD){
    variables.emplace_back( HistAxis("Track Start X Position (m)",  kXYBins, kTrkStartX),
                            "trkStartX", "Track start x position.  " );
    variables.emplace_back( HistAxis("Track Start Y Position (m)",  kXYBins, kTrkStartY),
                            "trkStartY", "Track start y position.  " );
    variables.emplace_back( HistAxis("Track Start Z Position (m)",  kZBins, kTrkStartZ),
                            "trkStartZ", "Track start z position.  " );
    variables.emplace_back( HistAxis("Track End X Position (m)",  kXYBins, kTrkEndX),
                            "trkEndX", "Track stop x position.  " );
    variables.emplace_back( HistAxis("Track End Y Position (m)",  kXYBins, kTrkEndY),
                            "trkEndY", "Track stop y position.  " );
    variables.emplace_back( HistAxis("Track End Z Position (m)",  kZBins, kTrkEndZ),
                            "trkEndZ", "Track stop z position.  " );
  }
  else if(isFD){
    variables.emplace_back( HistAxis("Track Start X Position (m)",  kFDXYBins, kTrkStartX),
                            "trkStartX", "Track start x position.  " );
    variables.emplace_back( HistAxis("Track Start Y Position (m)",  kFDXYBins, kTrkStartY),
                            "trkStartY", "Track start y position.  " );
    variables.emplace_back( HistAxis("Track Start Z Position (m)",  kFDZBins, kTrkStartZ),
                            "trkStartZ", "Track start z position.  " );
    variables.emplace_back( HistAxis("Track End X Position (m)",  kFDXYBins, kTrkEndX),
                            "trkEndX", "Track stop x position.  " );
    variables.emplace_back( HistAxis("Track End Y Position (m)",  kFDXYBins, kTrkEndY),
                            "trkEndY", "Track stop y position.  " );
    variables.emplace_back( HistAxis("Track End Z Position (m)",  kFDZBins, kTrkEndZ),
                            "trkEndZ", "Track stop z position.  " );
  }

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

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

  for(const auto& sel_base:selections_base){
    for(const auto& variable:variables){
      hists.emplace_back(sel_base, variable, 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;
}