FDDataMCSystBandLoad.C

Go to the documentation of this file.

// Macro to plot FD Data/MC spectra with systematic error band

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "NuXAna/Systs/NusSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSSystFunctions.h"
#include "NuXAna/macros/Ana01/Ana01Data.h"

#include "TFile.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "NuXAna/Vars/HistAxes.h"

#include <iostream>
#include <string>
#include <vector>

using namespace ana;

void PlotSysts(IPrediction* pred,
               const std::vector<const ISyst*>& systs,
               osc::IOscCalc* calc,
               TDirectory* out, FILE* text, strings strs,
               Spectrum* sData);

void PlotSysts(const Spectrum& nom,
               const std::vector<Spectrum>& upShifts,
               const std::vector<Spectrum>& downShifts,
               TDirectory* out, FILE* text, strings strs,
               Spectrum* sData);

void PlotSysts(IPrediction* pred,
               const ISyst* syst,
               osc::IOscCalc* calc,
               TDirectory* out, strings strs);

void PlotSysts(Spectrum& nom, Spectrum& sp1, Spectrum& sm1,
               TDirectory* out, strings strs, bool split = true);

TH1* GetHist(strings strs, int sigma);

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

  Loaders loaders123b;
  loaders123b.SetLoaderPath(fnamenear_concat, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders123b.SetLoaderFiles(fFDMC_non,        caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders123b.SetLoaderFiles(fFDMC_swp,        caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders123b.SetLoaderFiles(fFDMC_tau,        caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders123b.SetLoaderPath(fnameneardata_concat, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  Loaders loaders3c3d;
  loaders3c3d.SetLoaderPath(fnamenear_concat, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders3c3d.SetLoaderPath(fFDMC_non_3c,     caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders3c3d.SetLoaderPath(fFDMC_swp_3c,     caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders3c3d.SetLoaderPath(fFDMC_tau_3c,     caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders3c3d.SetLoaderPath(fnameneardata_concat, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  loaders123b.SetSpillCut(kStandardSpillCuts);
  loaders3c3d.SetSpillCut(kStandardSpillCuts);

  SpectrumLoader loaderFDdata(fnamefardata_unblind);
  loaderFDdata.SetSpillCut(kStandardSpillCuts);

  // Create vector of systematics
  std::vector<const ISyst*> systs = getAllNusSysts();

  SterileGenerator gen(kNCAxis, kNCBinsNumuCCAxis,
                       kNusFD, kNusND, kNumuND, kNoShift, kTuftsWeightCC);
  osc::OscCalcSterile* calc = DefaultSterileCalc(3);
  PredictionInterp pred123b(systs, calc, gen, loaders123b);
  PredictionInterp pred3c3d(systs, calc, gen, loaders3c3d);

  Spectrum* sCosmic = new Spectrum(loaderFDdata, kNCAxis,
                                   kInTimingSideband && kNusFD,
                                   kNoShift, kTimingSidebandWeight);

  Spectrum* sCosOOT = new Spectrum(loaderFDdata, kNCAxis,
                                   kInTimingSideband && kNusFD);

  Spectrum* sData = new Spectrum(loaderFDdata, kNCAxis, kInBeamSpill && kNusFD);

  loaders123b.Go();
  loaders3c3d.Go();
  loaderFDdata.Go();

  const double pot_p1p2e3b =
    kSecondAnaPeriod1POT
    + kSecondAnaPeriod2POT
    + kSecondAnaEpoch3bPOT;
  const double pot_e3ce3d  = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
  
  PredictionCombinePeriods pred({std::make_pair(&pred123b, pot_p1p2e3b),
                                 std::make_pair(&pred3c3d, pot_e3ce3d)});


  TFile* rootF = new TFile("Ana01_FDDataMCSystBandLoad.root", "RECREATE");

  TDirectory* tmp = gDirectory;
  TDirectory* saveDir = gDirectory;

  // Save nominal spectra
  saveDir = rootF->mkdir("sCosmic");
  sCosmic->SaveTo(saveDir);
  saveDir = rootF->mkdir("sCosOOT");
  sCosOOT->SaveTo(saveDir);
  saveDir = rootF->mkdir("sData");
  sData->SaveTo(saveDir);
  saveDir = rootF->mkdir("pred123b");
  pred123b.SaveTo(saveDir);
  saveDir = rootF->mkdir("pred3c3d");
  pred3c3d.SaveTo(saveDir);
  saveDir = rootF->mkdir("pred");
  pred.SaveTo(saveDir);

  rootF->Close();
}