DataMCLoad.C

Go to the documentation of this file.

// A macro for ND data/MC comparisons
// Aaron Markowitz

#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/Vars.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/Vars/HistAxes.h"

using namespace ana;

std::pair<Spectrum*, CheatDecomp*> make_pair(
    SpectrumLoaderBase& loader_data,
    SpectrumLoaderBase& loader_mc,
    HistAxis* axis,
    Cut* cut,
    const SystShifts& shift,
    const Var& wei
);

std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
    SpectrumLoaderBase& loader_data,
    SpectrumLoaderBase& loader_mc,
    std::map<std::string, HistAxis* > axes,
    std::map<std::string, Cut* > cuts,
    const SystShifts& shift,
    const Var& wei
);

std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
    Loaders& loader,
    std::map<std::string, HistAxis*> axes,
    std::map<std::string, Cut*> cuts,
    const SystShifts& shift,
    const Var& wei
);

void DataMCLoad()
{
  TH1::AddDirectory(0);

  // Strings to file paths
  //const std::string fneardata = "defname: prod_caf_R16-03-03-prod2reco.d_nd_numi_fhc_full_v1_goodruns with limit 10";
  //const std::string fneardatafull = "prod_caf_R16-03-03-prod2reco.d_nd_numi_fhc_full_v1_goodruns";
  const std::string fneardatafull = "defname: prod_caf_R16-11-12-feature_caf_size.a_nd_numi_fhc_full_v1_goodruns with limit 1";
  //const std::string fnamenearcaftest = "defname: prod_caf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_v1 with limit 10";
  //const std::string fnearcaffull = "prod_caf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_v1";
  const std::string fnearcaffull = "defname: prod_caf_R16-11-12-feature_caf_size.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_v1 with limit 1";

  // Set up a Loaders object
  // For full SA data and MC filesets, use 'fneardatafull' and 'fnamenearcaf' respectively
  Loaders loaders;
  loaders.SetLoaderPath(fneardatafull, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnearcaffull,  caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap);

  // Axis definitions
  std::string labelHpP     = "Hits per Plane";
  std::string labelGap     = "Primary Shower Gap to Vertex (cm)";
  std::string labelVtxX    = "X Vertex Position (cm)";
  std::string labelVtxY    = "Y Vertex Position (cm)";
  std::string labelVtxZ    = "Z Vertex Position (cm)";
  std::string labelCont    = "Minimum Distance to Any Detector Face (cm)";
  std::string labelHits    = "Number of Hits";
  std::string labelCVN     = "CVN NC Classifier";
  std::string labelLID     = "LID";
  std::string labelRemID   = "RemID";
  std::string labelNCP     = "NumucontPID";
  std::string labelPTP     = "p_{T}/p";
  std::string labelIsMuon  = "Is Muon";
  std::string labelEPerHit = "Average E/Hit (GeV)";
  std::string labelTop     = "Distance to Top of Detector (cm)";

  Binning kNCAllEBins   = Binning::Simple(24, 0, 6);
  Binning kGapBins      = Binning::Simple(25, 0., 200.);
  Binning kXYBinsND     = Binning::Simple(100, -250., 250.);
  Binning kZBinsND      = Binning::Simple(410, -50., 2000.);
  Binning kNearEdgeBins = Binning::Simple(160, 0., 800.);
  Binning kFinePIDBins  = Binning::Simple(100, 0., 1.);
  Binning kHitBins      = Binning::Simple(800, 0., 800.);
  Binning kEDepBins     = Binning::Simple(200, 0., 0.2);
  Binning kHpPBins      = Binning::Simple(10, 0., 10.);

  const Var kGap(
               [](const caf::SRProxy* sr){
                  if(sr->vtx.elastic.fuzzyk.nshwlid > 0) {
                    return sr->vtx.elastic.fuzzyk.gap;
                    }
                  float output = 999.;
                  return output;
                 });
  const Var kHpP = SIMPLEVAR(sel.nuecosrej.hitsperplane);

  const Var kRecoVtxX(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.f; }
      return sr->vtx.elastic.vtx.X();
    });
  const Var kRecoVtxY(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.f; }
      return sr->vtx.elastic.vtx.Y();
    });
  const Var kRecoVtxZ(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.f; }
      return sr->vtx.elastic.vtx.Z();
    });

  const Var kMinDistP(
    [](const caf::SRProxy* sr)
    {
      const caf::SRNueCosRej& cr = sr->sel.nuecosrej;
      float mindist = 10000.;
      mindist = std::min(mindist, std::min(cr.starteast,   cr.stopeast));
      mindist = std::min(mindist, std::min(cr.startwest,   cr.stopwest));
      mindist = std::min(mindist, std::min(cr.starttop,    cr.stoptop));
      mindist = std::min(mindist, std::min(cr.startbottom, cr.stopbottom));
      mindist = std::min(mindist, std::min(cr.startfront,  cr.stopfront));
      mindist = std::min(mindist, std::min(cr.startback,   cr.stopback));
      return (mindist > 9999. ? (float)-1. : mindist);
    });

  const Var kNCP = SIMPLEVAR(sel.cosrej.numucontpid);
  const Var kPartPTP(
    [](const caf::SRProxy* sr)
    {
      if(std::isnan(sr->sel.nuecosrej.partptp)) { return -5.f; }
      return sr->sel.nuecosrej.partptp;
    });
  const Var kIsMuon([](const caf::SRProxy* sr)
                    {
                      if(sr->vtx.elastic.fuzzyk.nshwlid < 1) { return -1; }
                      return sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon;
                    });
  const Var kEPerHit(
    [](const caf::SRProxy* sr)
    {
      double nhit = (double)sr->slc.nhit;
      if(nhit <= 0.) { return 0.; }
      return sr->slc.calE/nhit;
    });
  const Var kDistTop(
    [](const caf::SRProxy* sr)
    {
      return std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop);
    });

  const HistAxis kNCAllEAxis("Calorimetric Energy (GeV)", kNCAllEBins, kCaloE);
  const HistAxis kAxisHpP( labelHpP,     kHpPBins,      kHpP);
  const HistAxis kAxisGap( labelGap,     kGapBins,      kGap);
  const HistAxis kAxisVtxX(labelVtxX,    kXYBinsND,     kRecoVtxX);
  const HistAxis kAxisVtxY(labelVtxY,    kXYBinsND,     kRecoVtxY);
  const HistAxis kAxisVtxZ(labelVtxZ,    kZBinsND,      kRecoVtxZ);
  const HistAxis kAxisCont(labelCont,    kNearEdgeBins, kMinDistP);
  const HistAxis kAxisCVN( labelCVN,     kFinePIDBins,  kCVNnc);
  const HistAxis kAxisLID( labelLID,     kFinePIDBins,  kElecID);
  const HistAxis kAxisRem( labelRemID,   kRemidBinning, kRemID);
  const HistAxis kAxisNHit(labelHits,    kHitBins,      kNHit);
  const HistAxis kAxisNCP( labelNCP,     kFinePIDBins,  kNCP);
  const HistAxis kAxisPTP( labelPTP,     kFinePIDBins,  kPartPTP);
  const HistAxis kAxisIsMu(labelIsMuon,  kBooleanBins,  kIsMuon);
  const HistAxis kAxisEpH( labelEPerHit, kEDepBins,     kEPerHit);
  const HistAxis kAxisTop( labelTop,     kNearEdgeBins, kDistTop);

  // Two intermediate cuts, created from cuts in NusCuts.h
  const Cut kEQFid = kNusEventQuality && kNusNDFiducial;
  const Cut kPreselNCSel = kNusNDPresel && kNusNCSel;

  const Cut kNM1Fid    = kNusEventQuality && kNusNDContain  && kNusCosRejMod && kNusNCSel;
  const Cut kNM1Cont   = kNusEventQuality && kNusNDFiducial && kNusCosRejMod && kNusNCSel;
  const Cut kNM1NCSel  = kNusNDPresel && kNusCosRejMod;
  const Cut kNM1CosRej = kNusNDPresel && kNusNCSel;

  const Cut kFidMinusX(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) return false;

      const TVector3 vtx = sr->vtx.elastic.vtx;
      if(vtx.Y() < -100.0) return false;
      if(vtx.Y() >  100.0) return false;
      if(vtx.Z() <  200.0) return false;
      if(vtx.Z() > 1000.0) return false;
      return true;
    });
  const Cut kFidMinusY(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) return false;

      const TVector3 vtx = sr->vtx.elastic.vtx;
      if(vtx.X() < -100.0) return false;
      if(vtx.X() >  100.0) return false;
      if(vtx.Z() <  200.0) return false;
      if(vtx.Z() > 1000.0) return false;
      return true;
    });
  const Cut kFidMinusZ(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid ) return false;

      const TVector3 vtx = sr->vtx.elastic.vtx;
      if(vtx.X() < -100.0) return false;
      if(vtx.X() >  100.0) return false;
      if(vtx.Y() < -100.0) return false;
      if(vtx.Y() >  100.0) return false;
      return true;
    });

  const Cut kNCSelMinusHit(
    [](const caf::SRProxy* sr)
    {
      if(sr->sel.cvn.ncid < 0.2) return false;
      return true;
    });
  const Cut kNCSelMinusCVN(
    [](const caf::SRProxy* sr)
    {
      if(sr->slc.nhit >= 200) return false;
      if(sr->slc.nhit <  20) return false;
      return true;
    });

  const Cut kCosRejMinusPTP(
    [](const caf::SRProxy* sr)
    {
      if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
      if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon == 1) return false;
      double nhit = (double)sr->slc.nhit;
      if(nhit <= 0.) return false;
      if(sr->slc.calE/nhit <= 0.009) return false;
      return true;
    });
  const Cut kCosRejMinusIsM(
    [](const caf::SRProxy* sr)
    {
      double partptp = sr->sel.nuecosrej.partptp;
      if(partptp >= 0.8) return false;
      double nhit = (double)sr->slc.nhit;
      if(nhit <= 0.) return false;
      if(sr->slc.calE/nhit <= 0.009) return false;
      return true;
    });
  const Cut kCosRejMinusEpH(
    [](const caf::SRProxy* sr)
    {
      double partptp = sr->sel.nuecosrej.partptp;
      if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
      if(partptp >= 0.8) return false;
      if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon == 1) return false;
      return true;
    });

  // set up the axes map
  std::map<std::string, HistAxis*> axes;
  axes["ECal"]      = new HistAxis(kNCAxis);
  axes["ECalLong"]  = new HistAxis(kNCAllEAxis);
  axes["HitpPlane"] = new HistAxis(kAxisHpP);
  axes["Gap"]       = new HistAxis(kAxisGap);
  axes["VtxX"]      = new HistAxis(kAxisVtxX);
  axes["VtxY"]      = new HistAxis(kAxisVtxY);
  axes["VtxZ"]      = new HistAxis(kAxisVtxZ);
  axes["Cont"]      = new HistAxis(kAxisCont);
  axes["CVN"]       = new HistAxis(kAxisCVN);
  axes["RemID"]     = new HistAxis(kAxisRem);
  axes["LID"]       = new HistAxis(kAxisLID);
  axes["NHit"]      = new HistAxis(kAxisNHit);
  axes["EpH"]       = new HistAxis(kAxisEpH);
  axes["PTP"]       = new HistAxis(kAxisPTP);
  axes["IsMu"]      = new HistAxis(kAxisIsMu);
  axes["DistTop"]   = new HistAxis(kAxisTop);

  // set up the cuts map
  std::map<std::string, Cut*> cuts;
  cuts["DQ"]          = new Cut(kNoCut);
  cuts["EQ"]          = new Cut(kNusEventQuality);
  cuts["EQFid"]       = new Cut(kEQFid);
  cuts["Presel"]      = new Cut(kNusNDPresel);
  cuts["PreselNCSel"] = new Cut(kPreselNCSel);
  cuts["NusND"]       = new Cut(kNusND);
  cuts["AllE"]        = new Cut(kNusNDAllE);

  // make a map containing data spectra and MC cheat decomps for each axis/cut pair
  std::map<std::string, std::pair<ana::Spectrum*, ana::CheatDecomp*>> loaded_specs = MakeMap(
    loaders, axes, cuts, kNoShift, kTuftsWeightCC
  );

  std::map<std::string, Cut*> nm1s;
  nm1s["VtxX"] = new Cut(kNM1Fid && kFidMinusX);
  nm1s["VtxY"] = new Cut(kNM1Fid && kFidMinusY);
  nm1s["VtxZ"] = new Cut(kNM1Fid && kFidMinusZ);
  nm1s["Cont"] = new Cut(kNM1Cont);
  nm1s["CVN"]  = new Cut(kNM1NCSel && kNCSelMinusCVN);
  nm1s["NHit"] = new Cut(kNM1NCSel && kNCSelMinusHit);
  nm1s["PTP"]  = new Cut(kNM1CosRej && kCosRejMinusPTP);
  nm1s["EpH"]  = new Cut(kNM1CosRej && kCosRejMinusEpH);
  nm1s["IsMu"] = new Cut(kNM1CosRej && kCosRejMinusIsM);
  for(const auto& nm1 : nm1s) {
    std::string strcut = nm1.first;
    std::string strNM1 = "NM1";

    loaded_specs[strcut + strNM1] = make_pair(
      loaders.GetLoader(caf::kNEARDET, Loaders::kData, kBeam, Loaders::kNonSwap),
      loaders.GetLoader(caf::kNEARDET, Loaders::kMC, kBeam, Loaders::kNonSwap),
      axes[strcut], nm1.second, kNoShift, kTuftsWeightCC
    );
  }

  // Go loaders
  loaders.SetSpillCut(kStandardSpillCuts);
  loaders.Go();

  // Open output file
  TFile* rootOut = new TFile("/nova/ana/steriles/Ana01/DataMC.root", "RECREATE");

  for(const auto& loaded_spec : loaded_specs) {
    loaded_spec.second.first->SaveTo(rootOut, (loaded_spec.first + "Data").c_str());
    loaded_spec.second.second->SaveTo(rootOut, (loaded_spec.first + "MC").c_str());
  }

  rootOut->Close();
}

