BlessedPlotsLoad.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Loaders.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/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"

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

using namespace ana;

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

  // Set up some calculators:
  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc3f = DefaultSterileCalc(3);

  // Set up paths to files for loading
  Loaders loaders;

  loaders.SetLoaderPath(fnamenear_concat, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnamefar_concat,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnameswap_concat, caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(fnametau_concat,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders.SetLoaderPath(fnameneardata_concat, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  // Only consider beam with 14 diblocks
  loaders.SetSpillCut(kOnly14DB && kStandardSpillCuts);

  SpectrumLoader floadcos (fnamecos_concat, kCosmic);
  floadcos.SetSpillCut(kOnly14DB && kStandardSpillCuts);

  // Set up ND decompositon for ND spectra plots
  ProportionalDecomp decompNC(
    loaders, kNCAxis, kNusND,
    kNoShift, kNoShift, kTuftsWeightCC
  ); // NC Decomp

  // Set up the sterile prediction
  SterileGenerator gen(
    kNCAxis, kNCBinsNumuCCAxis,
    kNusFD, kNusND, kNumuND,
    kNoShift, kTuftsWeightCC
  );

  // Set up the prediction object
  PredictionSterile* pred = (PredictionSterile*)gen.Generate(
    loaders
  ).release();

  // Set up the cosmic spectrum
  Spectrum sCosmic(floadcos, kNCAxis, kInCosmicTimingWindow && kNusFD);

  // Create predictions necessary for calculating FD purity and efficiency
  PredictionNoExtrap pFDPurEffNum(
    loaders, kNCAxis, kNusFD && kIsNC,
    kNoShift, kTuftsWeightCC
  );
  PredictionNoExtrap pFDPurityDen(
    loaders, kNCAxis, kNusFD,
    kNoShift, kTuftsWeightCC
  );
  PredictionNoExtrap pFDEfcncyDen(
    loaders, kNCAxis, kNusEventQuality && kNusFDFiducial && kIsNC,
    kNoShift, kTuftsWeightCC
  );

  // Create Spectra necessary for calculating ND purity and efficiency
  Spectrum sNDPurEffNum(
    loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
    kNCAxis, kNusND && kIsNC,
    kNoShift, kTuftsWeightCC
  );
  Spectrum sNDPurityDen(
    loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
    kNCAxis, kNusND,
    kNoShift, kTuftsWeightCC
  );
  Spectrum sNDEfcncyDen(
    loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
    kNCAxis, kNusEventQuality && kNusNDFiducial && kIsNC,
    kNoShift, kTuftsWeightCC
  );

  // Fill the spectra!
  loaders.Go();
  floadcos.Go();

  // Two scale factors for scaling spectra, 1 data yr and 3 data yr
  double pot1yr = 6e20;
  double pot3yr = 18e20;

  // Create 2 cosmic spectra, one scaled to each of the above scale factors
  sCosmic.OverridePOT(sCosmic.Livetime()*2.5e18);
  Spectrum sCosmic1y = sCosmic.FakeData(pot1yr);
  Spectrum sCosmic3y = sCosmic.FakeData(pot3yr);

  // Create mock data predictions
  Spectrum sFDPred = pred->Predict(calc3f); // Create an FD spectrum assuming 3 flavor osc
  // Scale the FD beam to 1 data year, add 1 data yr cosmics
  Spectrum sFDData1y = sFDPred.FakeData(pot1yr) + sCosmic1y;
  Spectrum sFDData3y = sFDPred.FakeData(pot3yr) + sCosmic3y; // Repeat for 3 data yr

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

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

  // Save all of the objects
  saveDir = rootF->mkdir("decompNC");
  decompNC    .SaveTo(saveDir);
  saveDir = rootF->mkdir("prediction");
  pred       ->SaveTo(saveDir);
  saveDir = rootF->mkdir("sData1Yr");
  sFDData1y   .SaveTo(saveDir);
  saveDir = rootF->mkdir("sData3Yr");
  sFDData3y   .SaveTo(saveDir);
  saveDir = rootF->mkdir("sCos1Yr");
  sCosmic1y   .SaveTo(saveDir);
  saveDir = rootF->mkdir("sCos3Yr");
  sCosmic3y   .SaveTo(saveDir);
  saveDir = rootF->mkdir("pFDPENum");
  pFDPurEffNum.SaveTo(saveDir);
  saveDir = rootF->mkdir("pFDPDen");
  pFDPurityDen.SaveTo(saveDir);
  saveDir = rootF->mkdir("pFDEDen");
  pFDEfcncyDen.SaveTo(saveDir);
  saveDir = rootF->mkdir("sNDPENum");
  sNDPurEffNum.SaveTo(saveDir);
  saveDir = rootF->mkdir("sNDPDen");
  sNDPurityDen.SaveTo(saveDir);
  saveDir = rootF->mkdir("sNDEDen");
  sNDEfcncyDen.SaveTo(saveDir);

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