SystsAcceptanceLoad17.C

Go to the documentation of this file.

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSSystFunctions.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "NuXAna/Vars/HistAxes.h"

using namespace ana;

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

  Loaders loaders;
   const Var kWgtP31 = kXSecCVWgt2017*kPPFXFluxCVWgt;

/*     std::string fnear="defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with limit 50";
  std::string ffar ="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1 with limit 50";
  std::string fswap="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1 with limit 50";
  std::string ftau="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1 with limit 50";
   std::string fnameneardata_con="defname: prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns with limit 50";
*/
   std::string fnear="prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
   std::string ffar ="prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
   std::string fswap="prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
   std::string ftau="prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";
   std::string fnameneardata_con="prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";



  loaders.SetLoaderPath(fnear, caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(ffar,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fswap, caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(ftau,  caf::kFARDET,  Loaders::kMC, ana::kBeam, Loaders::kTauSwap);
  loaders.SetLoaderPath(fnameneardata_con, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  loaders.SetSpillCut(kStandardSpillCuts);

  // A map to set up systematics, starting with a label for the systematic,
  // and finishing with a vector for items necessary to set up the Syst object
  std::map<std::string, std::string> shiftlabels;
  shiftlabels["ND1"] = "ND1";
  shiftlabels["ND2"] = "ND2";
  shiftlabels["ND3"] = "ND3";
  shiftlabels["ND4"] = "ND4";

  const Cut kNDFirstBlock(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return false; }
      if(sr->vtx.elastic.vtx.X() > NDR) { return false; }
      if(sr->vtx.elastic.vtx.X() < NDL ) { return false; }
      if(sr->vtx.elastic.vtx.Y() > NDT) { return false; }
      if(sr->vtx.elastic.vtx.Y() <0) { return false; }
      if(sr->vtx.elastic.vtx.Z() < 0) { return false; }
       if(sr->vtx.elastic.vtx.Z() > 635) { return false; }
      return true;
    });

    const Cut kNDSecondBlock(
    [](const caf::SRProxy* sr)
    {
      if( !sr->vtx.elastic.IsValid) { return false; }
      if(sr->vtx.elastic.vtx.X() > NDR) { return false; }
      if(sr->vtx.elastic.vtx.X() < NDL ) { return false; }
      if(sr->vtx.elastic.vtx.Y() > NDT) { return false; }
      if(sr->vtx.elastic.vtx.Y() <0) { return false; }
      if(sr->vtx.elastic.vtx.Z() >1270) { return false; }
       if(sr->vtx.elastic.vtx.Z() < 635) { return false; }
      return true;
    });
   
     const Cut kNDThirdBlock(
    [](const caf::SRProxy* sr)
    {
      if( !sr->vtx.elastic.IsValid) { return false; }
      if(sr->vtx.elastic.vtx.X() > NDR) { return false; }
      if(sr->vtx.elastic.vtx.X() < NDL ) { return false; }
      if(sr->vtx.elastic.vtx.Y() > 0) { return false; }
      if(sr->vtx.elastic.vtx.Y() <NDB) { return false; }
      if(sr->vtx.elastic.vtx.Z() >1270) { return false; }
       if(sr->vtx.elastic.vtx.Z() < 635) { return false; }
      return true;
    });

     const Cut kNDFourthBlock(
    [](const caf::SRProxy* sr)
    {
      if(!sr->vtx.elastic.IsValid) { return false; }
      if(sr->vtx.elastic.vtx.X() > NDR) { return false; }
      if(sr->vtx.elastic.vtx.X() < NDL ) { return false; }
      if(sr->vtx.elastic.vtx.Y() > 0) { return false; }
      if(sr->vtx.elastic.vtx.Y() <NDB) { return false; }
      if(sr->vtx.elastic.vtx.Z() >635) { return false; }
       if(sr->vtx.elastic.vtx.Z() < 0) { return false; }
      return true;
    });



  // A map of systematic objects, indexed by the same label as above
  std::map<std::string, Cut*> systematics;
  systematics["ND1"] = new Cut(kNDFirstBlock);
   systematics["ND2"] = new Cut(kNDSecondBlock);
  systematics["ND3"] = new Cut(kNDThirdBlock);
  systematics["ND4"] = new Cut(kNDFourthBlock);



  // Define a map of samples and cuts/axes to set them up
  // Each sample, essentially, is a way to set up the analysis
  std::string labelRecoE = "Energy Deposited in Scintillator (GeV)";
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  std::pair<Cut*, Cut*> Ana01(new Cut(kNus17ND), new Cut(kNus17FDAlt));
  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
    new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kNusNDEnergy17),
    Ana01
  );

 HistAxis* FDAxis = new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kNusFDEnergy17);
  // Set up maps to the decompositions/predictions (each flavor component)
  // Nominal maps are indexed purely by sample label
  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
  std::map<std::string,                                     IDecomp*>     decomps_nominal;
  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
  std::map<std::string,                                     PredictionNoExtrap*>     predsNE_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
  std::map<std::string,                                     PredictionSterile*>      predsSt_nominal;
  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >  predsSt_shifted;

  // Set up the actual decompositions/predictions
  // sample.first = the sample label
  // shiftlabel.first = the systematic shift label
  // sample.second.first = the sample HistAxis
  // sample.second.second = the sample Cut
  // syst.first = the integer number of sigmas for the shift
  // syst.second = the systematic being shifted
  for(const auto& sample : cut_samples) {
    SterileGenerator gen(
      *FDAxis, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first,
      kNumuND, kNoShift, kWgtP31
    );

    // Set up the nominal decompositions
    // Only one is needed per sample
    decomps_nominal[sample.first] = new CheatDecomp(
      loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
      *sample.second.first, *sample.second.second.first,
      kNoShift, kWgtP31
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      loaders,
      *FDAxis, *sample.second.second.second,
      kNoShift, kWgtP31
    );
    predsSt_nominal[sample.first] = (PredictionSterile*)gen.Generate(
      loaders
    ).release();

    // Set up the shifted decompositions
    // One is needed for every sample, systematic, and sigma level
    for(const auto& shiftlabel : shiftlabels) {
      decomps_shifted[sample.first][shiftlabel.first][1] = new CheatDecomp(
        loaders.GetLoader(caf::kNEARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
        *sample.second.first, *sample.second.second.first && *systematics[shiftlabel.first],
        kNoShift, kWgtP31
      );
      predsNE_shifted[sample.first][shiftlabel.first][1] = new PredictionNoExtrap(
        loaders,
        *FDAxis, *sample.second.second.second,
        kNoShift, kWgtP31
      );

      SterileGenerator gen(
        *FDAxis, kNCBinsNumuCCAxis,
        *sample.second.second.second,
        *sample.second.second.first && *systematics[shiftlabel.first],
        kNumuND && *systematics[shiftlabel.first],
        kNoShift, kWgtP31
      );
      predsSt_shifted[sample.first][shiftlabel.first][1] =
        (PredictionSterile*)gen.Generate(loaders).release();
    }
  }

  loaders.Go();

  std::string folder = "./";
  std::string filenm = "SystsNDAcceptance17";

  std::string fullLocation = folder + filenm + ".root";
  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");

  SaveMaps(rootF,
    decomps_nominal, decomps_shifted,
    predsNE_nominal, predsNE_shifted,
    predsSt_nominal, predsSt_shifted
  );
}