MakeNus17Prediction.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "NuXAna/Cuts/NusCuts17.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 "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 "CAFAna/Vars/PPFXWeights.h"
#include "NuXAna/macros/Nus17/NusLoadProd3.h"

#include "OscLib/OscCalcSterile.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include <utility>

using namespace ana;

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

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

  
  // ND DATA + MC
  SpectrumLoader loaderNDMC(fNDmcDECAF);
  SpectrumLoader loaderNDdata(fNDdataDECAF);
  loaderNDMC.SetSpillCut(kStandardSpillCuts);
  loaderNDdata.SetSpillCut(kStandardSpillCuts);
  
  SpectrumLoader* loaderFDMC_non = new SpectrumLoader(fFDmcDECAF);
  SpectrumLoader* loaderFDMC_swp = new SpectrumLoader(fFDswapDECAF);
  SpectrumLoader* loaderFDMC_tau = new SpectrumLoader(fFDtauDECAF);
  loaderFDMC_non->SetSpillCut(kStandardSpillCuts);
  loaderFDMC_swp->SetSpillCut(kStandardSpillCuts);
  loaderFDMC_tau->SetSpillCut(kStandardSpillCuts);

  // FD DATA COSMICS
  SpectrumLoader loaderFDCOS(fFDcosmicdataDECAF, kCosmic);
  loaderFDCOS.SetSpillCut(kStandardSpillCuts);

  const Binning kNCDisappearanceEnergyBinningExt1 = Binning::Simple(100, 0, 10);
  const Binning kNCDisappearanceEnergyBinningExt2 = Binning::Simple(40, 0, 4);

  const HistAxis kNCNDAxisExtended1("Visible Energy (GeV)",
                                   kNCDisappearanceEnergyBinningExt1, kNus17Energy);
  const HistAxis kNCFDAxisExtended1("Visible Energy (GeV)",
                                   kNCDisappearanceEnergyBinningExt1, kNus17Energy);
  const HistAxis kNCBinsNumuCCAxisExt1("Reconstructed Neutrino Energy (GeV)",
                                       kNCDisappearanceEnergyBinningExt1, kCCE );

  const HistAxis kNCNDAxisExtended2("Visible Energy (GeV)",
                                   kNCDisappearanceEnergyBinningExt2, kNus17Energy);
  const HistAxis kNCFDAxisExtended2("Visible Energy (GeV)",
                                   kNCDisappearanceEnergyBinningExt2, kNus17Energy);
  const HistAxis kNCBinsNumuCCAxisExt2("Reconstructed Neutrino Energy (GeV)",
                                       kNCDisappearanceEnergyBinningExt2, kCCE );


  
  // Set up NC decompositon for ND spectra plots
  ProportionalDecomp decompNC1(
    loaderNDMC, loaderNDdata, kNCNDAxisExtended1, kNus17ND,
    kNoShift, kNoShift, kReweight
  ); // NC Decomp

  // Set up ND decompositon for Numu
  ProportionalDecomp decompNumu1(
    loaderNDMC, loaderNDdata, kNCBinsNumuCCAxisExt1, kNumuND,
    kNoShift, kNoShift, kReweight
  ); // NC Decomp

  ModularExtrapSterile extrap1 = ModularExtrapSterile::NCDisappearance(
    loaderNDMC, *loaderFDMC_swp, *loaderFDMC_non, *loaderFDMC_tau,
    decompNC1, decompNumu1,
    kNCFDAxisExtended1, kNCBinsNumuCCAxisExt1,
    kNus17FDAlt, kNus17ND, kNumuND,
    kNoShift, kReweight
  );

  // Test what 0.5 GeV to 4 GeV looks (for energy scaled)
  const Cut kECut_scale = kNus17Energy > 0.31884 && kNus17Energy < 2.5509;
  
  ProportionalDecomp decompNC2(
    loaderNDMC, loaderNDdata, kNCNDAxisExtended2, kNus17ND && kECut_scale,
    kNoShift, kNoShift, kReweight
  ); // NC Decomp

  // Set up ND decompositon for Numu
  ProportionalDecomp decompNumu2(
    loaderNDMC, loaderNDdata, kNCBinsNumuCCAxisExt2, kNumuND && kECut_scale,
    kNoShift, kNoShift, kReweight
  ); // NC Decomp

  ModularExtrapSterile extrap2 = ModularExtrapSterile::NCDisappearance(
    loaderNDMC, *loaderFDMC_swp, *loaderFDMC_non, *loaderFDMC_tau,
    decompNC2, decompNumu2,
    kNCFDAxisExtended2, kNCBinsNumuCCAxisExt2,
    kNus17FDAlt && kECut_scale, kNus17ND && kECut_scale, kNumuND && kECut_scale,
    kNoShift, kReweight
  );

  // Set up the prediction object
  PredictionSterile* pred1 = new PredictionSterile(&extrap1);
  PredictionSterile* pred2 = new PredictionSterile(&extrap2);

  // Set up the cosmic spectrum
  Spectrum sCosmic1(loaderFDCOS, kNCFDAxisExtended1, kInCosmicTimingWindow && kNus17FDAlt && kECut);
  Spectrum sCosmic2(loaderFDCOS, kNCFDAxisExtended2, kInCosmicTimingWindow && kNus17FDAlt && kECut_scale);

  // Fill the spectra!
  loaderNDMC.Go();
  loaderNDdata.Go();
  loaderFDMC_non->Go();
  loaderFDMC_swp->Go();
  loaderFDMC_tau->Go();
  loaderFDCOS.Go();

  // 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("decompNC1");
  decompNC1.SaveTo(saveDir);
  saveDir = rootF->mkdir("decompNumu1");
  decompNumu1.SaveTo(saveDir);
  saveDir = rootF->mkdir("extrap1");
  extrap1.SaveTo(saveDir);
  saveDir = rootF->mkdir("prediction1");
  pred1->SaveTo(saveDir);
  saveDir = rootF->mkdir("sCos1");
  sCosmic1.SaveTo(saveDir);


  saveDir = rootF->mkdir("decompNC2");
  decompNC2.SaveTo(saveDir);
  saveDir = rootF->mkdir("decompNumu2");
  decompNumu2.SaveTo(saveDir);
  saveDir = rootF->mkdir("extrap2");
  extrap2.SaveTo(saveDir);
  saveDir = rootF->mkdir("prediction2");
  pred2->SaveTo(saveDir);
  saveDir = rootF->mkdir("sCos2");
  sCosmic2.SaveTo(saveDir);

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