caf_numu_fd_validation_MC_no_tau.C

Go to the documentation of this file.

#ifdef __CINT__
void caf_numu_fd_validation_MC_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/NumuVars.h"

#include "TH1.h"
#include "TFile.h"

using namespace ana;

const double pot = 36e20;

// And some binning definitions
const Binning kXYBins = Binning::Simple(12,-9,9);
const Binning kZBins = Binning::Simple(12,0,60);
const Binning kEnergyBinning = Binning::Simple(20,0,5);

// Cut on energy
const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"},
                [](const caf::StandardRecord* sr)
                {
                   float muE = sr->energy.mutrkE[sr->sel.remid.bestidx].E;
                   float hadE = sr->energy.numusimp.trkccE - muE;
                   return ( (sr->hdr.ismc && sr->energy.numusimp.trkccE < 5.0) 
                            || (!sr->hdr.ismc && muE+1.14f*hadE < 5.0) );
                });

// Cut for blinding
const Cut kBlindedRegion({"slc.calE", "slc.meantime", "hdr.run"},
     [](const caf::StandardRecord* sr)
     {
       return sr->slc.calE > 0.8 && sr->slc.calE < 3.0
         && ( (sr->slc.meantime > 207000 && sr->slc.meantime < 237000)
        || (sr->slc.meantime > 271000 && sr->slc.meantime < 301000 && sr->hdr.run <= 17500) );
     });

// Dumb oscillations
const Var kWOscDumb({"mc.nnu", "mc.nu.woscdumb"},
                  [](const caf::StandardRecord* sr)
                  {
                    if(sr->mc.nnu == 0) return 0.f;
                    return sr->mc.nu[0].woscdumb;
                  });

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

const Var kSlcMaxX({"slc.boxmax.fX"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.X()/100;
      });
const Var kSlcMaxY({"slc.boxmax.fY"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.Y()/100;
      });
const Var kSlcMaxZ({"slc.boxmax.fZ"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.Z()/100;
      });

