FitSystEffectsLoad.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Loaders.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "NuXAna/Systs/NusSysts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/Ana01/Ana01Data.h"

#include "OscLib/OscCalcSterile.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "NuXAna/Vars/HistAxes.h"

#include <string>
#include <utility>

using namespace ana;

void FitSystEffectsLoad(std::string outfile)
{
  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);

  SterileGenerator gen(kNCAxis, kNCBinsNumuCCAxis,
                       kNusFD, kNusND, kNumuND, kNoShift, kTuftsWeightCC);

  std::map<std::string, std::vector<const ISyst*> > systmap;
  systmap["Beam"]     = std::vector<const ISyst*>();
  systmap["Birks"]    = std::vector<const ISyst*>();
  systmap["Calib"]    = std::vector<const ISyst*>();
  systmap["DataMC"]   = std::vector<const ISyst*>();
  systmap["GENIE"]    = std::vector<const ISyst*>();
  systmap["NDCont"]   = std::vector<const ISyst*>();
  systmap["NDRock"]   = std::vector<const ISyst*>();
  systmap["Norm"]     = std::vector<const ISyst*>();
  systmap["OscParam"] = std::vector<const ISyst*>();
  systmap["Stats"]    = std::vector<const ISyst*>();

  systmap["Beam"]    .push_back(&kNusBeamSysts);
  systmap["Birks"]   .push_back(&kNusBirksSyst);
  systmap["Calib"]   .push_back(&kNusCalibFlatSyst);
  systmap["Calib"]   .push_back(&kNusCalibRelSyst);
  systmap["Calib"]   .push_back(&kNusCalibSlopeXSyst);
  systmap["Calib"]   .push_back(&kNusCalibSlopeYSyst);
  systmap["DataMC"]  .push_back(&kNusNueCCSyst);
  systmap["DataMC"]  .push_back(&kNusNumuCCSyst);
  systmap["GENIE"]   .push_back(&kNusGENIESmallSysts);
  systmap["NDCont"]  .push_back(&kNusNDContSyst);
  systmap["NDRock"]  .push_back(&kNusNDRockSyst);
  systmap["Norm"]    .push_back(&kNusNormSyst);
  systmap["OscParam"].push_back(&NusOscParamSyst);
  systmap["Stats"]   .push_back(&kNusMCStatsSyst);

  systmap["All"] = getAllNusSysts();

  osc::OscCalcSterile* calc = DefaultSterileCalc(3);

  std::map<std::string, PredictionInterp*> pred123bs;
  std::map<std::string, PredictionInterp*> pred3c3ds;
  std::map<std::string, PredictionCombinePeriods*> preds;

  for(const auto& syst : systmap) {
    std::string axlabel = syst.first;

    pred123bs[axlabel] = new PredictionInterp(syst.second, calc, gen, loaders123b);
    pred3c3ds[axlabel] = new PredictionInterp(syst.second, calc, gen, loaders3c3d);
  }

  PredictionSterile* pred123b = (PredictionSterile*)gen.Generate(loaders123b).release();
  PredictionSterile* pred3c3d = (PredictionSterile*)gen.Generate(loaders3c3d).release();

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

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

  // Two scale factors for scaling spectra
  const double pot_p1p2e3b = kSecondAnaPeriod1POT+kSecondAnaPeriod2POT+kSecondAnaEpoch3bPOT;
  const double pot_e3ce3d  = kSecondAnaEpoch3cPOT+kSecondAnaEpoch3dPOT;

  for(const auto& syst : systmap) {
    std::string axlabel = syst.first;

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

  PredictionCombinePeriods* pred = new PredictionCombinePeriods({
    std::make_pair(pred123b, pot_p1p2e3b),
    std::make_pair(pred3c3d, 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& syst : systmap) {
    std::string axlabel = syst.first;

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

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

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

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

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