validation_numu_nd.C

Go to the documentation of this file.

// Fills spectra for systematic error band

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Vars/HistAxes.h"
//#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TGaxis.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TList.h"
#include "TGaxis.h"
#include "TAttLine.h"
#include "TAttMarker.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include <cmath>
#include <iostream>
#include <vector>
#include <list>
#include <sstream>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>


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


std::string GetOutputFileName(std::string& name){
  int pos = 0;
  for (size_t i=0; i<name.size(); i++){
    if(name[i] == '/') pos = i;
  }
  std::string tmp = name.substr(pos+1, name.size());
  return tmp.substr(0, tmp.size()-5) + "_web.root";

}

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 sel, TangibleAxis tanAxis, SpectrumLoader& loaderMC, const Cut& bkg=kNoCut):
    fSel(sel), 
    fAxis(tanAxis),
    fHist(loaderMC, fAxis.fObj, fSel.fObj),
    fName(tanAxis.fShortName + "_" + sel.fShortName){};

  Selection fSel;
  TangibleAxis fAxis; 
  Spectrum fHist;
  std::string fName;

};


void validation_numu_nd()
{


  std::cout << "\n\n\n"
            << "This script appears to be quite old. I've hacked some stuff to remove qepid, but I doubt it will run out of the box."
            << "\n\n"
            << "If you choose to run it, best of luck to you..."
            << "\n\n\n";


  gStyle->SetMarkerStyle(kFullCircle);
  TGaxis::SetMaxDigits(3);
  //std::string fnameMC ("/nova/ana/users/tamsett/data/sam_definitions/hadd_rocco_decaf_S15-05-22_nd_genie_fhc_nonswap_genierw_numu_contain.root");
  std::string fnameMC ("/nova/ana/nu_e_ana/concat/prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_fdfluxv08_numu_contain.root");
  std::string path = "numu/Validation/";
  std::string fnameOutput = path + GetOutputFileName(fnameMC);
  
  SpectrumLoader loaderMC(fnameMC);

  TH1D* spillHist = new TH1D("spills", ";Detector;Spills", 3, 0, 3);
  SpillVar spillRun([](const caf::SRSpill* spill) {return spill->det;});

  loaderMC.SetSpillCut(kStandardSpillCuts);
  loaderMC.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);

  const Binning kXYBins = Binning::Simple(55,-2.20,2.20);
  const Binning kZBins = Binning::Simple(32,  0,16.00);
  const Binning kEnergyBinning = Binning::Simple(50,0,5);
  const Binning kHadronicEnergyBinning = Binning::Simple(50,0,5);
  
  std::vector<Selection> selections;

  selections.emplace_back(kNumuContainND && kNumuQuality, "NumuContainND",
    "Selected events pass numu containment and slice quality cuts.  ");
  selections.emplace_back(kNumuContainND && kNumuQuality && kNumuNCRej, "NumuND", 
    "Selected events pass numu containment, slice quality, and ReMId cuts.  ");

  std::vector<TangibleAxis> variables;

  variables.emplace_back(
    HistAxis("Reconstructed Neutrino Energy [GeV]", kEnergyBinning, kCCE), 
      "numuE", "CC energy estimator.  ");

  variables.emplace_back(
    HistAxis("Slice N_{Hit}",  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("Reconstructed Hadronic Energy [GeV]",  kHadronicEnergyBinning, kHadE), 
      "hadE", "Hadronic enery, i.e. numu energy estimate minus muon track energy.  " );

  variables.emplace_back(
    HistAxis("Muon Track Energy [GeV]",  Binning::Simple(50, 0, 5), kMuE),
      "numuTrackE", "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("Average 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}",  Binning::Simple(50, 0, 100), kHadNHit),
      "hadNHit", "Number of hits in hadronic cluster.  ");

  variables.emplace_back(
    HistAxis("Off-track Calorimetric Hadronic Energy [GeV]", kHadronicEnergyBinning, kNumuHadCalE),
      "hadcalE", "Sum of calibrated energy deposits in hadronic cluster.  " );

  variables.emplace_back(
    HistAxis("On-track Hadronic Energy [GeV]",  Binning::Simple(50, 0, 0.5), kNumuHadTrkE),
      "hadTrkE", "Hadronic Energy in active region, on the track.  " );

  variables.emplace_back(
    HistAxis("Visible On-track Hadronic Energy [GeV]",  kHadronicEnergyBinning, kNumuHadVisE),
      "hadVisE", "Hadronic Energy in active region, on the track.  " );
  
  variables.emplace_back(
    HistAxis("Slice Maximum Y [m]",  kXYBins, kSlcMaxY),
      "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("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.  " );
  
  variables.emplace_back(
    HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,500), 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("Muon Track Length [m]",  Binning::Simple(50,0,16), kTrkLength),
      "trkLength", "Primary track length.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: Scattering Log-likelihood",  Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
      "scatLL", "ReMId scattering log log-likelihood for primary track.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: dE/dx Log-likelihood",  Binning::Simple(50,-3,1), kReMIdDEDxLLH),
      "dedxLL", "ReMId dE/dx log log-likelihood for primary track.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: Non-hadronic Plane Fraction", Binning::Simple(50,0,1), kReMIdMeasFrac),
      "nonHadPlaneFrac", "ReMId Non-hadronic plane fraction.  " );
  
  variables.emplace_back(
    HistAxis("ReMId",  kRemidBinning, kRemID),
      "remid", "ReMId kNN score.  " );

  variables.emplace_back(
    HistAxis("Number of Uncalibrated Hits in Slice",  Binning::Simple(50, 0, 50), kSlcUnCalibNHit),
      "slcUnCalibNHit", "Uncalibrated slice nhit.  " );

  variables.emplace_back(
    HistAxis("Slice Mean Time [#mus]",  Binning::Simple(60, -50, 550), kSlcMeanTime),
      "slcMeanTime", "Slice mean time TNS in microseconds.  " );

  variables.emplace_back(
    HistAxis("Slice Start Time [#mus]",  Binning::Simple(60, -50, 550), kSlcStartTime),
      "slcStartTime", "Slice start time TNS in microseconds.  " );

  variables.emplace_back(
    HistAxis("Slice End Time [#mus]",  Binning::Simple(60, -50, 550), kSlcEndTime),
      "slcEndTime", "Slice end time TNS in microseconds.  " );

  variables.emplace_back(
    HistAxis("Slice Minimum X [m]",  kXYBins, kSlcMinX),
      "slcMinX", "Slice Min X [m].  " );

  variables.emplace_back(
    HistAxis("Slice Minimum Y [m]",  kXYBins, kSlcMinY),
      "slcMinY", "Slice Min Y [m].  " );

  variables.emplace_back(
    HistAxis("Slice Minimum Z [m]",  kZBins, kSlcMinZ),
      "slcMinZ", "Slice Min Z [m].  " );

  variables.emplace_back(
    HistAxis("Slice Maximum X [m]",  kXYBins, kSlcMaxX),
      "slcMaxX", "Slice Max X [m].  " );

  variables.emplace_back(
    HistAxis("Slice Maximum Y [m]",  kXYBins, kSlcMaxY),
      "slcMaxY", "Slice Max Y [m].  " );

  variables.emplace_back(
    HistAxis("Slice Maximum Z [m]",  kZBins, kSlcMaxZ),
      "slcMaxZ", "Slice Max Z [m].  " );

  variables.emplace_back(
    HistAxis("Slice Extent X [m]",  Binning::Simple(55, 0, 4.4), kSlcExtentX),
      "slcExtentX", "Slice Extent X [m].  " );

  variables.emplace_back(
    HistAxis("Slice Extent Y [m]",  Binning::Simple(55, 0, 4.4), kSlcExtentY),
      "slcExtentY", "Slice Extent Y [m].  " );

  variables.emplace_back(
    HistAxis("Slice Extent Z [m]",  Binning::Simple(18, 0, 16), kSlcExtentZ),
      "slcExtentZ", "Slice Extent Z [m].  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{X}",  Binning::Simple(100, -1, 1), kDirX),
      "dirX", "X-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{Y}",  Binning::Simple(100, -1, 1), kDirY),
      "dirY", "Y-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{Z}",  Binning::Simple(50, 0, 1), kDirZ),
      "dirZ", "Z-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{NuMI}",  Binning::Simple(50, 0, 1), kCosNumi),
      "cosNumi", "Beam direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Number of Missing Planes in Kalman Track",  Binning::Simple(50, 0, 50), 
                                                                                                                  kTrkNPlaneGap),
      "trkNPlaneGap", "Track N Plane Gap.  " );

  variables.emplace_back(
    HistAxis("Visible Slice Energy Per Hit [GeV]",  Binning::Simple(40, 0, 0.04),
                                                          kSlcCalEPerNHit),
      "slcCalEPerNHit", "Slice Energy Per Slice NHit.  " );



  std::vector<UsefulHist> hists;
  hists.reserve(selections.size() * variables.size());
  
  for(const auto& sel:selections){
    for(const auto& variable:variables){
      hists.emplace_back(sel, variable, loaderMC);
    }
  }

  loaderMC.Go();
  double pot = hists[0].fHist.POT();

  TFile  inputfile(fnameMC.c_str(),     "READ");
  TFile outputfile(fnameOutput.c_str(), "RECREATE");

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

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

  inputfile.Close();
  outputfile.Close();

}