specprod_systematics.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////
2 /// \brief Macro to create a given xsec systematic for the
3 /// numucc inc analysis
4 /// \author Connor Johnson
5 /// \date June 2019
6 //////////////////////////////////////////////////////////////////
7 
8 #include "CAFAna/Vars/Vars.h"
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Systs/Systs.h"
18 #include "CAFAna/Systs/BeamSysts.h"
19 #include "CAFAna/Systs/EnergySysts2018.h"
21 
24 
31 
32 #include "TFile.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 
36 #include <iostream>
37 #include <cmath>
38 #include <sstream>
39 #include <string>
40 #include <fstream>
41 #include <iomanip>
42 #include <map>
43 #include <algorithm>
44 
45 //////////////////////////////////////////////////////////////////
46 /// List of systematics to allow studies over. Each section treated differently
47 //////////////////////////////////////////////////////////////////
48 
49 /// Systematics defined by a separate hist axis
50 std::map<string, const vector<const ana::HistAxis* > > histaxis_systematics = {
55 };
56 
57 /// Systematics defined by a BeamSyst
58 std::map<string, const ana::BeamSyst*> beamSystematics = {
59  {"2kA", &ana::kBeamHornCurrent},
60  {"02mmBeamSpotSize", &ana::kBeamSpotSize},
61  {"1mmBeamShiftX", &ana::kBeamPosX},
62  {"1mmBeamShiftY", &ana::kBeamPosY},
63  {"3mmHorn1X", &ana::kBeamH1PosX},
64  {"3mmHorn1Y", &ana::kBeamH1PosY},
65  {"3mmHorn2X", &ana::kBeamH2PosX},
66  {"3mmHorn2Y", &ana::kBeamH2PosY},
67  {"7mmTargetZ", &ana::kBeamTarget},
68  {"MagneticFieldinDecayPipe", &ana::kBeamMagField},
69  {"1mmHornWater", &ana::kBeamGeomWater}
70 };
71 
72 /// Other pre-defined systematics
73 std::map<string, const ana::ISyst*> otherSystematics = {
75 };
76 
77 /// Systematics that run on the nominal dataset
78 std::vector<std::string> standardVars =
79 {
80  "cv",
81  "nominal",
82  "lightdown",
83  "lightup",
84  "ckv",
85  "calibneg",
86  "calibpos",
87  "calibshift",
88  "neutron_nom",
89  "nom_good_seventh"
90 };
91 
92 /// Systematics that use a different dataset
93 std::map<string, const std::string> dataset_to_samdef = {
94  {"nominal", "nd_fhc_remid-hotfix_caf_minus_muonid_training_minus_fakedata"},
95  {"nom_good_seventh", "nd_fhc_remid-hotfix_goodseventh_common_subruns"},
96  {"lightdown", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1_good_seventh"},
97  {"lightup", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1_good_seventh"},
98  {"ckv", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1"},
99  {"calibneg", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1_good_seventh"},
100  {"calibpos", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1_good_seventh"},
101  {"calibshift", "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1_good_seventh"},
102  {"neutron", "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1"},
103  {"neutron_nom", "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1"}
104 };
105 
106 //////////////////////////////////////////////////////////////////
107 /// Override folder naming in certain cases
108 //////////////////////////////////////////////////////////////////
109 std::map<string, const std::string> syst_to_folder = {
110  {"muoneup", "MuESUp"},
111  {"muonedw", "MuESDw"}
112 };
113 
114 //////////////////////////////////////////////////////////////////
115 /// Define the weights. Build off of each other, and default to "1".
116 //////////////////////////////////////////////////////////////////
117 
122 
123 //////////////////////////////////////////////////////////////////
124 /// Process all \a systs input, and store output in \a outfile
125 //////////////////////////////////////////////////////////////////
126 void specprod_systematics(const std::vector<std::string> systs, const std::string outfile = "fSystSpec.root"){
127  std::map<std::string, ana::SpectrumLoaderBase*> vloaders; ///< Map of loaders
128  std::map<std::string, std::map<std::string, ana::xsec::UnfoldingVariable*>> vUnfoldingVars; ///< Nested map of UnfoldingVariables ([systName][variableName])
129 
130  // Loop through systematics
131  for (unsigned int isyst = 0; isyst < systs.size(); isyst++){
132  std::string syst = systs[isyst];
134  if (dataset_to_samdef.find(syst) != dataset_to_samdef.end())
135  dataset = syst;
136  else
137  dataset = "nominal";
138 
139  cout << "Creating spectrums for spec: " << syst << " on dataset: " << dataset << endl;
140 
141  /// Create loader, only one per dataset
143  if (vloaders.find(dataset) != vloaders.end())
144  loader = vloaders[dataset];
145  else{
146  loader = new ana::SpectrumLoader(dataset_to_samdef[dataset]);
148  vloaders[dataset] = loader;
149  }
150 
151  // // CV, and the systematic datasets
152  if (std::find(standardVars.begin(), standardVars.end(), syst) != standardVars.end()){
153  vUnfoldingVars[syst]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
156 
157  vUnfoldingVars[syst]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
160 
161  vUnfoldingVars[syst]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
164  }
165 
166  // HistAxis systematics (muon energy, angle systematics)
167  else if (histaxis_systematics.find(syst) != histaxis_systematics.end()){
168  // const ana::HistAxis * eavail_histaxis = histaxis_systematics[syst][0];
169  const ana::HistAxis * enu_histaxis = histaxis_systematics[syst][1];
170  const ana::HistAxis * q2_histaxis = histaxis_systematics[syst][2];
171  vUnfoldingVars[syst]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
174 
175  vUnfoldingVars[syst]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
176  &ana::kTrueEStandardAxisST, enu_histaxis,
178 
179  vUnfoldingVars[syst]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
180  &ana::kTrueQ2StandardAxisST, q2_histaxis,
182  }
183 
184  // Beam Systematics
185  else if (beamSystematics.find(syst) != beamSystematics.end()){
186  std::string bs_up = syst + "_Up";
187  std::string bs_dw = syst + "_Dw";
188 
189  // First the up-shifts
190  const ana::BeamSyst* bshift = beamSystematics[syst];
191  ana::SystShifts * bsystshiftup = new ana::SystShifts(bshift, +1);
192  vUnfoldingVars[bs_up]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
194  &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, bsystshiftup, &wgtST);
195  vUnfoldingVars[bs_up]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
197  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftup, &wgtST);
198  vUnfoldingVars[bs_up]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
200  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftup, &wgtST);
201 
202  // Next the down-shifts
203  ana::SystShifts * bsystshiftdw = new ana::SystShifts(bshift, -1);
204  vUnfoldingVars[bs_dw]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
206  &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, bsystshiftdw, &wgtST);
207  vUnfoldingVars[bs_dw]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
209  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftdw, &wgtST);
210  vUnfoldingVars[bs_dw]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
212  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, bsystshiftdw, &wgtST);
213  }
214  else if (otherSystematics.find(syst) != otherSystematics.end()){
215  std::string syst_up = syst + "_Up";
216  std::string syst_dw = syst + "_Dw";
217 
218  // First the up-shifts
219  const ana::ISyst* isyst = otherSystematics[syst];
220  ana::SystShifts * systshiftup = new ana::SystShifts(isyst, +1);
221  vUnfoldingVars[syst_up]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
223  &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, systshiftup, &wgtST);
224  vUnfoldingVars[syst_up]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
226  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftup, &wgtST);
227  vUnfoldingVars[syst_up]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
229  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftup, &wgtST);
230 
231  // Next the down-shifts
232  ana::SystShifts * systshiftdw = new ana::SystShifts(isyst, -1);
233  vUnfoldingVars[syst_dw]["EAvail"] = new ana::xsec::UnfoldingVariable(loader,
235  &ana::kAllNumuCCCuts, &ana::kIsTrueSigST, systshiftdw, &wgtST);
236  vUnfoldingVars[syst_dw]["ENu"] = new ana::xsec::UnfoldingVariable(loader,
238  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftdw, &wgtST);
239  vUnfoldingVars[syst_dw]["Q2"] = new ana::xsec::UnfoldingVariable(loader,
241  &ana::kAllNumuCC1DCuts, &ana::kIsTrueSig1DST, systshiftdw, &wgtST);
242  }
243  else {
244  cerr << "Attempt to run on unknown systematic. Syst: " << syst << "Not found. Skipping." << endl;
245  exit(1);
246  }
247  } // end systematic loop
248 
249  // Call go on all loaders (no particular order, could be dumb. Hmm...)
250  for(pair<std::string, ana::SpectrumLoaderBase*> loader : vloaders){
251  loader.second->Go();
252  }
253 
254  // Store it all.
255  TFile* fOut = new TFile(outfile.c_str(), "RECREATE");
256  fOut->cd();
257  TDirectory* dir;
258 
259  // Loop through all systematics
260  for (std::pair<std::string, std::map<std::string, ana::xsec::UnfoldingVariable*> > systVars : vUnfoldingVars){
261  std::string syst = systVars.first;
262  std::map<std::string, ana::xsec::UnfoldingVariable*> unfoldVars = systVars.second;
263 
264  std::string foldername = "";
265  if (syst_to_folder.find(syst) != syst_to_folder.end())
266  foldername = syst_to_folder[syst];
267  else
268  foldername = syst;
269 
270  dir = fOut->mkdir(syst.c_str());
271  dir->cd();
272  for (std::pair<std::string, ana::xsec::UnfoldingVariable*> unfoldVar : unfoldVars){
273  std::string varName = unfoldVar.first;
274  ana::xsec::UnfoldingVariable * var = unfoldVar.second;
275 
276  TDirectory * varDir = dir->mkdir(varName.c_str());
277  var->SaveSpectrums(varDir);
278  }
279  }
280  fOut->Close();
281  return;
282 }
283 
284 //////////////////////////////////////////////////////////////////
285 /// Process a single \a syst and store it in \a outfile
286 //////////////////////////////////////////////////////////////////
288  if (outfile.empty())
289  outfile = "fSystSpec_" + syst + ".root";
290  if (syst == "all")
292  "cv", "muoneup", "muonedw", "angleup", "angledw", "2kA",
293  "02mmBeamSpotSize",
294  "1mmBeamShiftX",
295  "1mmBeamShiftY",
296  "3mmHorn1X",
297  "3mmHorn1Y",
298  "3mmHorn2X",
299  "3mmHorn2Y",
300  "7mmTargetZ",
301  "MagneticFieldinDecayPipe",
302  "1mmHornWater",
303  "neutron",
304  "lightdown", "lightup", "ckv",
305  "calibneg", "calibpos", "calibshift",
306  "neutron_nom", "nom_good_seventh"
307  }, outfile);
308  else
309  specprod_systematics(std::vector<string> {syst}, outfile);
310  return;
311 }
const HistAxis kRecoEStandardAxis_MuonEUp("Reconstructed Neutrino Energy (GeV)", enubins, kRecoEUp)
const NuTruthCut kIsTrueSigST
const HistAxis kRecoQ2StandardAxis("Reco Q2 (GeV)", q2bins, kRecoq2)
const HistAxis kRecoQ2StandardAxis_MuonEDw("Reco Q2 (GeV)", q2bins, kRecoq2_MuonEDw)
void SaveSpectrums(TDirectory *d)
Save each necessary Spectrum to its own subfolder of dir.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const BeamSyst kBeamTarget((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"7mmTargetZ","+/- 7mm Target z Position")
Target z position shift +/-7mm.
Definition: BeamSysts.h:113
Beam systematics:
Definition: BeamSysts.h:37
const NuTruthCut kIsTrueSig1DST
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const BeamSyst kBeamHornCurrent((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"2kA","2kAHornCurrent","+/- 2kA Horn Current")
Horn Current +/-2kA.
Definition: BeamSysts.h:97
OStream cerr
Definition: OStream.cxx:7
const ana::NuTruthVar wgtST
void SetSpillCut(const SpillCut &cut)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const HistAxis kRecoQ2StandardAxis_AngleDw("Reco Q2 (GeV)", q2bins, kRecoq2_AngleDw)
const BeamSyst kBeamGeomWater((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmHornWater","+/- 1mm water on Horn 1")
Water layer on horn 1: +/- 1mm.
Definition: BeamSysts.h:128
std::map< string, const std::string > syst_to_folder
Override folder naming in certain cases.
std::map< string, const std::string > dataset_to_samdef
Systematics that use a different dataset.
const HistAxis kRecoQ2StandardAxis_AngleUp("Reco Q2 (GeV)", q2bins, kRecoq2_AngleUp)
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleDw("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailAngleDw)
const NuTruthHistAxis kTrueQ2StandardAxisST("True Q2 (GeV)", q2bins, kTrueQ2_NT)
const BeamSyst kBeamH1PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1X","+/- 3mm Horn 1 X Position")
Horn 1 and 2 position +-3mm in X and Y separately.
Definition: BeamSysts.h:107
const ana::NuTruthVar wgt_xsST
Define the weights. Build off of each other, and default to "1".
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisUp("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailUp)
std::map< string, const ana::BeamSyst * > beamSystematics
Systematics defined by a BeamSyst.
const HistAxis kRecoEStandardAxis_MuonEDw("Reconstructed Neutrino Energy (GeV)", enubins, kRecoEDw)
std::map< string, const vector< const ana::HistAxis * > > histaxis_systematics
Macro to create a given xsec systematic for the numucc inc analysis.
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:48
loader
Definition: demo0.py:10
const ana::Var wgt
std::map< string, const ana::ISyst * > otherSystematics
Other pre-defined systematics.
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const BeamSyst kBeamPosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftY","Beam Position Y")
Definition: BeamSysts.h:104
Base class for the various types of spectrum loader.
const Cut kAllNumuCC1DCuts
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
const BeamSyst kBeamH1PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1Y","+/- 3mm Horn 1 Y Position")
Definition: BeamSysts.h:108
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
const BeamSyst kBeamSpotSize((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"0.2mmBeamSpotSize","0p2mmBeamSpotSize"," 1.3 +/- 0.2 mm Spot Size")
Beam Spot Size 1.3 +/- 0.2 mm both XY.
Definition: BeamSysts.h:100
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisAngleUp("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailAngleUp)
const BeamSyst kBeamH2PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2X","+/- 3mm Horn 2 X Position")
Definition: BeamSysts.h:109
const ana::Var wgt_xs
exit(0)
void specprod_systematics(const std::vector< std::string > systs, const std::string outfile="fSystSpec.root")
Process all systs input, and store output in outfile.
const BeamSyst kBeamPosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftX","Beam Position X")
Beam position on target +-1 mm, X/Y separately.
Definition: BeamSysts.h:103
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
std::vector< std::string > standardVars
Systematics that run on the nominal dataset.
const HistAxis kRecoMuKEVsCosVsEavailStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavail)
const BeamSyst kBeamH2PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2Y","+/- 3mm Horn 2 Y Position")
Definition: BeamSysts.h:110
const BeamSyst kBeamMagField((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"Magnetic Field in Decay Pipe","MagnFieldDecayPipe","Magnetic Field in Decay Pipe")
Constant magnetic field in decay pipe.
Definition: BeamSysts.h:116
Template for Var and SpillVar.
const HistAxis kRecoMuKEVsCosVsEavailStandardAxisDw("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevseavailbins, kRecoMuKEVsCosVsEavailDw)
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 Cut kAllNumuCCCuts
FILE * outfile
Definition: dump_event.C:13
enum BeamMode string
const HistAxis kRecoQ2StandardAxis_MuonEUp("Reco Q2 (GeV)", q2bins, kRecoq2_MuonEUp)