mcTruthPredictions.C

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////
/// \brief Macro to create the Monte Carlo truth histograms on
///        all useful weights
/// \author Connor Johnson
/// \date Jan 2020
//////////////////////////////////////////////////////////////////

#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "CAFAna/Systs/EnergySysts2018.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"

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

#include "NDAna/numucc_inc/NumuCCIncAnalysis.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 "TH1.h"
#include "TH2.h"

#include <iostream>
#include <cmath>
#include <sstream>
#include <string>
#include <fstream>
#include <iomanip>
#include <map>
#include <algorithm>

using namespace ana;

const NuTruthVar wgt_xsST = kXSecCVWgt2018_smallerDISScale_NT;
const NuTruthVar wgt_ppfxST = kPPFXFluxCVWgtST;
const NuTruthVar wgtST    = wgt_ppfxST * wgt_xsST;

// Bins
const Binning mukebins_fine = Binning::Simple(200, 0.5, 2.5);
const Binning angvsmukevseavailbins_fine = Binning::Simple(
  mukebins_fine.NBins()*
  angbinsCustom.NBins()*
  eavailbins.NBins(),
  0,
  mukebins_fine.NBins()*
  angbinsCustom.NBins()*
  eavailbins.NBins());
const Binning angvsmukebins_fine = Binning::Simple(
  mukebins_fine.NBins()*
  angbinsCustom.NBins(),
  0,
  mukebins_fine.NBins()*
  angbinsCustom.NBins());
const Binning q2bins_fine = Binning::Simple(500, 0.0, 5.0);
const Binning enubins_fine = Binning::Simple(400, 0.0, 10.0);

const NuTruthVar kTrueMuKEVsCosVsEavail_FineST = Var3D(
  kTrueMuCosthetaST,
  angbinsCustom,
  kTrueMuKEST,
  mukebins_fine,
  kTrueEavailST,
  eavailbins);

const NuTruthVar kTrueMuKEVsCos_FineST = Var2D(
  kTrueMuCosthetaST,
  angbinsCustom,
  kTrueMuKEST,
  mukebins_fine);

const NuTruthHistAxis kTrueMuKEVsCosVsEavail_FineStandardAxisST(
  "True Muon Kinematics",
  angvsmukevseavailbins_fine, 
  kTrueMuKEVsCosVsEavail_FineST
  );

const NuTruthHistAxis kTrueMuKEVsCos_FineStandardAxisST(
  "True Muon Kinematics",
  angvsmukebins_fine,
  kTrueMuKEVsCos_FineST
  );

const NuTruthHistAxis kTrueE_FineStandardAxisST(
  "True Neutrino Energy (GeV)",
  enubins_fine,
  kTrueEST);

const NuTruthHistAxis kTrueQ2_FineStandardAxisST(
  "True Q2 (GeV)",
  q2bins_fine,
  kTrueQ2_NT
  );

map<string, NuTruthVar> weights {
  {"unweight", kNuTruthUnweighted},
  {"ppfx",     wgt_ppfxST},
  {"xsec",     wgt_xsST},
  {"both",     wgtST}
};

void mcTruthPredictions(string filename = "mcPredictions.root")
{
  // Define the spectrum
  SpectrumLoader * loader = new SpectrumLoader("dataset_def_name_newest_snapshot nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata");

  map<string, map<string, Spectrum* > > specmap;

  for (pair<string, NuTruthVar> weight : weights)
  {
    map<string, Spectrum*> specs;
    specs["EAvail"] = new Spectrum(*loader, kTrueMuKEVsCosVsEavail_FineStandardAxisST, kIsTrueSigST,   kNoShift, weight.second);
    specs["MuKin"]  = new Spectrum(*loader, kTrueMuKEVsCos_FineStandardAxisST,         kIsTrueSigST,   kNoShift, weight.second);
    specs["ENu"]    = new Spectrum(*loader, kTrueE_FineStandardAxisST,                 kIsTrueSig1DST, kNoShift, weight.second);
    specs["Q2"]     = new Spectrum(*loader, kTrueQ2_FineStandardAxisST,                kIsTrueSig1DST, kNoShift, weight.second);

    specmap[weight.first] = specs;
  }

  // Go
  loader->Go();

  // Story the outputs
  TFile * outFile = TFile::Open(filename.c_str(), "RECREATE");
  TDirectory * predDir = outFile->mkdir("geniePrediction");
  for(pair<string, map<string, Spectrum* > > weightSpecs : specmap)
  {
    TDirectory * weightDir = predDir->mkdir(weightSpecs.first.c_str());
    for(pair<string, Spectrum*> spec : weightSpecs.second)
      spec.second->SaveTo(weightDir->mkdir(spec.first.c_str()));
  }
  outFile->Close();
}