//..............................................................................
std::pair<Spectrum*, CheatDecomp*> make_pair(SpectrumLoaderBase& loader_data,
                                             SpectrumLoaderBase& loader_mc,
                                             HistAxis* axis,
                                             Cut* cut,
                                             const SystShifts& shift,
                                             const Var& wei)
{
  std::pair<Spectrum*, CheatDecomp*> output;
  output.first = new Spectrum(
    loader_data,
    *axis, *cut,
    shift, wei
  );
  output.second = new CheatDecomp(
    loader_mc,
    *axis, *cut,
    shift, wei
  );

  return output;
}

//..............................................................................
std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
    SpectrumLoaderBase& loader_data,
    SpectrumLoaderBase& loader_mc,
    std::map<std::string, HistAxis* > axes,
    std::map<std::string, Cut* > cuts,
    const SystShifts& shift,
    const Var& wei
)
{
  std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > output;
  for(const auto& cut : cuts) {
    for(const auto& axis : axes) {
      output[axis.first + cut.first] = make_pair(
        loader_data, loader_mc,
        axis.second, cut.second, shift, wei
      );
    }
  }

  return output;
};

//..............................................................................
std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
    Loaders& loader,
    std::map<std::string, HistAxis*> axes,
    std::map<std::string, Cut*> cuts,
    const SystShifts& shift,
    const Var& wei
)
{
  return MakeMap(
    loader.GetLoader(caf::kNEARDET, Loaders::kData, kBeam, Loaders::kNonSwap),
    loader.GetLoader(caf::kNEARDET, Loaders::kMC, kBeam, Loaders::kNonSwap),
    axes, cuts, shift, wei
  );
};