void caf_numu_fd_validation_MC_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
{
  std::cout << "\nrun  : ---- Running CAF NuMU FD Validation - MC (no Tau).\n";
  std::cout << "run  :     input non-swap file name:    " << kInputNonSwapFileName << std::endl;
  std::cout << "run  :         input swap file name:    " << kInputSwapFileName    << std::endl;
  std::cout << "run  :             output file name:    " << kOutputFileName       << std::endl;

  SpectrumLoader loaderNonswap(kInputNonSwapFileName);
  SpectrumLoader loaderSwap(kInputSwapFileName);

  loaderNonswap.SetSpillCut(kStandardSpillCuts);
  loaderSwap.SetSpillCut(kStandardSpillCuts);

  struct Selection
  {
    std::string fname;
    Cut cut;
  };

  const int kNumCuts = 10;

  Selection vCut[kNumCuts] = 
    {
      {"nocut", (!kIsDytmanMEC) && kNoCut},
      // {"nocut", kNoCut},
      {"nocut_blinded", (!kIsDytmanMEC) && !kBlindedRegion},
      {"DQ", (!kIsDytmanMEC) && kNumuQuality},
      {"DQ_blinded", (!kIsDytmanMEC) && kNumuQuality && !kBlindedRegion},
      {"DQ_containment", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD},
      {"DQ_containment_blinded", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && !kBlindedRegion},
      {"DQ_containment_NCRej", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && kNumuNCRej},
      {"DQ_containment_NCRej_blinded", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && kNumuNCRej && !kBlindedRegion},
      {"all_numu_cuts", (!kIsDytmanMEC) && kNumuFD && kCCEshiftBelow5},
      {"all_numu_cuts_blinded", (!kIsDytmanMEC) && kNumuFD && kCCEshiftBelow5 && !kBlindedRegion}
    };

  struct Plot
  {
    std::string label;
    std::string fname;
    Binning bins;
    Var var;
  };

  const int kNumPlots = 34;

  Plot plots[kNumPlots] = {
    {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccEshift", kEnergyBinning, kCCEshift},
    {";#Deltat from t_{0} (#mus);Events / (12 #mus)", "tus", Binning::Simple(35,26.125,446.125), kSliceTime},
    {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccE", kEnergyBinning, kCCE},
    {";Calorimetric Energy (GeV);Events / (0.25 GeV)", "calE", kEnergyBinning, SIMPLEVAR(slc.calE)},
    {";Hadronic Energy (GeV);Events / (0.25 GeV)", "hadEshift", kEnergyBinning, kHadEshift},
    {";Calorimetric Hadronic Energy (GeV);Events / (0.25 GeV)", "hadcalE", kEnergyBinning, kNumuHadCalE},
    {";Slice Maximum x (m);Events / (1.5 m)", "maxX", kXYBins, kSlcMaxX},
    {";Slice Maximum y (m);Events / (1.5 m)", "maxY", kXYBins, kSlcMaxY},
    {";Slice Maximum z (m);Events / (5 m)", "maxZ", kZBins, kSlcMaxZ},
    {";Number of Tracks in Slice;Events / bin", "nkal", Binning::Simple(12, 1, 13), kNKalman},
    {";Total Number of Hits;Events / bin", "nhit", Binning::Simple(8, 0, 400), kNHit},
    {";Cosmic Rejection BDT Output;Events / bin", "bdt", Binning::Simple(10, 0.53, 0.73), kNumuContPID},
    {";Track Start X Position (m);Events / (1.5 m)", "trkStartX", kXYBins, kTrkStartX},
    {";Track Start Y Position (m);Events / (1.5 m)", "trkStartY", kXYBins, kTrkStartY},
    {";Track Start Z Position (m);Events / (5 m)", "trkStartZ", kZBins, kTrkStartZ},
    {";Track End X Position (m);Events / (1.5 m)", "trkEndX", kXYBins, kTrkEndX},
    {";Track End Y Position (m);Events / (1.5 m)", "trkEndY", kXYBins, kTrkEndY},
    {";Track End Z Position (m);Events / (5 m)", "trkEndZ", kZBins, kTrkEndZ},
    {";Number of Planes to Front;Events / bin", "nPlanesFront", Binning::Simple(15,0,900), kNplanesToFront},
    {";Number of Planes to Back;Events / bin", "nPlanesBack", Binning::Simple(15,0,900), kNplanesToBack},
    {";Number of Cells to Edge;Events / bin", "nCellsFromEdge", Binning::Simple(18,0,180), kNcellsFromEdge},
    {";Basic Track Forward Distance to Edge (m);Events / (5 m)", "cosFwdDist", kZBins, kCosFwdDist},
    {";Basic Track Backward Distance to Edge (m);Events / (5 m)", "cosBackDist", kZBins, kCosBackDist},
    {";Kalman Track Forward Distance to Edge (m);Events / (5 m)", "kalFwdDist", kZBins, kKalmanFwdDist},
    {";Kalman Track Backward Distance to Edge (m);Events / (5 m)", "kalBackDist", kZBins, kKalmanBackDist},
    {";Slice Duration (ns);Events / (200 ns)", "sliceDuration", Binning::Simple(15,0,3000), kSliceDuration},
    {";Number of Hits in Main Track;Events / bin", "trkNhits", Binning::Simple(20,0,500), kTrkNhits},
    {";Number of Hits not in Main Track;Events / bin", "hadNhit", Binning::Simple(8, 0, 160), kHadNHit},
    {";Length of Main Track (m);Events / (2 m)", "trkLength", Binning::Simple(15,0,30), kTrkLength},
    {";Scattering Log-Likelihood;Events / bin", "scatLL", Binning::Simple(11,-0.2,0.3), kReMIdScatLLH},
    {";dE/dx Log-Likelihood;Events / bin", "dedxLL", Binning::Simple(12,-0.1,0.5), kReMIdDEDxLLH},
    {";cos#theta_{beam};Events / bin", "angkal", Binning::Simple(12, 0.3, 1), kNumuCosRejAngleKal},
    {";cos#theta_{Z};Events / bin", "dirZ", Binning::Simple(12, 0.3, 1), kDirZ},
    {";Muon ID;Events / bin", "remid", Binning::Simple(14,0.75-1./160,1.+1./160), kRemID}
  };

  Spectrum* hNonswap[kNumPlots][kNumCuts];
  Spectrum* hSwap[kNumPlots][kNumCuts];

  for(int j = 0; j < kNumCuts; ++j){

    Selection s = vCut[j];

    for(int i = 0; i < kNumPlots; ++i){
  
      Plot p = plots[i];
          
      hNonswap[i][j] = new Spectrum(p.label, p.bins, loaderNonswap, p.var, s.cut, kNoShift, kWOscDumb);
      hSwap[i][j] = new Spectrum(p.label, p.bins, loaderSwap, p.var, s.cut, kNoShift, kWOscDumb);

    }
  }

  // GO!
  std::cout << "\nrun  : --- run loaders.\n";
  loaderNonswap.Go();
  loaderSwap.Go();
  std::cout << "\nrun  : --- done.\n";

  std::cout << "\nrun  : --- save output.\n";
  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");

  for(int j = 0; j < kNumCuts; ++j){
    for(int i = 0; i < kNumPlots; ++i){

      std::string fullName = "reco-"+plots[i].fname+"-"+vCut[j].fname;

      TH1* NonSwap = hNonswap[i][j]->ToTH1(pot);
      TH1* FluxSwap = hSwap[i][j]->ToTH1(pot);

      NonSwap->SetName(fullName.c_str());
      NonSwap->Add(FluxSwap);

      NonSwap->Write();

    } // End loop over i
  } // End loop over j

  TH1D* TotalPOT = new TH1D("meta-TotalPOT",";;Total POT", 1, 0, 1);
  TotalPOT->SetBinContent(1, pot);
  TotalPOT->Write();

  f->Write();
  f->Close();

  std::cout << "\nrun  : --- done.\n";

} // End of function

#endif