Nus18SystsBeamTranspLoad.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 // Beam Transport systematics for Nus18
3 //
4 // Author: Sijith Edayath
5 // Modified by Mike Wallbank (wallbank@fnal.gov)
6 /////////////////////////////////////////////////////////////////////
7 
9 #include "CAFAna/Core/Loaders.h"
10 #include "CAFAna/Core/SystShifts.h"
11 #include "CAFAna/Core/Utilities.h"
12 #include "CAFAna/Cuts/Cuts.h"
13 #include "NuXAna/Cuts/NusCuts18.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
23 #include "CAFAna/Systs/BeamSysts.h"
24 #include "CAFAna/Systs/Systs.h"
26 #include "CAFAna/Vars/HistAxes.h"
27 #include "NuXAna/Vars/HistAxes.h"
29 
32 
33 using namespace ana;
34 
35 /// Get Beam Transport systematics labels
36 std::map<std::string, std::vector<std::string> > GetBeamTranspShiftLabels();
37 
38 /// Run systematics for Beam Transport
40 
41  TH1::AddDirectory(0);
42  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;
43 
44  // Get options
45  ECAFType cafType;
46  Loaders::FluxType fluxType;
47  bool extrap = false;
48  if (opt.Contains("caf", TString::kIgnoreCase))
49  cafType = ECAFType::kFullCAF;
50  else if (opt.Contains("concat", TString::kIgnoreCase))
51  cafType = ECAFType::kNusConcat;
52  else {
53  std::cout << "ECAFType required: caf or concat" << std::endl;
54  abort();
55  }
56  if (opt.Contains("fhc", TString::kIgnoreCase))
57  fluxType = Loaders::FluxType::kFHC;
58  else if (opt.Contains("rhc", TString::kIgnoreCase))
59  fluxType = Loaders::FluxType::kRHC;
60  else {
61  std::cout << "FluxType required: fhc or rhc" << std::endl;
62  abort();
63  }
64  if (opt.Contains("extrap", TString::kIgnoreCase))
65  extrap = true;
66 
67  // Set up the loaders
68  Loaders* loaders = new Prod4NomLoaders(cafType, fluxType);
70 
71  // Vector of sigmas, levels to set each systematic
72  const std::vector<int> sigmas{-3, -2, -1, +1, +2, +3};
73 
74  // Get the labels for all Beam Transport systematics
75  std::map<std::string, std::vector<std::string> > shiftLabels = GetBeamTranspShiftLabels();
76 
77  // Define map of systematic objects, indexed by same labels as previous map
78  std::map<std::string, std::map<int, SystShifts*> > systematics;
79 
80  // Make SystShifts
81  std::string rootfile = FindCAFAnaDir() + "/data/beam/TABeamSyst.root";
82  for (const auto& shiftLabel : shiftLabels)
83  for(const auto& sigma : sigmas)
84  systematics[shiftLabel.first][sigma] =
85  new SystShifts(new BeamSyst(rootfile, shiftLabel.second[0], shiftLabel.second[1]), sigma);
86 
87  // Define a map of samples and cuts/axes to set them up
88  // Each sample, essentially, is a way to set up the analysis
89  HistAxis AxisRecoE = kNus18AxisE;
90  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
91  cut_samples["Nus18"] = std::make_pair(&AxisRecoE, std::make_pair(new Cut(kNus18ND), new Cut(kNus18FD)));
92 
93  // Set up maps to the decompositions/predictions (each flavor component)
94  // Nominal maps are indexed purely by sample label
95  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
96  std::map<std::string, NDPredictionSterile*> predsND_nominal;
97  std::map<std::string, FDPredictionSterile*> predsFD_nominal;
98  std::map<std::string, PredictionSterile*> predsExtrap_nominal;
99  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted;
100  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted;
101  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsExtrap_shifted;
102 
103  // Set up the actual decompositions/predictions
104  // sample.first = the sample label
105  // shiftLabel.first = the systematic shift label
106  // sample.second.first = the sample HistAxis
107  // sample.second.second = the sample Cut
108  // syst.first = the integer number of sigmas for the shift
109  // syst.second = the systematic being shifted
110  for (const auto& sample : cut_samples) {
111 
112  const NDPredictionGenerator gennd(*sample.second.first, *sample.second.second.first, kNoShift, kReweight);
113  const FDPredictionGenerator genfd(*sample.second.first, *sample.second.second.second, kNoShift, kReweight);
114  SterileGenerator* genExtrap;
115  if (extrap)
116  genExtrap = new SterileGenerator(*sample.second.first, kNus18BinsNumuCCAxis,
117  *sample.second.second.second, *sample.second.second.first,
118  kNumuCutND2018, kNoShift, kReweight);
119 
120  // Set up the nominal decompositions
121  // Only one is needed per sample
122  predsND_nominal[sample.first] = (NDPredictionSterile*)gennd.Generate(*loaders).release();
123  predsFD_nominal[sample.first] = (FDPredictionSterile*)genfd.Generate(*loaders).release();
124  if (extrap)
125  predsExtrap_nominal[sample.first] = (PredictionSterile*)genExtrap->Generate(*loaders).release();
126 
127  // Set up the shifted decompositions
128  // One is needed for every sample, systematic, and sigma level
129  for(const auto& shiftLabel : shiftLabels) {
130  for(const auto& syst : systematics[shiftLabel.first]) {
131  predsND_shifted[sample.first][shiftLabel.first][syst.first]
132  = (NDPredictionSterile*)gennd.Generate(*loaders, *syst.second).release();
133  predsFD_shifted[sample.first][shiftLabel.first][syst.first]
134  = (FDPredictionSterile*)genfd.Generate(*loaders, *syst.second).release();
135  if (extrap)
136  predsExtrap_shifted[sample.first][shiftLabel.first][syst.first]
137  = (PredictionSterile*)genExtrap->Generate(*loaders, *syst.second).release();
138  }
139  }
140 
141  if (extrap)
142  delete genExtrap;
143 
144  }
145 
146  loaders->Go();
147 
148  std::string folder = "./";
149  std::string filenm = "SystsBeam18";
150 
151  std::string fullLocation = folder + filenm + ".root";
152  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");
153  SaveMaps(rootF, predsND_nominal, predsND_shifted);
154  SaveMaps(rootF, predsFD_nominal, predsFD_shifted);
155  if (extrap)
156  SaveMaps(rootF, predsExtrap_nominal, predsExtrap_shifted);
157  rootF->Close();
158 
159  // delete loaders;
160 
161  return;
162 
163 }
164 
165 std::map<std::string, std::vector<std::string> > GetBeamTranspShiftLabels() {
166 
167  // Map to set up systematics
168  std::map<std::string, std::vector<std::string> > shiftLabels;
169 
170  // First element is a label for the systematic
171  shiftLabels["BeamHornCurrent"] = std::vector<std::string>();
172  shiftLabels["BeamSpotSize"] = std::vector<std::string>();
173  shiftLabels["BeamposX"] = std::vector<std::string>();
174  shiftLabels["BeamposY"] = std::vector<std::string>();
175  shiftLabels["BeamH1PosX"] = std::vector<std::string>();
176  shiftLabels["BeamH1PosY"] = std::vector<std::string>();
177  shiftLabels["BeamH2PosX"] = std::vector<std::string>();
178  shiftLabels["BeamH2PosY"] = std::vector<std::string>();
179  shiftLabels["TargetPosZ"] = std::vector<std::string>();
180  shiftLabels["BeamExp"] = std::vector<std::string>();
181  shiftLabels["HornWater"] = std::vector<std::string>();
182  shiftLabels["totErr"] = std::vector<std::string>();
183 
184  // Second element is vector for items required to set up Syst object
185  shiftLabels["BeamHornCurrent"].push_back("2kA");
186  shiftLabels["BeamHornCurrent"].push_back("Horn Current");
187  shiftLabels["BeamSpotSize"] .push_back("0.2mmBeamSpotSize");
188  shiftLabels["BeamSpotSize"] .push_back("Spot Size");
189  shiftLabels["BeamposX"] .push_back("1mmBeamShiftX");
190  shiftLabels["BeamposX"] .push_back("Beam Position X");
191  shiftLabels["BeamposY"] .push_back("1mmBeamShiftY");
192  shiftLabels["BeamposY"] .push_back("Beam Position Y");
193  shiftLabels["BeamH1PosX"] .push_back("3mmHorn1X");
194  shiftLabels["BeamH1PosX"] .push_back("Horn 1 X Position");
195  shiftLabels["BeamH1PosY"] .push_back("3mmHorn1Y");
196  shiftLabels["BeamH1PosY"] .push_back("Horn 1 Y Position");
197  shiftLabels["BeamH2PosX"] .push_back("3mmHorn2X");
198  shiftLabels["BeamH2PosX"] .push_back("Horn 2 X Position");
199  shiftLabels["BeamH2PosY"] .push_back("3mmHorn2Y");
200  shiftLabels["BeamH2PosY"] .push_back("Horn 2 Y Position");
201  shiftLabels["TargetPosZ"] .push_back("7mmTargetZ");
202  shiftLabels["TargetPosZ"] .push_back("Target Z position");
203  shiftLabels["BeamExp"] .push_back("Magnetic Field in Decay Pipe");
204  shiftLabels["BeamExp"] .push_back("Magnet Field in Decay Pipe");
205  shiftLabels["HornWater"] .push_back("1mmHornWater");
206  shiftLabels["HornWater"] .push_back("New Horn Geometry and 1mm water");
207  shiftLabels["totErr"] .push_back("totErr");
208  shiftLabels["totErr"] .push_back("Combined Beam Transport Systematics");
209 
210  return shiftLabels;
211 
212 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Beam systematics:
Definition: BeamSysts.h:37
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
void Nus18SystsBeamTranspLoad(TString opt)
Run systematics for Beam Transport.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Generates extrapolated NC predictions using ProportionalDecomp.
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const Cut kNus18FD
Definition: NusCuts18.h:100
const HistAxis kNus18AxisE("Energy Deposited in Scintillator (GeV)", kNus18EnergyBinning, kNus18Energy)
Axes used in Nus18 analysis by nus group.
Definition: HistAxes.h:16
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
std::map< std::string, std::vector< std::string > > GetBeamTranspShiftLabels()
Get Beam Transport systematics labels.
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
Take the output of an extrapolation and oscillate it as required.
double sigma(TH1F *hist, double percentile)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
bool systematics
Definition: fnexvscaf.C:31
A prediction object compatible with sterile oscillations.
const HistAxis kNus18BinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNus18EnergyBinning, kCCE)
Definition: HistAxes.h:17
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Take the output of an extrapolation and oscillate it as required.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
Generates Near Detector predictions.
ECAFType
Definition: Loaders.h:19
enum BeamMode string