SystsCalibAbsLoad17.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 SystsCalibAbsLoad17()
{
  TH1::AddDirectory(0);


   const Var kWgtP3 = kXSecCVWgt2017*kPPFXFluxCVWgt;   
  

  const std::string fnameneardata_conc="prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";
  SpectrumLoader floadneardata(fnameneardata_conc);
  floadneardata.SetSpillCut(kStandardSpillCuts);

  std::string fnear_ns="defname: prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_no-shift_v1";
  std::string ffar_ns="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_lightmodel-noshift_v2";
  std::string fswap_ns="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_lightmodel-noshift_v2";
  std::string ftau_ns="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_lightmodel-noshift_v2";
  
  SpectrumLoader floadnear_ns(fnear_ns);
   SpectrumLoader floadfar_ns(ffar_ns);
   SpectrumLoader floadswap_ns(fswap_ns);
   SpectrumLoader floadtau_ns(ftau_ns);
  floadnear_ns    .SetSpillCut(kStandardSpillCuts);
  floadfar_ns     .SetSpillCut(kStandardSpillCuts);
  floadswap_ns    .SetSpillCut(kStandardSpillCuts);
  floadtau_ns     .SetSpillCut(kStandardSpillCuts);

  



   std::string Llevelu="defname:  prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1";
  std::string Lleveld="defname:  prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1";
    std::string cherenkov="defname:  prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1";



 std::string Llevelu_fd="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v2";
  std::string Llevelu_swp="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v2";
  std::string Llevelu_tau="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_lightmodel-lightup-calibdown_v2";

  std::string Lleveld_fd="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v2";
  std::string Lleveld_swp="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v2";
  std::string Lleveld_tau="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_lightmodel-lightdown-calibup_v2";

  
  std::string cherenkovfd="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v2";
  std::string cherenkovswp="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_ckv-proton-shift-down_v2";
  std::string cherenkovtau="defname:  prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_ckv-proton-shift-down_v2";



  std::map<std::string, std::vector<std::string> > shiftlabels;
  shiftlabels["ULLevel"] = std::vector<std::string>();
  shiftlabels["DLLevel"] = std::vector<std::string>();
  shiftlabels["Cherenkov"] = std::vector<std::string>();



  shiftlabels["ULLevel"].push_back(Llevelu);
  shiftlabels["ULLevel"].push_back(Llevelu_fd);
  shiftlabels["ULLevel"].push_back(Llevelu_swp);
  shiftlabels["ULLevel"].push_back(Llevelu_tau);

  shiftlabels["DLLevel"].push_back(Lleveld);
  shiftlabels["DLLevel"].push_back(Lleveld_fd);
  shiftlabels["DLLevel"].push_back(Lleveld_swp);
  shiftlabels["DLLevel"].push_back(Lleveld_tau);

  shiftlabels["Cherenkov"].push_back(cherenkov);
  shiftlabels["Cherenkov"].push_back(cherenkovfd);
  shiftlabels["Cherenkov"].push_back(cherenkovswp);
  shiftlabels["Cherenkov"].push_back(cherenkovtau);

/*  shiftlabels["Noshift"].push_back(fnear_ns);
  shiftlabels["Noshift"].push_back(ffar_ns);
  shiftlabels["Noshift"].push_back(fswap_ns);
  shiftlabels["Noshift"].push_back(ftau_ns);
*/

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

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

  //std::cout << "Hey I am here "<<std::cout;

  // 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) {
    decomps_nominal[sample.first] = new ProportionalDecomp(
      floadnear_ns, floadneardata,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kNoShift, kWgtP3
    );
    
     predsNE_nominal[sample.first] = new PredictionNoExtrap(
      floadfar_ns, floadswap_ns, floadtau_ns,
      labelRecoE, kNCDisappearanceEnergyBinning,kNusNDEnergy17,*sample.second.second.second,
      kNoShift, kWgtP3
    );



    ProportionalDecomp* decompnc_n   = new ProportionalDecomp(
      floadnear_ns, floadneardata,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kNoShift, kWgtP3
    );
    ProportionalDecomp* decompnumu_n = new ProportionalDecomp(
      floadnear_ns, floadneardata,
      kNCBinsNumuCCAxis, kNumuND,
      kNoShift, kNoShift, kWgtP3
    );
    ModularExtrapSterile* extrap_n = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
      floadnear_ns, floadswap_ns, floadfar_ns, floadtau_ns, *decompnc_n, *decompnumu_n,
      *sample.second.first, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first, kNumuND,
      kNoShift, kWgtP3
    ));
    predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);

    for(const auto& systlabel : systematics) {
      int ndpos = 0;
      int fdpos = 1;
      int fspos = 2;
      int ftpos=3;
      decomps_shifted[sample.first][systlabel.first][1] = new ProportionalDecomp(
        *systlabel.second[ndpos],floadneardata,
        *sample.second.first, *sample.second.second.first,
        kNoShift, kNoShift, kWgtP3
      );
      predsNE_shifted[sample.first][systlabel.first][1] = new PredictionNoExtrap(
        *systlabel.second[fdpos], *systlabel.second[fspos], *systlabel.second[ftpos],
        labelRecoE, kNCDisappearanceEnergyBinning ,
        kNusNDEnergy17, *sample.second.second.second,
        kNoShift, kWgtP3
      );

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

  for(const auto& systlabel : systematics) {
    for(const auto& loader : systlabel.second) {
      loader->Go();
    }
  }
  floadneardata.Go();
  floadnear_ns.Go();
  floadfar_ns .Go();
  floadswap_ns.Go();
  floadtau_ns.Go();

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

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