GenieMultiverseSyst.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <boost/multi_index_container.hpp>
4 #include <boost/multi_index/member.hpp>
5 #include <boost/multi_index/ordered_index.hpp>
6 #include <set>
7 #include <string>
8 
9 #include "CAFAna/Core/Cut.h"
10 #include "CAFAna/Core/FwdDeclare.h"
11 #include "CAFAna/Core/HistAxis.h"
12 #include "CAFAna/Core/Var.h"
13 #include "CAFAna/Core/SystShifts.h"
14 #include "TGraphAsymmErrors.h"
15 
16 // multi_index abbreviation
17 using boost::multi_index_container;
18 using namespace boost::multi_index;
19 
20 // Forward declare TH2
21 class TH2;
22 
23 namespace ana
24 {
25  /// forward declarations
26  class SpectrumLoaderBase;
27 
28  /// Obtain GENIE systematic band by randomly varying the number of sigmas
29  /// pulled out from a normal distribution for all the enabled GENIE tunable parameters.
30  /// Originally the json format was chosen as the format of knob configuration file.
31  /// Due to this bug: https://svn.boost.org/trac/boost/ticket/6785,
32  /// turned to pure text file instead.
34  {
35  public:
37  {
42  kBandFromMedian
43  };
47  const HistAxis& axis,
48  const Cut& cut,
49  const SystShifts& shift,
50  const Var& wei = kUnweighted,
51  std::string config_pathname = "knob_config.txt");
52 
53  GenieMultiverseSpectra(unsigned int nuniverses,
54  SpectrumLoaderBase& loader,
55  const HistAxis& axis,
56  const Cut& cut,
57  const SystShifts& shift,
58  const Var& wei,
59  std::vector<const ISyst*> systs);
60 
61  /// Multi-universe constructor for spill truch variables.
62  GenieMultiverseSpectra(unsigned int nuniverses,
63  SpectrumLoaderBase& loader,
64  const NuTruthHistAxis& axis,
65  const NuTruthCut& cut,
66  const SystShifts& shift,
67  const NuTruthVar& wei = kNuTruthUnweighted,
68  std::string config_pathname = "knob_config.txt");
69 
70  /// Multi-universe constructor for spill truch variables.
71  GenieMultiverseSpectra(unsigned int nuniverses,
72  SpectrumLoaderBase& loader,
73  const NuTruthHistAxis& axis,
74  const NuTruthCut& cut,
75  const SystShifts& shift,
76  const NuTruthVar& wei,
77  std::vector<const ISyst*> systs);
78 
79  GenieMultiverseSpectra(unsigned int nuniverses,
80  std::vector<std::unique_ptr<ana::Spectrum>>& spectra,
81  bool fromfile = false);
82 
83  GenieMultiverseSpectra(std::vector<std::unique_ptr<ana::Spectrum>>& spectra);
84 
85 
87 
88  /// return the nominal histogram with the RMS of all other universes as
89  /// error band
90  TH1* ToTH1() const;
91  TH1* ToTH1(double pot) const;
92  TGraphAsymmErrors* ToAreaNormalizedTH1(double pot,
93  int col = -1,
94  int errCol = -1,
95  float headroom = 1.3,
96  bool newaxis = true) const;
97 
98  TH2* ToTH2() const;
99  TH2* ToTH2(double pot) const;
100 
101  /// return the nominal universe
102  const Spectrum* Nominal() const;
103  /// return all universes
104  std::vector<const Spectrum*> AllUniverses() const {return {fSpectra.begin(), fSpectra.end()};}
105  /// return RMS extremes
106  const Spectrum* UpperRMS() const;
107  const Spectrum* LowerRMS() const;
108  const Spectrum* Mean() const;
109  const Spectrum* UpperSigma(BandOptions opt = kBandFromNominal) const;
110  const Spectrum* LowerSigma(BandOptions opt = kBandFromNominal) const;
111 
112  /// Save the Spectrum structure of all universe to a TFile directory.
113  void SaveTo(TDirectory* dir, const std::string& name) const;
114 
115  /// Load the Spectrum structure of all universe from a TFile directory.
116  static std::unique_ptr<GenieMultiverseSpectra> LoadFrom(TDirectory* dir, const std::string& name);
117 
118  /// Get the number of universes in an object of this class.
119  unsigned int GetNUniverses(){return fNUniverses;};
120 
121  /// Calculate error band given 2D histograms
122  static std::pair<TH2*, TH2*> GetSigmaHistograms(std::vector<TH2*>, BandOptions opt = kBandFromNominal);
123 
124  const auto & GetShiftTable() const { return fShiftTable; };
125 
126  protected:
127  ///
128  /// spectra of all universes
129  ///
130  std::vector<Spectrum*> fSpectra;
131 
132  int fUniverseSeed = 1001; ///< Random seed sequence.
133  ///< Need to have each universe with the same physical parameters.
134  ///< Start from 1001, giving way to the PPFX multi-universe.
135 
136  mutable BandOptions fCurBandOpt = kNoBand;
137 
138  ///
139  /// function for building knob configuration table
140  ///
141  void BuildKnobConfigTable(std::string config_pathname);
142 
143  ///
144  /// function to create one universe
145  ///
146  Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
147  const HistAxis& axis,
148  const Cut& cut,
149  const SystShifts& shift = kNoShift,
150  const Var& wei = kUnweighted);
151 
152  ///
153  /// Function to create one universe
154  ///
155  Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
156  const HistAxis& axis,
157  const Cut& cut,
158  std::vector<const ISyst*> systs,
159  const Var& wei = kUnweighted);
160 
161  ///
162  /// Function to create one universe for NuTruth variables.
163  ///
164  Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
165  const NuTruthHistAxis& axis,
166  const NuTruthCut& cut,
167  const SystShifts& shift = kNoShift,
168  const NuTruthVar& wei = kNuTruthUnweighted);
169 
170  ///
171  /// Function to create one universe for NuTruth variables.
172  ///
173  Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
174  const NuTruthHistAxis& axis,
175  const NuTruthCut& cut,
176  std::vector<const ISyst*> systs,
177  const NuTruthVar& wei = kNuTruthUnweighted);
178 
179  private:
180  /// find the boundaries from RMS/sigma of all universes and store results to
181  /// fUpperRMS, fLowerRMS, fUpperSigma, fLowerSigma, etc.
182  void FindBandBoundaries(BandOptions opt = kBandFromNominal) const;
183 
184  /// number of universes to generate
185  unsigned int fNUniverses;
186 
187  bool fLoadFromFile = false;
188  /// modes of valid configuration.
189  std::set<int> fModes = {0, 1, 2};
190 
191  std::vector<std::map<const ISyst*, double>> fShiftTable;
192 
193  mutable Spectrum* fUpperRMS = NULL; ///< spectra for upper rms
194  mutable Spectrum* fLowerRMS = NULL; ///< spectra for lower rms
195  mutable Spectrum* fMean = NULL; ///< spectra for mean
196  mutable Spectrum* fUpperSigma = NULL; ///< spectra for upper sigma
197  mutable Spectrum* fLowerSigma = NULL; ///< spectra for lower sigma
198  mutable bool fBoundariesFound = false;
199 
200  /// Find +/- 1 sigma around the mean for a bin of the Spectrum.
201  /// A sigma is defined as enclosing 34% of area with the mean.
202  static float BinSigma(std::vector<float> events, float nsigmas, float mean = -1);
203 
204  /// declarations for knob configuration table
205  // use boost::multi_index to store a lookup table for knob sampling modes
207  {
208  int id;
210  int mode;
211 
212  knob_sampling_mode(int id_,std::string name_,int mode_):id(id_),name(name_),mode(mode_){}
213  };
214 
215  /* tags for accessing the corresponding indices of employee_set */
216  struct id{};
217  struct name{};
218 
219  typedef multi_index_container<
221  indexed_by<
222  ordered_unique<
223  tag<id>, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,int,id)>,
224  ordered_unique<
225  tag<name>,BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,std::string,name)> >
227  protected:
228  ///
229  /// knob configuration table
230  ///
232 
233  };
234  /// typedef for backward compatibility
236 
238  {
239  public:
240  /// This class generates the parameters of the multi-universes. No more, no less.
242  std::string config_pathname = "knob_config.txt");
243  GenieMultiverseParameters(unsigned int nuniverses,
244  std::vector<const ISyst*> systs);
245  /// Return the tables of knob shift values.
246  std::vector<std::map<const ISyst*, double>> ShiftTables() {return fShiftTables;};
247  std::vector<std::map<const ISyst*, double>> NuTruthShiftTables() {return fNuTruthShiftTables;};
248  std::vector<SystShifts> GetSystShifts() {return fSystShifts;};
249  private:
250  int fUniverseSeed; ///< Random seed sequence.
251  ///< Need to have each universe with the same physical parameters.
252  ///< Start from 1001, giving way to the PPFX multi-universe.
253  //int fUniverseSeedTruth;
254  bool fIsNuTruth = false;
255 
256  std::vector<std::map<const ISyst*, double>> fShiftTables; ///< Container of shift tables for all universes.
257  std::vector<std::map<const ISyst*, double>> fNuTruthShiftTables; ///< Container of neutrino truth shift tables for all universes.
258  std::vector<SystShifts> fSystShifts;
259 
260  void BuildKnobConfigTable(std::string config_pathname); ///< Function for building knob configuration table.
261 
263  {
264  int id;
266  int mode;
267 
268  knob_sampling_mode(int id_,std::string name_,int mode_):id(id_),name(name_),mode(mode_){}
269  };
270  /* tags for accessing the corresponding indices of employee_set */
271  struct id{};
272  struct name{};
273 
274  typedef multi_index_container<
276  indexed_by<
277  ordered_unique<
278  tag<id>, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,int,id)>,
279  ordered_unique<
280  tag<name>,BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,std::string,name)> >
282  /// knob configuration table
284  };
285 
286  ///
287  /// Class for shape-only systematic effect.
288  ///
290  {
291  public:
294  const HistAxis& axis,
295  const Cut& cut,
296  const SystShifts& shift,
297  const Var& wei = kUnweighted,
298  std::string config_pathname = "knob_config.txt");
299 
300  GenieMultiverseNormalizedSpectra(unsigned int nuniverses,
301  SpectrumLoaderBase& loader,
302  const HistAxis& axis,
303  const Cut& cut,
304  const SystShifts& shift,
305  const Var& wei,
306  std::vector<const ISyst*> systs);
307 
309  unsigned int nuniverses,
310  std::vector<std::unique_ptr<ana::Spectrum>> & spectra,
311  bool fromfile = false) : GenieMultiverseSpectra(nuniverses, spectra, fromfile){};
312 
313  const Spectrum* UpperSigma(BandOptions opt = kBandFromNominal);
314  const Spectrum* LowerSigma(BandOptions opt = kBandFromNominal);
315 
316  /// Save the Spectrum structure of all universe to a TFile directory.
317  void SaveTo(TDirectory* dir, const std::string& name);
318 
319  /// Load the Spectrum structure of all universe from a TFile directory.
320  static std::unique_ptr<GenieMultiverseNormalizedSpectra> LoadFrom(TDirectory* dir, const std::string& name);
321 
322  private:
324 
325  /// A flag signaling whether the spectra are normalized already.
326  bool fSpectraNormalized = false;
327 
328  Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
329  const HistAxis& axis,
330  const Cut& cut,
331  const SystShifts& shift = kNoShift,
332  const Var& wei = kUnweighted);
333 
334  /// Normalize spectra. It has to be done AFTER SpectrumLoader::Go() is called.
335  void NormalizeSpectra();
336  };
337 }
multi_index_container< knob_sampling_mode, indexed_by< ordered_unique< tag< id >, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode, int, id)>, ordered_unique< tag< name >, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode, std::string, name)> > > knob_sampling_mode_set
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH2 * ToTH2(const TH1 *h, const HistAxis &xaxis, const HistAxis &yaxis)
std::vector< std::map< const ISyst *, double > > fShiftTables
Container of shift tables for all universes.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
GenieMultiverseNormalizedSpectra(unsigned int nuniverses, std::vector< std::unique_ptr< ana::Spectrum >> &spectra, bool fromfile=false)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
multi_index_container< knob_sampling_mode, indexed_by< ordered_unique< tag< id >, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode, int, id)>, ordered_unique< tag< name >, BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode, std::string, name)> > > knob_sampling_mode_set
std::vector< std::map< const ISyst *, double > > fShiftTable
knob_sampling_mode(int id_, std::string name_, int mode_)
const auto & GetShiftTable() const
declarations for knob configuration table
float BinSigma(std::vector< double > events, float nsigmas, float pivot)
Definition: doTemplateFit.C:64
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::vector< std::map< const ISyst *, double > > ShiftTables()
Return the tables of knob shift values.
std::vector< std::map< const ISyst *, double > > NuTruthShiftTables()
Int_t col[ntarg]
Definition: Style.C:29
#define pot
knob_sampling_mode_set fKnobConfigTable
loader
Definition: demo0.py:10
const SystShifts kNoShift
Definition: SystShifts.cxx:22
unsigned int GetNUniverses()
Get the number of universes in an object of this class.
Base class for the various types of spectrum loader.
std::vector< SystShifts > fSystShifts
const Cut cut
Definition: exporter_fd.C:30
knob_sampling_mode_set fKnobConfigTable
knob configuration table
std::vector< std::map< const ISyst *, double > > fNuTruthShiftTables
Container of neutrino truth shift tables for all universes.
TDirectory * dir
Definition: macro.C:5
GenieMultiverseSpectra MultiverseSpectra
typedef for backward compatibility
std::vector< SystShifts > GetSystShifts()
const NuTruthVar kNuTruthUnweighted
Definition: Var.h:100
std::vector< const Spectrum * > AllUniverses() const
return all universes
const int nuniverses
void events(int which)
Definition: Cana.C:52
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
unsigned int fNUniverses
number of universes to generate
knob_sampling_mode(int id_, std::string name_, int mode_)
enum BeamMode string