SystsCalibRelLoad17.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.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 "CAFAna/Vars/PPFXWeights.h"
#include "NuXAna/Vars/HistAxes.h"

using namespace ana;

void SystsCalibRelLoad17()
{
  TH1::AddDirectory(0);
  const Var kWgtP31 = kXSecCVWgt2017*kPPFXFluxCVWgt;


  std::string fnamenear_con="prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string fnameneardata_con="prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";
  std::string fFDMC_nons="prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string fFDMC_swps="prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
  std::string fFDMC_taus="prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";

/*    std::string fnamenear_con="defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with limit 1";
    std::string fnameneardata_con="defname: prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns with limit 1";
  std::string fFDMC_nons="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1 with limit 1";
  std::string fFDMC_swps="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1 with limit 1";
  std::string fFDMC_taus="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1 with limit 1";
*/


  SpectrumLoader floadnear(fnamenear_con);
  SpectrumLoader floadfar(fFDMC_nons);
  SpectrumLoader floadswap(fFDMC_swps);
  SpectrumLoader floadtau(fFDMC_taus);
  SpectrumLoader floadneardata(fnameneardata_con);

  floadnear    .SetSpillCut(kStandardSpillCuts);
  floadfar     .SetSpillCut(kStandardSpillCuts);
  floadswap    .SetSpillCut(kStandardSpillCuts);
  floadtau     .SetSpillCut(kStandardSpillCuts);
  floadneardata.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::vector<std::string> > shiftlabels;
  shiftlabels["CalFlat095ND"] = std::vector<std::string>();
  shiftlabels["CalFlat105ND"] = std::vector<std::string>();
  shiftlabels["CalFlat095FD"] = std::vector<std::string>();
  shiftlabels["CalFlat105FD"] = std::vector<std::string>();
  shiftlabels["CalFlat095ND"].push_back("prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-neg-offset_v1");
  shiftlabels["CalFlat095FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_calib-shift-fd-xyview-neg-offset_v1");
  shiftlabels["CalFlat095FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_calib-shift-fd-xyview-neg-offset_v1");
  shiftlabels["CalFlat095FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_calib-shift-fd-xyview-neg-offset_v1");


  shiftlabels["CalFlat105ND"].push_back("prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-pos-offset_v1");
  shiftlabels["CalFlat105FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_calib-shift-fd-xyview-pos-offset_v1");
  shiftlabels["CalFlat105FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_calib-shift-fd-xyview-pos-offset_v1");
  shiftlabels["CalFlat105FD"].push_back("prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_calib-shift-fd-xyview-pos-offset_v1");
  // A map of systematic objects, indexed by the same label as above
  std::map<std::string, std::vector<SpectrumLoader*> > systematics;

  for(const auto& shiftlabel : shiftlabels) {
    systematics[shiftlabel.first] = std::vector<SpectrumLoader*>();
    for(const auto& samdef : shiftlabel.second){
      systematics[shiftlabel.first].push_back(new SpectrumLoader(samdef));
      systematics[shiftlabel.first].back()->SetSpillCut(kStandardSpillCuts);
    }
  }

  // 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 = "Visible Energy (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
  for(const auto& sample : cut_samples) {
    // Set up the nominal decompositions
    // Only one is needed per sample
    decomps_nominal[sample.first] = new CheatDecomp(
      floadnear,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kWgtP31
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      floadfar, floadswap, floadtau,
      labelRecoE, kNCDisappearanceEnergyBinning,kNusFDEnergy17, *sample.second.second.second,
      kNoShift, kWgtP31
    );

    ProportionalDecomp* decompnc_n   = new ProportionalDecomp(
      floadnear, floadneardata,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kNoShift, kWgtP31
    );
    ProportionalDecomp* decompnumu_n = new ProportionalDecomp(
      floadnear, floadneardata,
      kNCBinsNumuCCAxis, kNumuND,
      kNoShift, kNoShift, kWgtP31
    );
    ModularExtrapSterile* extrap_n = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
      floadnear, floadswap, floadfar, floadtau, *decompnc_n, *decompnumu_n,
      *FDAxis, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first, kNumuND,
      kNoShift, kWgtP31
    ));
    predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);
   
    //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

    // Set up the shifted decompositions
    // One is needed for every sample, systematic, and sigma level
    //
    for(const auto& systlabel : systematics) {
      int fdpos = 0;
      int fspos = 1;
      int ftpos = 2;
      if(systlabel.first.find("ND") != std::string::npos) {
        decomps_shifted[sample.first][systlabel.first][1] = new CheatDecomp(
          *systlabel.second[0],
          *sample.second.first, *sample.second.second.first,
          kNoShift, kWgtP31
        );
        predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
          floadfar, floadswap, floadtau,
          labelRecoE, kNCDisappearanceEnergyBinning,kNusFDEnergy17, *sample.second.second.second,
          kNoShift, kWgtP31
        );

        ProportionalDecomp* decompnc_s   = new ProportionalDecomp(
          *systlabel.second[0], floadneardata,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kNoShift, kWgtP31
        );
        ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
          *systlabel.second[0], floadneardata,
          kNCBinsNumuCCAxis, kNumuND,
          kNoShift, kNoShift, kWgtP31
        );
        ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
          *systlabel.second[0], floadswap, floadfar,
          floadtau, *decompnc_s, *decompnumu_s,
          *FDAxis, kNCBinsNumuCCAxis,
          *sample.second.second.second, *sample.second.second.first, kNumuND,
          kNoShift, kWgtP31
        ));
        predsSt_shifted[sample.first][systlabel.first][1] = new PredictionSterile(
          extrap_s
        );
      }
      else {
        decomps_shifted[sample.first][systlabel.first][1] = new CheatDecomp(
          floadnear,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kWgtP31
        );
        predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
          *systlabel.second[fdpos], *systlabel.second[fspos], *systlabel.second[ftpos],
          labelRecoE, kNCDisappearanceEnergyBinning,kNusFDEnergy17, *sample.second.second.second,
          kNoShift, kWgtP31
        );

        ProportionalDecomp* decompnc_s   = new ProportionalDecomp(
          floadnear, floadneardata,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kNoShift, kWgtP31
        );
        ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
          floadnear, floadneardata,
          kNCBinsNumuCCAxis, kNumuND,
          kNoShift, kNoShift, kWgtP31
        );
        ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
          floadnear, *systlabel.second[fspos], *systlabel.second[fdpos],
          *systlabel.second[ftpos], *decompnc_s, *decompnumu_s,
          *FDAxis, kNCBinsNumuCCAxis,
          *sample.second.second.second, *sample.second.second.first, kNumuND,
          kNoShift, kWgtP31
        ));
        predsSt_shifted[sample.first][systlabel.first][1] = new PredictionSterile(
          extrap_s
        );
      }
    }
  }

  for(const auto& systlabel : systematics) {
    for(const auto& loader : systlabel.second) {
      loader->Go();
    }
  }
  floadnear.Go();
  floadfar .Go();
  floadswap.Go();
  floadtau .Go();
  floadneardata.Go();

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

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