DataVtxDistributionLoad.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 <string>
#include <utility>

using namespace ana;

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

  // 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)";

  Binning kXYBinsFD     = Binning::Simple(320,  -800.,  800.);
  Binning kZBinsFD      = Binning::Simple(1200,    0., 6000.);

  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 HistAxis kAxisVtxX(labelVtxX, kXYBinsFD, kRecoVtxX);
  const HistAxis kAxisVtxY(labelVtxY, kXYBinsFD, kRecoVtxY);
  const HistAxis kAxisVtxZ(labelVtxZ, kZBinsFD,  kRecoVtxZ);

  // set up the cuts map
  std::map<std::string, Cut*> cuts;
  cuts["All"] = new Cut(kNusFD);

  std::map<std::string, Spectrum*> sXYs;
  std::map<std::string, Spectrum*> sXZs;
  std::map<std::string, Spectrum*> sYZs;

  for(const auto& strcut : cuts) {
    std::string label = strcut.first;

    sXYs[label] = new Spectrum(loaderFDdata, kAxisVtxX, kAxisVtxY, kInBeamSpill && *strcut.second);
    sXZs[label] = new Spectrum(loaderFDdata, kAxisVtxZ, kAxisVtxX, kInBeamSpill && *strcut.second);
    sYZs[label] = new Spectrum(loaderFDdata, kAxisVtxZ, kAxisVtxY, kInBeamSpill && *strcut.second);
  }

  // Fill the spectra!
  loaderFDdata.Go();

  // 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& strcut : cuts) {
    std::string label = strcut.first;

    dir = "sXY" + sep + label;
    saveDir = rootF->mkdir(dir.c_str());
    sXYs[label]->SaveTo(saveDir);

    dir = "sXZ" + sep + label;
    saveDir = rootF->mkdir(dir.c_str());
    sXZs[label]->SaveTo(saveDir);

    dir = "sYZ" + sep + label;
    saveDir = rootF->mkdir(dir.c_str());
    sYZs[label]->SaveTo(saveDir);
  }

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