SystsCalibRelLoad.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 "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/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSSystFunctions.h"

using namespace ana;

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

  SpectrumLoader floadnear(fnamenear_concat);
  SpectrumLoader floadfar( fnamefar_concat);
  SpectrumLoader floadswap(fnameswap_concat);
  SpectrumLoader floadtau( fnametau_concat);
  SpectrumLoader floadneardata(fnameneardata_concat);

  floadnear    .SetSpillCut(kOnly14DB && kStandardSpillCuts);
  floadfar     .SetSpillCut(kOnly14DB && kStandardSpillCuts);
  floadswap    .SetSpillCut(kOnly14DB && kStandardSpillCuts);
  floadtau     .SetSpillCut(kOnly14DB && kStandardSpillCuts);
  floadneardata.SetSpillCut(kOnly14DB && 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_R16-03-03-prod2reco.a_nd_genie_nonswap_nogenierw_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1");
  shiftlabels["CalFlat095FD"].push_back("prod_caf_R16-03-03-prod2reco.a_fd_genie_nonswap_fhc_nova_v08_full_calib-shift-fd-xy-neg-offset_v1");
  shiftlabels["CalFlat095FD"].push_back("prod_caf_R16-03-03-prod2reco.a_fd_genie_fluxswap_fhc_nova_v08_full_calib-shift-fd-xy-neg-offset_v1");
  shiftlabels["CalFlat105ND"].push_back("prod_caf_R16-03-03-prod2reco.a_nd_genie_nonswap_nogenierw_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1");
  shiftlabels["CalFlat105FD"].push_back("prod_caf_R16-03-03-prod2reco.d_fd_genie_nonswap_fhc_nova_v08_full_calib-shift-fd-xy-pos-offset_v1");
  shiftlabels["CalFlat105FD"].push_back("prod_caf_R16-03-03-prod2reco.d_fd_genie_fluxswap_fhc_nova_v08_full_calib-shift-fd-xy-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(kOnly14DB && 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 = "Calorimetric Energy (GeV)";
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  std::pair<Cut*, Cut*> Ana01(new Cut(kNusND), new Cut(kNusFD));
  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
    new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kCaloE),
    Ana01
  );

  // 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, kTuftsWeightCC
    );
    predsNE_nominal[sample.first] = new PredictionNoExtrap(
      floadfar, floadswap, floadtau,
      sample.second.first->label, sample.second.first->bins,
      sample.second.first->var, *sample.second.second.second,
      kNoShift, kTuftsWeightCC
    );

    ProportionalDecomp* decompnc_n   = new ProportionalDecomp(
      floadnear, floadneardata,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kNoShift, kTuftsWeightCC
    );
    ProportionalDecomp* decompnumu_n = new ProportionalDecomp(
      floadnear, floadneardata,
      kNCBinsNumuCCAxis, kNumuND,
      kNoShift, kNoShift, kTuftsWeightCC
    );
    ModularExtrapSterile* extrap_n = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
      floadnear, floadswap, floadfar, floadtau, *decompnc_n, *decompnumu_n,
      *sample.second.first, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first, kNumuND,
      kNoShift, kTuftsWeightCC
    ));
    predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);

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

      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, kTuftsWeightCC
        );
        predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
          floadfar, floadswap, floadtau,
          sample.second.first->label, sample.second.first->bins,
          sample.second.first->var, *sample.second.second.second,
          kNoShift, kTuftsWeightCC
        );

        ProportionalDecomp* decompnc_s   = new ProportionalDecomp(
          *systlabel.second[0], floadneardata,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kNoShift, kTuftsWeightCC
        );
        ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
          *systlabel.second[0], floadneardata,
          kNCBinsNumuCCAxis, kNumuND,
          kNoShift, kNoShift, kTuftsWeightCC
        );
        ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
          *systlabel.second[0], floadswap, floadfar,
          floadtau, *decompnc_s, *decompnumu_s,
          *sample.second.first, kNCBinsNumuCCAxis,
          *sample.second.second.second, *sample.second.second.first, kNumuND,
          kNoShift, kTuftsWeightCC
        ));
        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, kTuftsWeightCC
        );
        predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
          *systlabel.second[fdpos], *systlabel.second[fspos], floadtau,
          sample.second.first->label, sample.second.first->bins,
          sample.second.first->var, *sample.second.second.second,
          kNoShift, kTuftsWeightCC
        );

        ProportionalDecomp* decompnc_s   = new ProportionalDecomp(
          floadnear, floadneardata,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kNoShift, kTuftsWeightCC
        );
        ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
          floadnear, floadneardata,
          kNCBinsNumuCCAxis, kNumuND,
          kNoShift, kNoShift, kTuftsWeightCC
        );
        ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
          floadnear, *systlabel.second[fspos], *systlabel.second[fdpos],
          floadtau, *decompnc_s, *decompnumu_s,
          *sample.second.first, kNCBinsNumuCCAxis,
          *sample.second.second.second, *sample.second.second.first, kNumuND,
          kNoShift, kTuftsWeightCC
        ));
        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 = "/nova/ana/steriles/Ana01/Systematics/";
  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
  );
}