caf_nue_data_mc.C

Go to the documentation of this file.

#ifdef __CINT__
void caf_nue_data_mc(){}
#else

#include "CAFAna/Core/EventList.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Vars/Vars.h"
using namespace ana;

#include "OscLib/func/IOscCalculator.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TLegend.h"
#include "TPad.h"

#define ELASTICVAR(VAR) \
  Var({"vtx.nelastic", "vtx.elastic.*"}, \
    [](const caf::StandardRecord* sr) \
  { \
    if(sr->vtx.nelastic < 1) return -1000.; \
    return double(sr->vtx.elastic[0].vtx.VAR); \
  });

const Var kVtxX = ELASTICVAR(X());
const Var kVtxY = ELASTICVAR(Y());
const Var kVtxZ = ELASTICVAR(Z());

const Var kEPerHit({"slc.calE","slc.nhit"},
                   [](const caf::StandardRecord* sr)
                   {
                     return sr->slc.calE/sr->slc.nhit;
                   });

#define SHWLIDVAR(VAR) \
  Var({"shw.nshwlid", "shw.shwlid.*"}, \
    [](const caf::StandardRecord* sr) \
  { \
    if(sr->shw.nshwlid < 1) return -1000.; \
    return double(sr->shw.shwlid[0].VAR); \
  });

const Var kShwStartX = SHWLIDVAR(start.X());
const Var kShwStartY = SHWLIDVAR(start.Y());
const Var kShwStartZ = SHWLIDVAR(start.Z());

const Var kShwStopX = SHWLIDVAR(stop.X());
const Var kShwStopY = SHWLIDVAR(stop.Y());
const Var kShwStopZ = SHWLIDVAR(stop.Z());

const Var kShwCalE = SHWLIDVAR(calE);

const Var kShwNHit = SHWLIDVAR(nhit);
const Var kShwWidth = SHWLIDVAR(width);
const Var kShwLen = SHWLIDVAR(len);
const Var kShwGap = SHWLIDVAR(gap);

const Var kHadCalE({"shw.nshwlid", "shw.shwlid.*","slc.calE"},
                   [](const caf::StandardRecord* sr)
                   {
                     if(sr->shw.nshwlid < 1) return sr->slc.calE;
                     return (sr->slc.calE -
                             sr->shw.shwlid[0].calE);
                   });

#define ELECIDVAR(VAR) \
  Var({"sel.elecid.nshwlid", "sel.elecid.shwlid.*"}, \
    [](const caf::StandardRecord* sr) \
  { \
    if(sr->sel.elecid.nshwlid < 1) return -1000.; \
    return double(sr->sel.elecid.shwlid[0].VAR); \
  });

const Var kCosTheta = ELECIDVAR(costheta);
const Var kPi0Mass = ELECIDVAR(pi0mass);
const Var kShwEFrac = ELECIDVAR(shwEFrac);

const Var kMuLLL = ELECIDVAR(mulll);
const Var kMuLLT = ELECIDVAR(mullt);
const Var kELLL = ELECIDVAR(elll);
const Var kELLT = ELECIDVAR(ellt);
const Var kEGLLL = ELECIDVAR(eglll);
const Var kEGLLT = ELECIDVAR(egllt);
const Var kEPLLL = ELECIDVAR(eplll);
const Var kEPLLT = ELECIDVAR(epllt);
const Var kENLLL = ELECIDVAR(enlll);
const Var kENLLT = ELECIDVAR(enllt);
const Var kEMuLLL = ELECIDVAR(emulll);
const Var kEMuLLT = ELECIDVAR(emullt);
const Var kEPi0LLL = ELECIDVAR(epi0lll);
const Var kEPi0LLT = ELECIDVAR(epi0llt);
const Var kEPiLLL = ELECIDVAR(epilll);
const Var kEPiLLT = ELECIDVAR(epillt);


const Var kConfusion({"sel.lem.pid", "sel.elecid.ann"},
                     [](const caf::StandardRecord* sr)
                     {
                       if(sr->sel.lem.pid < .8 && sr->sel.elecid.ann < .95) return .5;
                       if(sr->sel.lem.pid < .8 && sr->sel.elecid.ann >= .95) return 1.5;
                       if(sr->sel.lem.pid >= .8 && sr->sel.elecid.ann < .95) return 2.5;
                       if(sr->sel.lem.pid >= .8 && sr->sel.elecid.ann >= .95) return 3.5;
                       
                     });

