specprod_numuccinc.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////
2 /// \brief Macro for numucc_inc for running on NERSC
3 ///
4 /// \author Connor Johnson
5 /// \date Sept 2019
6 //////////////////////////////////////////////////////////////////
7 
9 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Vars/XsecTunes.h"
14 
17 
23 
24 #include "TFile.h"
25 
26 #include <iostream>
27 #include <string>
28 #include <map>
29 
30 const std::string NOMINAL_DATASET = "nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata";
31 const unsigned int N_GENIE_UNIVERSES = 1000;
32 const unsigned int N_PPFX_UNIVERSES = 100;
33 
35 const ana::Var wgt_xs = ana::VarFromNuTruthVar(wgt_xsST, 1);
38 
39 void specprod_numuccinc(std::string outfileName, int firstGenie, int lastGenie, int firstPPFX = -1, int lastPPFX = -1)
40 {
41  // Maps for the UnfoldingVariable class. Store in order of [(uint) iuniv][(string) "varname"]
42  std::map<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > > genieUnfoldingVars;
43  std::map<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > > ppfxUnfoldingVars;
44 
45  // Create the nominal spectrum loader
46  ana::SpectrumLoader * nominalLoader = new ana::SpectrumLoader(NOMINAL_DATASET.c_str());
47  nominalLoader->SetSpillCut(ana::kStandardSpillCuts);
48 
49  if ((firstGenie < 0 || lastGenie < firstGenie) && (firstPPFX < 0 || lastPPFX < firstPPFX))
50  {
51  cout << "Invalid choices for both genie and ppfx universes. Will not create spectrums.\nExiting..." << endl;
52  exit(1);
53  }
54 
55  // Load the genie, use firstGenie < 0 to skip them all
56  if (firstGenie >= 0 && lastGenie >= firstGenie)
57  {
58  cout << "Running " << to_string(lastGenie - firstGenie + 1) << " genie universes. (" << to_string(firstGenie) << "-" << to_string(lastGenie) << ")" << endl;
60  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
61 
62  for(int igenie = firstGenie; igenie <= lastGenie; ++igenie)
63  {
64  genieUnfoldingVars[igenie]["EAvail"] = new ana::xsec::UnfoldingVariable(nominalLoader,
66  &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, &genie_shifts[igenie], &wgtST);
67 
68  genieUnfoldingVars[igenie]["ENu"] = new ana::xsec::UnfoldingVariable(nominalLoader,
70  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &genie_shifts[igenie], &wgtST);
71 
72  genieUnfoldingVars[igenie]["Q2"] = new ana::xsec::UnfoldingVariable(nominalLoader,
74  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, &genie_shifts[igenie], &wgtST);
75  }
76  }
77 
78  // Load the ppfx, use firstPPFX < 0 to skip them all
79  if (firstPPFX >= 0 && lastPPFX >= firstPPFX)
80  {
81  cout << "Running " << to_string(lastPPFX - firstPPFX + 1) << " ppfx universes. (" << to_string(firstPPFX) << "-" << to_string(lastPPFX) << ")" << endl;
82  const std::vector<ana::NuTruthVar> vppfxST = ana::GetkPPFXFluxUnivWgtST();
83 
84  for(int ippfx = firstPPFX; ippfx <= lastPPFX; ++ippfx)
85  {
86  const ana::NuTruthVar * xsTimesPPFXi = new ana::NuTruthVar(wgt_xsST * vppfxST[ippfx]);
87 
88  ppfxUnfoldingVars[ippfx]["EAvail"] = new ana::xsec::UnfoldingVariable(nominalLoader,
91 
92  ppfxUnfoldingVars[ippfx]["ENu"] = new ana::xsec::UnfoldingVariable(nominalLoader,
95 
96  ppfxUnfoldingVars[ippfx]["Q2"] = new ana::xsec::UnfoldingVariable(nominalLoader,
99  }
100  }
101 
102  // Fill the spectrums
103  nominalLoader->Go();
104 
105  // Store it all.
106  TFile* fOut = new TFile(outfileName.c_str(), "RECREATE");
107  fOut->cd();
108  TDirectory* dir;
109 
110  // Save the genie universes (carefully grabbing by references for slight optimization)
111  for (const std::pair<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > >& genieUnivVars : genieUnfoldingVars)
112  {
113  unsigned int genieUniv = genieUnivVars.first;
114  std::string genieUnivStr = to_string(genieUniv);
115  std::string folderName = "genie" + genieUnivStr;
116  const std::map<std::string, ana::xsec::UnfoldingVariable*>& unfoldVars = genieUnivVars.second;
117 
118  dir = fOut->mkdir(folderName.c_str());
119  dir->cd();
120  for (const std::pair<std::string, ana::xsec::UnfoldingVariable*>& genieUnivVar : unfoldVars)
121  {
122  std::string varName = genieUnivVar.first;
123  ana::xsec::UnfoldingVariable * var = genieUnivVar.second;
124 
125  TDirectory * varDir = dir->mkdir(varName.c_str());
126  var->SaveSpectrums(varDir);
127  }
128  }
129 
130  // Save the ppfx universes
131  for (const std::pair<unsigned int, std::map<std::string, ana::xsec::UnfoldingVariable* > >& ppfxUnivVars : ppfxUnfoldingVars)
132  {
133  unsigned int ppfxUniv = ppfxUnivVars.first;
134  std::string ppfxUnivStr = to_string(ppfxUniv);
135  std::string folderName = "ppfx" + ppfxUnivStr;
136  const std::map<std::string, ana::xsec::UnfoldingVariable*>& unfoldVars = ppfxUnivVars.second;
137 
138  dir = fOut->mkdir(folderName.c_str());
139  dir->cd();
140  for (const std::pair<std::string, ana::xsec::UnfoldingVariable*>& ppfxUnivVar : unfoldVars)
141  {
142  std::string varName = ppfxUnivVar.first;
143  ana::xsec::UnfoldingVariable * var = ppfxUnivVar.second;
144 
145  TDirectory * varDir = dir->mkdir(varName.c_str());
146  var->SaveSpectrums(varDir);
147  }
148  }
149  fOut->Close();
150  return;
151 }
const NuTruthCut kIsTrueSigST
const HistAxis kRecoQ2StandardAxis("Reco Q2 (GeV)", q2bins, kRecoq2)
void SaveSpectrums(TDirectory *d)
Save each necessary Spectrum to its own subfolder of dir.
const NuTruthCut kIsTrueSig1DST
const ana::Var wgt_xs
void SetSpillCut(const SpillCut &cut)
const NuTruthHistAxis kTrueQ2StandardAxisST("True Q2 (GeV)", q2bins, kTrueQ2_NT)
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
virtual void Go() override
Load all the registered spectra.
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:48
const std::string NOMINAL_DATASET
Macro for numucc_inc for running on NERSC.
const ana::Var wgt
const unsigned int N_PPFX_UNIVERSES
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void specprod_numuccinc(std::string outfileName, int firstGenie, int lastGenie, int firstPPFX=-1, int lastPPFX=-1)
const Cut kAllNumuCC1DCuts
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
const ana::NuTruthVar wgtST
std::vector< SystShifts > GetSystShifts()
exit(0)
const NuTruthHistAxis kTrueMuKEVsCosVsEavailStandardAxisST("True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevseavailbins, kTrueMuKEVsCosVsEavailST)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const HistAxis kRecoMuKEVsCosVsEavailStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavail)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Template for Var and SpillVar.
const HistAxis kRecoEStandardAxis("Reconstructed Neutrino Energy (GeV)", enubins, kRecoE)
const NuTruthHistAxis kTrueEStandardAxisST("True Neutrino Energy (GeV)", enubins, kTrueEST)
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 unsigned int N_GENIE_UNIVERSES
const Cut kAllNumuCCCuts
const ana::NuTruthVar wgt_xsST