makeRealDataFluxes.C

Go to the documentation of this file.

// makeFakeData.C

#include "CAFAna/Core/Cut.h"
#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/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 "TFile.h"
#include "TH1.h"
#include "TH2.h"

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

using namespace ana;

std::map<std::string, ana::HistAxis> unfoldingVars{
  {"data_eavail", kRecoMuKEVsCosVsEavailStandardAxis},
  {"data_enu",    kRecoEStandardAxis},
  {"data_q2",     kRecoQ2StandardAxis}
};

std::map<std::string, ana::Cut> recoCuts{
  {"data_eavail", kAllNumuCCCuts},
  {"data_enu",    kAllNumuCC1DCuts},
  {"data_q2",     kAllNumuCC1DCuts}
};

const ana::NuTruthVar   ppfxXSecWeightsST = kPPFXFluxCVWgtST * kXSecCVWgt2018_smallerDISScale_NT;
const ana::Var          ppfxXSecWeights   = ana::VarFromNuTruthVar(ppfxXSecWeightsST, 1);
const vector<pair<const string, const ana::Var>> weights{
  {"unwgt",             kUnweighted},
  {"ppfxwgt",           kPPFXFluxCVWgt},
  {"ppfxxsecwgt",       ppfxXSecWeights}
};

void makeRealDataFluxes(std::string outputfile = "fRealDataFluxes.root"){

  std::string data_def      = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
  std::string nominal_def   = "def_snapshot nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata";

  SpectrumLoader * dataloader = new SpectrumLoader(data_def);
  SpectrumLoader * nomloader  = new SpectrumLoader(nominal_def);

  std::map<std::string, std::map<std::string, Spectrum*> > dataSpectrumsReco;
  ////////////// Fake data /////////////////////////////////////////////////
  for(std::pair<std::string, ana::HistAxis> unfoldingVar : unfoldingVars){
    for(std::pair<const string, const ana::Var> weight : weights){
      ana::HistAxis& recoAxis = unfoldingVar.second;
      Spectrum * dataReco = new Spectrum(*dataloader, recoAxis, recoCuts.at(unfoldingVar.first), kNoShift, weight.second);
      dataSpectrumsReco[unfoldingVar.first][weight.first] = dataReco;
    }
  }

  ////////////// Fluxes ////////////////////////////////////////////////////
  const int wanted_pdg = 14;
  const Binning FluxBinning  = Binning::Simple(480, 0.0, 120.0);

  // Vertex default to the definitions in NumuCCIncCuts.h
  TVector3 flux_vtxmin = vtxmin;
  TVector3 flux_vtxmax = vtxmax;

  Spectrum * fluxData  = DeriveFlux(*dataloader, FluxBinning, wanted_pdg, &flux_vtxmin, &flux_vtxmax, kNoShift, kPPFXFluxCVWgtST);
  Spectrum * fluxCV    = DeriveFlux(*nomloader, FluxBinning, wanted_pdg, &flux_vtxmin, &flux_vtxmax, kNoShift, kPPFXFluxCVWgtST);

  ////////////// Go() //////////////////////////////////////////////////////
  dataloader->Go();
  nomloader->Go();

  ////////////// Saving ////////////////////////////////////////////////////
  TFile * fOut = new TFile(outputfile.c_str(), "RECREATE");
  TDirectory * dir;

  for (std::pair<std::string, std::map<std::string, Spectrum*> > dataReco : dataSpectrumsReco){
    // Data Spectrums
    dir = fOut->mkdir(dataReco.first.c_str());
    dir = dir->mkdir("reco");
    dir->cd();

    std::map<std::string, Spectrum*>& weightSpecs = dataReco.second;
    for(std::pair<std::string, Spectrum*> weightSpec : weightSpecs)
      weightSpec.second->SaveTo(dir->mkdir(weightSpec.first.c_str()));
  }

  // Flux
  dir = fOut->mkdir("fluxes");
  dir->cd();
  fluxData->SaveTo(dir->mkdir("Data"));
  fluxCV->SaveTo(dir->mkdir("CV"));
}