void caf_nue_data_mc()
{
  SpectrumLoader loaderData("defname: prod_caf_S15-05-22a_nd_numi_goodruns with stride 100");
  SpectrumLoader loaderDataStagger("defname: prod_caf_S15-05-22a_nd_numi_stagger_goodruns with stride 100");
  SpectrumLoader loaderMC("defname: prod_caf_S15-05-22_nd_genie_fhc_nonswap with stride 100");

  loaderData.SetSpillCut(kStandardSpillCuts);
  loaderDataStagger.SetSpillCut(kStandardSpillCuts);
  loaderMC.SetSpillCut(kStandardDQCuts);

  const Var kPIDExp = SIMPLEVAR(sel.lem.pidexp);
  const Var kMeany = SIMPLEVAR(sel.lem.meanyexp);
  const Var kQFrac = SIMPLEVAR(sel.lem.meanqfracexp);
  const Var kEDiff = SIMPLEVAR(sel.lem.energydiffexp);
  const Var kEnrich = SIMPLEVAR(sel.lem.enrichfracexp);
  const Var kNElastic = SIMPLEVAR(vtx.nelastic);
  const Var kNShwLID = SIMPLEVAR(shw.nshwlid);
  const Var kLongestProng = SIMPLEVAR(sel.rvp.longtr);

  //  const Var kPtp  = SIMPLEVAR(sel.nuecosrej.photptp);
  const Var kPtp({"sel.nuecosrej.photptp"},
                 [](const caf::StandardRecord* sr)
                 {
                   // This happened in the MC
                   if(isnan(sr->sel.nuecosrej.photptp) ||
                      isinf(sr->sel.nuecosrej.photptp)) return -1.f;
                   return sr->sel.nuecosrej.photptp;
                 });

  const Var kHitsPerPlane = SIMPLEVAR(sel.nuecosrej.hitsperplane);

  const Var kElecId = SIMPLEVAR(sel.elecid.ann);


  struct HistDef
  {
    std::string name;
    HistAxis axis;
  };

  const int kNumVars = 52;
  const HistDef defs[kNumVars] = {
    {"lem", {"LEM", Binning::Simple(100, 0, 1), kLEM}},
    {"lem_zoom", {"LEM", Binning::Simple(20, .8, 1), kLEM}},
    {"nhits", {"Number of hits", Binning::Simple(50, 0, 250), kNHit}},
    {"lem_pidexp", {"f_{sig}", Binning::Simple(100, 0, 1), kPIDExp}},
    {"lem_meany", {"<y>", Binning::Simple(80, 0, .8), kMeany}},
    {"lem_qfrac", {"<f_{Q}>", Binning::Simple(80, 0, .8), kQFrac}},
    {"lem_ediff", {"D", Binning::Simple(100, -.02, +.02), kEDiff}},
    {"lem_enrich", {"f_{enr}", Binning::Simple(100, 0, 1), kEnrich}},
    {"calE", {"E_{cal}", Binning::Simple(100, 0, 5), kCaloE}},
    {"nvtx", {"Number of vertices", Binning::Simple(3, 0, 3), kNElastic}},
    {"vtxx", {"Vertex X", Binning::Simple(50, -200, +200), kVtxX}},
    {"vtxy", {"Vertex Y", Binning::Simple(50, -200, +200), kVtxY}},
    {"vtxz", {"Vertex Z", Binning::Simple(40, 0, +800), kVtxZ}},

    {"nshw", {"Number of showers", Binning::Simple(10, 0, 10), kNShwLID}},
    {"shwstartx", {"Shower X Start", Binning::Simple(50, -200, +200), kShwStartX}},
    {"shwstarty", {"Shower Y Start", Binning::Simple(50, -200, +200), kShwStartY}},
    {"shwstartz", {"Shower Z Start", Binning::Simple(50, 0, +1000), kShwStartZ}},
    {"shwstopx", {"Shower X Stop", Binning::Simple(50, -200, +200), kShwStopX}},
    {"shwstopy", {"Shower Y Stop", Binning::Simple(50, -200, +200), kShwStopY}},
    {"shwstopz", {"Shower Z Stop", Binning::Simple(70, 0, +1400), kShwStopZ}},
    {"longprong", {"Longest prong (cm)", Binning::Simple(60, 0, 600), kLongestProng}},
    {"e_per_hit", {"Calorimetric energy / hit (GeV)", Binning::Simple(40, 0, .04), kEPerHit}},
    {"shw_nhit", {"Number of hits in primary shower", Binning::Simple(50, 0, 300), kShwNHit}},
    {"shw_len", {"Length of primary shower (cm)", Binning::Simple(50, 0, 500), kShwLen}},
    {"shw_width", {"Width of primary shower (cm)", Binning::Simple(20, 0, 20), kShwWidth}},
    {"shw_gap", {"Primary shower gap to vertex (cm)", Binning::Simple(25, 0, 200), kShwGap}},
    {"shw_efrac", {"Fraction of slice energy in shower", Binning::Simple(52, -.02, 1.02), kShwEFrac}},

    {"shw_e", {"Primary shower calorimetric energy (GeV)", Binning::Simple(30, 0, 3), kShwCalE}},
    {"had_e", {"(Slice - Shower) calorimetric energy (GeV)", Binning::Simple(30, 0, 3), kHadCalE}},

    {"ptp", {"p_{T}/p", Binning::Simple(52, -.02, 1.02), kPtp}},
    {"costheta", {"cos#theta", Binning::Simple(25, 0, +1), kCosTheta}},
    {"hits_per_plane", {"Hits per plane", Binning::Simple(10, 0, 10), kHitsPerPlane}},
    {"pi0mass", {"#pi^{0} mass (GeV)", Binning::Simple(50, 0, .5), kPi0Mass}},

    {"lid", {"LID", Binning::Simple(120, -0.1, 1.1), kElecId}},
    {"lid_zoom", {"LID", Binning::Simple(20, 0.95, 1.05), kElecId}},

    {"mulll", {"#mu LLL", Binning::Simple(50, -1.2, 3), kMuLLL}},
    {"elll", {"e LLL", Binning::Simple(50, -1.2, 3), kELLL}},
    {"emulll", {"e/#mu LLL", Binning::Simple(50, -1.2, 3), kEMuLLL}},
    {"eglll", {"e/#gamma LLL", Binning::Simple(50, -1.2, 1.2), kEGLLL}},
    {"eplll", {"e/p LLL", Binning::Simple(50, -1.2, 1.2), kEPLLL}},
    {"enlll", {"e/n LLL", Binning::Simple(50, -1.2, 1.2), kENLLL}},
    {"epilll", {"e/#pi LLL", Binning::Simple(50, -1.2, 1.2), kEPiLLL}},
    {"epi0lll", {"e/#pi^{0} LLL", Binning::Simple(50, -1.2, 1.2), kEPi0LLL}},

    {"mullt", {"#mu LLT", Binning::Simple(80, 0, 8), kMuLLT}},
    {"ellt", {"e LLT", Binning::Simple(80, 0, 8), kELLT}},
    {"emullt", {"e/#mu LLT", Binning::Simple(30, -1, +2), kEMuLLT}},
    {"egllt", {"e/#gamma LLT", Binning::Simple(30, -1, +2), kEGLLT}},
    {"epllt", {"e/p LLT", Binning::Simple(30, -1, +2), kEPLLT}},
    {"enllt", {"e/n LLT", Binning::Simple(30, -1, +2), kENLLT}},
    {"epillt", {"e/#pi LLT", Binning::Simple(30, -1, +2), kEPiLLT}},
    {"epi0llt", {"e/#pi^{0} LLT", Binning::Simple(30, -1, +2), kEPi0LLT}},
    {"confuse", {"Neither, LID, LEM, Both", Binning::Simple(4, 0, 4), kConfusion}},
  };

  const int kNumSels = 4;
  const Cut sels[kNumSels] = {kNoCut, kNueNDFirstAnaPresel, kNueNDFirstAnaPreselPlusLEM, kNueNDFirstAnaPreselPlusLID};
  const std::string selNames[kNumSels] = {"nocut", "presel", "lem", "lid"};

  Spectrum* spects[kNumSels][kNumVars];
  Spectrum* spectsStagger[kNumSels][kNumVars];
  IPrediction* preds[kNumSels][kNumVars];

  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const HistAxis& axis = defs[varIdx].axis;
      spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx]);
      spectsStagger[selIdx][varIdx] = new Spectrum(loaderDataStagger, axis, sels[selIdx]);
      preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
                                                     axis.label, axis.bins, axis.var,
                                                     sels[selIdx]);
    }
  }

  loaderMC.Go();
  loaderData.Go();
  loaderDataStagger.Go();

  TFile* fout = new TFile("fout.root", "RECREATE");
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* name = defs[varIdx].name.c_str();
      spects[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_%s", name)));
      spectsStagger[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_%s_stagger", name)));
      preds[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_%s", name)));
    }
  }
}

#endif