AnaResultsLoad.C

Go to the documentation of this file.

// This macro creates spectra for OBSERVED
// and predicted event rates in analysis region
// This uses UNBLIND DATA. Proceed with extreme caution!

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/Ana01/Ana01Data.h"
#include "NuXAna/Vars/HistAxes.h"

#include <string>
#include <utility>

using namespace ana;

void AnaResultsLoad(std::string outfile)
{
  TH1::AddDirectory(0);

  // ND DATA + MC
  SpectrumLoader loaderNDMC(fnamenear_concat);
  SpectrumLoader loaderNDdata(fnameneardata_concat);
  loaderNDMC.SetSpillCut(kStandardSpillCuts);
  loaderNDdata.SetSpillCut(kStandardSpillCuts);
  
    // PERIOD 1, 2 + EPOCH 3B
  SpectrumLoader loaderFDMC_non(fFDMC_non);
  SpectrumLoader loaderFDMC_swp(fFDMC_swp);
  SpectrumLoader loaderFDMC_tau(fFDMC_tau);
  loaderFDMC_non.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_swp.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_tau.SetSpillCut(kStandardSpillCuts);

  // PERIOD 3C (and implied 3D)
  SpectrumLoader loaderFDMC_non_3c(fFDMC_non_3c);
  SpectrumLoader loaderFDMC_swp_3c(fFDMC_swp_3c);
  SpectrumLoader loaderFDMC_tau_3c(fFDMC_tau_3c);
  loaderFDMC_non_3c.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_swp_3c.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_tau_3c.SetSpillCut(kStandardSpillCuts);

  // FD DATA
  SpectrumLoader loaderFDdata(fnamefardata_unblind);
  loaderFDdata.SetSpillCut(kStandardSpillCuts);

  std::string labelVtxX    = "X Vertex Position (cm)";
  std::string labelVtxY    = "Y Vertex Position (cm)";
  std::string labelVtxZ    = "Z Vertex Position (cm)";
  std::string labelHits    = "Number of Hits";
  std::string labelCVN     = "CVN NCID";
  std::string labelNCP     = "NumucontPID";
  std::string labelPTP     = "PT/P";
  std::string labelEPerHit = "Average E/Hit (GeV)";
  std::string labelTop     = "Dist. to Det. Top (cm)";

  Binning kXYBinsFD     = Binning::Simple(64, -800., 800.);
  Binning kZBinsFD      = Binning::Simple(60,    0., 6000.);
  Binning kNearEdgeBins = Binning::Simple(64, 0., 1600.);
  Binning kFatPIDBins   = Binning::Simple(25, 0., 1.);
  Binning kHitBins      = Binning::Simple(50, 0., 200.);
  Binning kEDepBins     = Binning::Simple(25, 0., .1);

  const Var kRecoVtxX(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.; }
      return sr->vtx.elastic.vtx.x;
    });
  const Var kRecoVtxY(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.; }
      return sr->vtx.elastic.vtx.y;
    });
  const Var kRecoVtxZ(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return 9000.; }
      return sr->vtx.elastic.vtx.z;
    });

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

  std::map<std::string, HistAxis*> axes;
  axes["CalE"]    = new HistAxis(kNCAxis);
  axes["VtxX"]    = new HistAxis(labelVtxX,    kXYBinsFD,     kRecoVtxX);
  axes["VtxY"]    = new HistAxis(labelVtxY,    kXYBinsFD,     kRecoVtxY);
  axes["VtxZ"]    = new HistAxis(labelVtxZ,    kZBinsFD,      kRecoVtxZ);
  axes["CVN"]     = new HistAxis(labelCVN,     kFatPIDBins,   kCVNnc);
  axes["NHit"]    = new HistAxis(labelHits,    kHitBins,      kNHit);
  axes["CosBDT"]  = new HistAxis(labelNCP,     kFatPIDBins,   kNCP);
  axes["PTP"]     = new HistAxis(labelPTP,     kFatPIDBins,   kPartPTP);
  axes["EperHit"] = new HistAxis(labelEPerHit, kEDepBins,     kEPerHit);
  axes["DistTop"] = new HistAxis(labelTop,     kNearEdgeBins, kDistTop);

  std::map<std::string, IDecomp*> ncdecomps;
  std::map<std::string, IDecomp*> numudecomps;
  std::map<std::string, ModularExtrapSterile*> extrap123bs;
  std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
  std::map<std::string, PredictionSterile*> predst123bs;
  std::map<std::string, PredictionSterile*> predst3c3ds;
  std::map<std::string, PredictionCombinePeriods*> preds;
  std::map<std::string, Spectrum*> cosmics;
  std::map<std::string, Spectrum*> dataspecs;

  for(const auto& ax : axes) {
    std::string axlabel = ax.first;

    ncdecomps[axlabel] = new ProportionalDecomp(
      loaderNDMC, loaderNDdata,
      *ax.second, kNusND,
      kNoShift, kNoShift, kTuftsWeightCC
    );

    numudecomps[axlabel] = new ProportionalDecomp(
      loaderNDMC, loaderNDdata,
      kNCBinsNumuCCAxis, kNumuND,
      kNoShift, kNoShift, kTuftsWeightCC
    );

    extrap123bs[axlabel] = new ModularExtrapSterile(
      ModularExtrapSterile::NCDisappearance(
        loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
        *ncdecomps[axlabel], *numudecomps[axlabel],
        *ax.second, kNCBinsNumuCCAxis,
        kNusFD, kNusND, kNumuND,
        kNoShift, kTuftsWeightCC
      )
    );

    extrap3c3ds[axlabel] = new ModularExtrapSterile(
      ModularExtrapSterile::NCDisappearance(
        loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
        *ncdecomps[axlabel], *numudecomps[axlabel],
        *ax.second, kNCBinsNumuCCAxis,
        kNusFD, kNusND, kNumuND,
        kNoShift, kTuftsWeightCC
      )
    );

    predst123bs[axlabel] = new PredictionSterile(extrap123bs[axlabel]);
    predst3c3ds[axlabel] = new PredictionSterile(extrap3c3ds[axlabel]);

    cosmics[axlabel] = new Spectrum(loaderFDdata, *ax.second,
                                    kInTimingSideband && kNusFD,
                                    kNoShift, kTimingSidebandWeight);

    dataspecs[axlabel] = new Spectrum(loaderFDdata, *ax.second, kInBeamSpill && kNusFD);
  }

  // Fill the spectra!
  loaderNDMC.Go();
  loaderNDdata.Go();
  loaderFDMC_non.Go();
  loaderFDMC_swp.Go();
  loaderFDMC_tau.Go();
  loaderFDMC_non_3c.Go();
  loaderFDMC_swp_3c.Go();
  loaderFDMC_tau_3c.Go();
  loaderFDdata.Go();

  // Two scale factors for scaling spectra, 1 data yr and 3 data yr
  const double pot_p1p2e3b = kSecondAnaPeriod1POT+kSecondAnaPeriod2POT+kSecondAnaEpoch3bPOT;
  const double pot_e3ce3d  = kSecondAnaEpoch3cPOT+kSecondAnaEpoch3dPOT;

  for(const auto& ax : axes) {
    std::string axlabel = ax.first;

    preds[axlabel] = new PredictionCombinePeriods({
      std::make_pair(predst123bs[axlabel], pot_p1p2e3b),
      std::make_pair(predst3c3ds[axlabel], pot_e3ce3d)
    });
  }

  // Set up output filename
  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");

  // Save all of the objects
  TDirectory* tmp = gDirectory;
  TDirectory* saveDir = gDirectory;
  std::string dir;
  std::string sep = "__";

  for(const auto& ax : axes) {
    std::string axlabel = ax.first;

    dir = "ncdecomp" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    ncdecomps[axlabel]->SaveTo(saveDir);

    dir = "numudecomp" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    numudecomps[axlabel]->SaveTo(saveDir);

    dir = "extrap123b" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    extrap123bs[axlabel]->SaveTo(saveDir);

    dir = "extrap3c3d" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    extrap3c3ds[axlabel]->SaveTo(saveDir);

    dir = "predst123b" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    predst123bs[axlabel]->SaveTo(saveDir);

    dir = "predst3c3d" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    predst3c3ds[axlabel]->SaveTo(saveDir);

    dir = "pred" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    preds[axlabel]->SaveTo(saveDir);

    dir = "cosmic" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    cosmics[axlabel]->SaveTo(saveDir);

    dir = "dataspec" + sep + axlabel;
    saveDir = rootF->mkdir(dir.c_str());
    dataspecs[axlabel]->SaveTo(saveDir);
  }

  tmp->cd();
  rootF->Close(); // Close the file
}