caf_numu_nuenergy_vs_xy.C

Go to the documentation of this file.

// Uses truth information. Only run on MC

/*                                                                                                                                           
*/
#ifdef __CINT__
void caf_numu_nuenergy_vs_xy(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 "TH3.h"
#include "TProfile2D.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:

  UsefulHist(Selection selA, TangibleAxis tanXAxis, TangibleAxis tanYAxis, TangibleAxis tanZAxis, 
             SpectrumLoader& loader, const Cut& cutMEC = !kIsDytmanMEC, const Cut& bkg=kNoCut):
    fSelA(selA), 
    fXAxis(tanXAxis),
    fYAxis(tanYAxis),
    fZAxis(tanZAxis),
    fHist(loader, fXAxis.fObj, fYAxis.fObj, fZAxis.fObj, fSelA.fObj && cutMEC),
    fName("reco-" + fXAxis.fShortName + "_vs_" + fYAxis.fShortName + "_vs_" + fZAxis.fShortName + "-" + selA.fShortName ){};
  
  Selection fSelA;
  TangibleAxis fXAxis; 
  TangibleAxis fYAxis; 
  TangibleAxis fZAxis; 
  Spectrum fHist;
  std::string fName;

};


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

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

 
  // 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);
                     });
  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
 

  std::vector<TangibleAxis> variables;
  variables.emplace_back(kNumuCCAxis, "numuE",
                         "CC energy estimator (GeV)  ");


  TangibleAxis tanTrkStartX (HistAxis("Track Start X Position (m)", kXYBins, kTrkStartX), "trkStartX", "Track start x position. ");
  TangibleAxis tanTrkStartY (HistAxis("Track Start Y Position (m)", kXYBins, kTrkStartY), "trkStartY", "Track start y 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, tanTrkStartX, tanTrkStartY, 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 outputfile(kOutputFileName.c_str(), "RECREATE");

  for(const auto& hist:hists){
    TH3 *temp = hist.fHist.ToTH3(pot);
    std::string tempName = hist.fName;
    TProfile2D* tempProj = (TProfile2D*)temp -> Project3DProfile("yx");
    //temp->Write(tempName.c_str());
    tempProj->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();
  std::cout << "\nrun  : --- done.\n";



  // dictionary to translate histogram names 
  /*

  */

  return;
}

#endif