specprod_numuccinc.C

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////
/// \brief Macro for numucc_inc for running on NERSC
///        
/// \author Connor Johnson
/// \date Sept 2019
//////////////////////////////////////////////////////////////////

#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"

#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "NDAna/numucc_inc/NumuCCIncBins.h"
#include "NDAna/numucc_inc/NumuCCIncVars.h"
#include "NDAna/numucc_inc/NumuCCIncCuts.h"
#include "NDAna/numucc_inc/NDXSecMuonPID.h"
#include "NDAna/numucc_inc/UnfoldingVariable.h"

#include "TFile.h"

#include <iostream>
#include <string>
#include <map>

const std::string NOMINAL_DATASET = "nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata";
const unsigned int N_GENIE_UNIVERSES = 1000;
const unsigned int N_PPFX_UNIVERSES = 100;

const ana::NuTruthVar wgt_xsST = ana::kXSecCVWgt2018_smallerDISScale_NT;
const ana::Var wgt_xs          = ana::VarFromNuTruthVar(wgt_xsST, 1);
const ana::NuTruthVar wgtST    = ana::kPPFXFluxCVWgtST*wgt_xsST;
const ana::Var wgt             = ana::VarFromNuTruthVar(wgtST, 1);

void specprod_numuccinc(std::string outfileName, int firstGenie, int lastGenie, int firstPPFX = -1, int lastPPFX = -1)
{
  // Maps for the UnfoldingVariable class. Store in order of [(uint) iuniv][(string) "varname"]
  std::map<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > > genieUnfoldingVars;
  std::map<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > > ppfxUnfoldingVars;

  // Create the nominal spectrum loader
  ana::SpectrumLoader * nominalLoader = new ana::SpectrumLoader(NOMINAL_DATASET.c_str());
  nominalLoader->SetSpillCut(ana::kStandardSpillCuts);

  if ((firstGenie < 0 || lastGenie < firstGenie) && (firstPPFX < 0 || lastPPFX < firstPPFX))
  {
    cout << "Invalid choices for both genie and ppfx universes. Will not create spectrums.\nExiting..." << endl;
    exit(1);
  }

  // Load the genie, use firstGenie < 0 to skip them all
  if (firstGenie >= 0 && lastGenie >= firstGenie)
  {
    cout << "Running " << to_string(lastGenie - firstGenie + 1) << " genie universes. (" << to_string(firstGenie) << "-"  << to_string(lastGenie) << ")" << endl;
    ana::GenieMultiverseParameters verse(N_GENIE_UNIVERSES, ana::getAllXsecSysts_2018_RPAFix());
    std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();

    for(int igenie = firstGenie; igenie <= lastGenie; ++igenie)
    {
      genieUnfoldingVars[igenie]["EAvail"] = new ana::xsec::UnfoldingVariable(nominalLoader,
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, &genie_shifts[igenie], &wgtST);

      genieUnfoldingVars[igenie]["ENu"]    = new ana::xsec::UnfoldingVariable(nominalLoader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &genie_shifts[igenie], &wgtST);

      genieUnfoldingVars[igenie]["Q2"]     = new ana::xsec::UnfoldingVariable(nominalLoader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &genie_shifts[igenie], &wgtST);
    }
  }

  // Load the ppfx, use firstPPFX < 0 to skip them all
  if (firstPPFX >= 0 && lastPPFX >= firstPPFX)
  {
    cout << "Running " << to_string(lastPPFX - firstPPFX + 1) << " ppfx universes. (" << to_string(firstPPFX) << "-" << to_string(lastPPFX) << ")" << endl;
    const std::vector<ana::NuTruthVar> vppfxST = ana::GetkPPFXFluxUnivWgtST();

    for(int ippfx = firstPPFX; ippfx <= lastPPFX; ++ippfx)
    {
      const ana::NuTruthVar * xsTimesPPFXi = new ana::NuTruthVar(wgt_xsST * vppfxST[ippfx]);

      ppfxUnfoldingVars[ippfx]["EAvail"] = new ana::xsec::UnfoldingVariable(nominalLoader,
        &ana::kTrueMuKEVsCosVsEavailStandardAxisST, &ana::kRecoMuKEVsCosVsEavailStandardAxis,
        &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, &ana::kNoShift, xsTimesPPFXi);

      ppfxUnfoldingVars[ippfx]["ENu"]    = new ana::xsec::UnfoldingVariable(nominalLoader, 
        &ana::kTrueEStandardAxisST, &ana::kRecoEStandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, xsTimesPPFXi);

      ppfxUnfoldingVars[ippfx]["Q2"]     = new ana::xsec::UnfoldingVariable(nominalLoader,
        &ana::kTrueQ2StandardAxisST, &ana::kRecoQ2StandardAxis,
        &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &ana::kNoShift, xsTimesPPFXi);
    }
  }

  // Fill the spectrums
  nominalLoader->Go();

  // Store it all.
  TFile* fOut = new TFile(outfileName.c_str(), "RECREATE"); 
  fOut->cd();
  TDirectory* dir;

  // Save the genie universes (carefully grabbing by references for slight optimization)
  for (const std::pair<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > >& genieUnivVars : genieUnfoldingVars)
  {
    unsigned int genieUniv    = genieUnivVars.first;
    std::string  genieUnivStr = to_string(genieUniv);
    std::string  folderName   = "genie" + genieUnivStr;
    const std::map<std::string, ana::xsec::UnfoldingVariable*>& unfoldVars = genieUnivVars.second;

    dir = fOut->mkdir(folderName.c_str());
    dir->cd();
    for (const std::pair<std::string, ana::xsec::UnfoldingVariable*>& genieUnivVar : unfoldVars)
    {
      std::string varName     = genieUnivVar.first;
      ana::xsec::UnfoldingVariable * var = genieUnivVar.second;

      TDirectory * varDir = dir->mkdir(varName.c_str());
      var->SaveSpectrums(varDir);
    }
  }

  // Save the ppfx universes
  for (const std::pair<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > >& ppfxUnivVars : ppfxUnfoldingVars)
  {
    unsigned int ppfxUniv    = ppfxUnivVars.first;
    std::string  ppfxUnivStr = to_string(ppfxUniv);
    std::string  folderName  = "ppfx" + ppfxUnivStr;
    const std::map<std::string, ana::xsec::UnfoldingVariable*>& unfoldVars = ppfxUnivVars.second;

    dir = fOut->mkdir(folderName.c_str());
    dir->cd();
    for (const std::pair<std::string, ana::xsec::UnfoldingVariable*>& ppfxUnivVar : unfoldVars)
    {
      std::string varName     = ppfxUnivVar.first;
      ana::xsec::UnfoldingVariable * var = ppfxUnivVar.second;

      TDirectory * varDir = dir->mkdir(varName.c_str());
      var->SaveSpectrums(varDir);
    }
  }
  fOut->Close();
  return;
}