caf_numu_vars.C

Go to the documentation of this file.

// Uses truth information. Only run on MC

/*                                                                                                                                           
*/
#ifdef __CINT__
void caf_numu_vars(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/Vars/HistAxes.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;

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);
                       });

class UsefulHist{

public:

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

};


void caf_numu_vars(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 ////////////////////////
  std::vector<Selection> selections_base;

  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");
  ////////////////////////////////////////////////////////////////////////////////

  // variables
  const Var kRecoMuE({"energy.mutrkE.E","sel.remid.bestidx","energy.mutrkE"},
                     [](const caf::StandardRecord* sr)
                     {
                       if(sr->sel.remid.bestidx < 0.f) return 0.f;
                       if(sr->sel.remid.bestidx == 999) return 0.f;
                       if(sr->energy.mutrkE.size() == 0) return 0.f;
                       return (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
                     });


  // 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);

  // Axes
  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
  /*
    Variables not included yet:
    kSliceTime, 
    kCosBackDist, kCosFwdDist, kKalmanFwdDist, kKalmanBackDist, 
    kNKalman, kNumuCosRejAngleKal, kNumuContPID, kNplanesToFront, kNplanesToBack, kNcellsFromEdge
   */

  std::vector<TangibleAxis> variables;
  variables.emplace_back(kNumuCCShiftAxis, "reco-numuEShift",
                         "Shifted CC energy estimator (GeV)  ");
  variables.emplace_back(kNumuCCAxis, "reco-numuE",
                         "CC energy estimator (GeV)  ");
  variables.emplace_back(axRecoMuE, "reco-muE",
                         "Muon energy (GeV)");
  variables.emplace_back( HistAxis("Slice N_{Hit} (GeV)",  Binning::Simple(50, 0, 500), kNHit),
                         "reco-slcNHit", "Number of hits in slice.  " );
  variables.emplace_back( HistAxis("Slice Calorimetric Energy (GeV)",  kEnergyBinning, kCaloE),
                         "reco-calE", "Calorimetric energy of slice.  " );
  variables.emplace_back( HistAxis("Hadronic Energy (GeV)",  kHadronicEnergyBinning, kHadE),
                         "reco-hadE",
                         "Hadronic enery, i.e. numu energy estimate minus muon track energy.  " );
  variables.emplace_back( HistAxis("Hadronic Energy (GeV)",  kHadronicEnergyBinning, kHadEshift),
                         "reco-hadEShift",
                         "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),
                         "reco-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),
                         "reco-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),
                         "reco-hadNHit", "Number of hits in hadronic cluster.  ");
  variables.emplace_back( HistAxis("Calorimetric Hadronic Energy (GeV)",
                                  kHadronicEnergyBinning, kNumuHadCalE),
                         "reco-hadcalE", "Sum of calibrated energy deposits in hadronic cluster.  " );
  variables.emplace_back( HistAxis("Slice Maximum Y (m)",  kXYBins, kMaxY),
                         "reco-maxy", "Maximum Y position of slice.  " );
  variables.emplace_back( HistAxis("Number of Tracks in Slice",  Binning::Simple(14, 1, 15), kNKalman),
                         "reco-nkal", "Number of tracks in slice.  " );
  variables.emplace_back( HistAxis("Number of Hits in Slice",  Binning::Simple(50, 0, 500), kNHit),
                         "reco-nhit", "Number of hits in slice.  " );
  variables.emplace_back( HistAxis("Muon Track cos(#theta_{Z})",  Binning::Simple(50, 0,1), kDirZ),
                         "reco-dirZ", "Z-direction of muon track.  " );
  variables.emplace_back( HistAxis("Track Start X Position (m)",  kXYBins, kTrkStartX),
                         "reco-trkStartX", "Track start x position.  " );
  variables.emplace_back( HistAxis("Track Start Y Position (m)",  kXYBins, kTrkStartY),
                         "reco-trkStartY", "Track start y position.  " );
  variables.emplace_back( HistAxis("Track Start Z Position (m)",  kZBins, kTrkStartZ),
                         "reco-trkStartZ", "Track start z position.  " );
  variables.emplace_back( HistAxis("Track End X Position (m)",  kXYBins, kTrkEndX),
                         "reco-trkEndX", "Track stop x position.  " );
  variables.emplace_back( HistAxis("Track End Y Position (m)",  kXYBins, kTrkEndY),
                         "reco-trkEndY", "Track stop y position.  " );
  variables.emplace_back( HistAxis("Track End Z Position (m)",  kZBins, kTrkEndZ),
                         "reco-trkEndZ", "Track stop z position.  " );
  variables.emplace_back( HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,3000), kSliceDuration),
                         "reco-sliceDuration", "Slice duration.  " );
  variables.emplace_back( HistAxis("Number of Hits in Primary Track",  Binning::Simple(50,0,500), kTrkNhits),
                         "reco-trkNhits", "Number of hits on primary track.  ");
  variables.emplace_back( HistAxis("Length of Primary Track (m)",  Binning::Simple(50,0,16), kTrkLength),
                         "reco-trkLength", "Primary track length.  " );
  variables.emplace_back( HistAxis("Scattering Log-likelihood",  Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
                         "reco-scatLL", "ReMId scattering log log-likelihood for primary track.  " );
  variables.emplace_back( HistAxis("dE/dx Log-likelihood",  Binning::Simple(50,-3,1), kReMIdDEDxLLH),
                         "reco-dedxLL", "ReMId dE/dx log log-likelihood for primary track.  " );
  variables.emplace_back( HistAxis("Non-hadronic Plane Fraction",
                                  Binning::Simple(50,0,1), kReMIdMeasFrac),
                         "reco-nonHadPlaneFrac", "ReMId Non-hadronic plane fraction.  " );
  variables.emplace_back( HistAxis("Muon ID",  kRemidBinning, kRemID),
                         "reco-remid", "ReMId kNN score.  " );
  ////////////////////////////////////////////////////////////////////////////////

  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;
}

#endif