makeFakeDataFluxes.C
Go to the documentation of this file.
1 // makeFakeData.C
2 
3 #include "CAFAna/Vars/Vars.h"
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Systs/Systs.h"
14 
17 
23 
24 #include "TFile.h"
25 #include "TH1.h"
26 #include "TH2.h"
27 
28 #include <iostream>
29 #include <cmath>
30 #include <sstream>
31 #include <string>
32 #include <fstream>
33 #include <iomanip>
34 #include <map>
35 
36 using namespace ana;
37 
38 std::map<std::string, std::pair<ana::HistAxis, ana::HistAxis>> unfoldingVars{
40  {"fakedata_enu", {kRecoEStandardAxis, kTrueEStandardAxis}},
41  {"fakedata_q2", {kRecoQ2StandardAxis, kTrueQ2StandardAxis}}
42 };
43 
44 std::map<std::string, ana::Cut> fakeDataCuts{
45  {"fakedata", kAllNumuCCCuts},
46  {"fakedata_enu", kAllNumuCC1DCuts},
47  {"fakedata_q2", kAllNumuCC1DCuts}
48 };
49 
51 const ana::Var ppfxXSecWeights = ana::VarFromNuTruthVar(ppfxXSecWeightsST, 1);
52 const vector<pair<const string, const ana::Var>> weights{
53  {"unwgt", kUnweighted},
54  {"ppfxwgt", kPPFXFluxCVWgt},
55  {"ppfxxsecwgt", ppfxXSecWeights}
56 };
57 
58 void makeFakeDataFluxes(std::string outputfile = "fSpecsFakeData.root"){
59 
60  std::string fakedata_def = "dataset_def_name_newest_snapshot nd_fhc_remid-hotfix_caf_fake_data";
61  std::string nominal_def = "dataset_def_name_newest_snapshot nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata";
62 
63  SpectrumLoaderBase * fdloader = new SpectrumLoader(fakedata_def);
64  SpectrumLoaderBase * nomloader = new SpectrumLoader(nominal_def);
65 
66  std::map<std::string, std::map<std::string, Spectrum*> > fakeDataSpectrumsReco;
67  std::map<std::string, std::map<std::string, Spectrum*> > fakeDataSpectrumsTrue;
68  ////////////// Fake data /////////////////////////////////////////////////
69  for(std::pair<std::string, std::pair<ana::HistAxis, ana::HistAxis>> unfoldingVar : unfoldingVars){
70  for(std::pair<const string, const ana::Var> weight : weights){
71  ana::HistAxis& recoAxis = unfoldingVar.second.first;
72  ana::HistAxis& trueAxis = unfoldingVar.second.second;
73  Spectrum * fakeDataReco = new Spectrum(*fdloader, recoAxis, fakeDataCuts.at(unfoldingVar.first), kNoShift, weight.second);
74  Spectrum * fakeDataTrue = new Spectrum(*fdloader, trueAxis, fakeDataCuts.at(unfoldingVar.first) && kIsTrueSig, kNoShift, weight.second);
75  fakeDataSpectrumsReco[unfoldingVar.first][weight.first] = fakeDataReco;
76  fakeDataSpectrumsTrue[unfoldingVar.first][weight.first] = fakeDataTrue;
77  }
78  }
79 
80  ////////////// Fluxes ////////////////////////////////////////////////////
81  const int wanted_pdg = 14;
82  const Binning FluxBinning = Binning::Simple(480, 0.0, 120.0);
83 
84  // Vertex default to the definitions in NumuCCIncCuts.h
85  TVector3 flux_vtxmin = vtxmin;
86  TVector3 flux_vtxmax = vtxmax;
87 
88  Spectrum * fluxFD = DeriveFlux(*fdloader, FluxBinning, wanted_pdg, &flux_vtxmin, &flux_vtxmax, kNoShift, kPPFXFluxCVWgtST);
89  Spectrum * fluxCV = DeriveFlux(*nomloader, FluxBinning, wanted_pdg, &flux_vtxmin, &flux_vtxmax, kNoShift, kPPFXFluxCVWgtST);
90 
91  ////////////// Go() //////////////////////////////////////////////////////
92  fdloader->Go();
93  nomloader->Go();
94 
95  ////////////// Saving ////////////////////////////////////////////////////
96  TFile * fOut = new TFile(outputfile.c_str(), "RECREATE");
97  TDirectory * dir;
98 
99  for (std::pair<std::string, std::map<std::string, Spectrum*> > fakeDataReco : fakeDataSpectrumsReco){
100  // Fake Data
101  dir = fOut->mkdir(fakeDataReco.first.c_str());
102  dir = dir->mkdir("reco");
103  dir->cd();
104 
105  std::map<std::string, Spectrum*>& weightSpecs = fakeDataReco.second;
106  for(std::pair<std::string, Spectrum*> weightSpec : weightSpecs)
107  weightSpec.second->SaveTo(dir->mkdir(weightSpec.first.c_str()));
108  }
109 
110  for (std::pair<std::string, std::map<std::string, Spectrum*> > fakeDataTrue : fakeDataSpectrumsTrue){
111  // Fake Data
112  dir = fOut->GetDirectory(fakeDataTrue.first.c_str(), true);
113  dir = dir->mkdir("true");
114  dir->cd();
115 
116  std::map<std::string, Spectrum*>& weightSpecs = fakeDataTrue.second;
117  for(std::pair<std::string, Spectrum*> weightSpec : weightSpecs)
118  weightSpec.second->SaveTo(dir->mkdir(weightSpec.first.c_str()));
119  }
120 
121  // Flux
122  dir = fOut->mkdir("fluxes");
123  dir->cd();
124  fluxFD->SaveTo(dir->mkdir("FakeData"));
125  fluxCV->SaveTo(dir->mkdir("CV"));
126 }
void makeFakeDataFluxes(std::string outputfile="fSpecsFakeData.root")
const ana::Var ppfxXSecWeights
const HistAxis kTrueEStandardAxis
const HistAxis kRecoQ2StandardAxis("Reco Q2 (GeV)", q2bins, kRecoq2)
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
const TVector3 vtxmin(-130,-176, 225)
std::vector< double > Spectrum
Definition: Constants.h:746
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
virtual void Go()=0
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:47
const HistAxis kTrueMuKEVsCosVsEavailStandardAxis
const SystShifts kNoShift
Definition: SystShifts.cxx:22
Base class for the various types of spectrum loader.
const Cut kIsTrueSig
const Cut kAllNumuCC1DCuts
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:8
const ana::NuTruthVar ppfxXSecWeightsST
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:58
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
std::map< std::string, ana::Cut > fakeDataCuts
const HistAxis kRecoMuKEVsCosVsEavailStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavail)
const TVector3 vtxmax(160, 160, 1000)
Template for Vars applied to any type of object.
Definition: FwdDeclare.h:12
const HistAxis kTrueQ2StandardAxis
const HistAxis kRecoEStandardAxis("Reconstructed Neutrino Energy (GeV)", enubins, kRecoE)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
Definition: Var.h:14
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
const Cut kAllNumuCCCuts
std::map< std::string, std::pair< ana::HistAxis, ana::HistAxis > > unfoldingVars
const vector< pair< const string, const ana::Var > > weights
enum BeamMode string