RandomUniverseUtility.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 4/20/20.
3 //
4 
5 // Framework includes
7 #include "art_root_io/TFileService.h"
8 #include "art_root_io/TFileDirectory.h"
10 #include "cetlib_except/exception.h"
11 
15 
16 #include "TH1.h"
17 
18 namespace cmf{
19 
20  static RandomUniverseUtility *gRUU = nullptr;
21 
22  //----------------------------------------------------------------------------
24  {
25  if(gRUU == nullptr) gRUU = new RandomUniverseUtility();
26 
27  return gRUU;
28  }
29 
30  //----------------------------------------------------------------------------
32  {
33  }
34 
35  //----------------------------------------------------------------------------
37  {
38  }
39 
40  //----------------------------------------------------------------------------
42  {
43  // set the seed for the random number to be what is configured, plus the
44  // start universe value in case lots are jobs are submitted with the same seed
45  // by accident
46  fRand.SetSeed(pset.get<unsigned long int>("RandomSeed", 0) +
47  pset.get<unsigned int >("StartUniverse", 0));
48 
49  fNumUniverses = pset.get<unsigned int>("NumUniverses", 1);
50  fPoissonFluctuate = pset.get<bool >("PoissonFluctuate", false);
51  fPoissonByHist = pset.get<bool >("PoissonDrawWithHistogram", false);
52  fPoissonByBin = pset.get<bool >("PoissonDrawWithBin", true );
53 
54  // fill histogram with systematic names and which values were used
56 
57  fSystDrawHist = tfs->make<TH2D>("SystematicParameterDraws",
58  ";Universe;Uncertainty Draw;",
60  0,
62  2,
63  0,
64  2);
65  fSystDrawHist->SetCanExtend(TH1::kAllAxes);
66  fSystDrawHist->SetStats(false);
67 
68  fOscParamDrawHist = tfs->make<TH2D>("OscillationParameterDraws",
69  ";Unvierse;Parameter Draw",
71  0,
73  2,
74  0,
75  2);
76  fOscParamDrawHist->SetCanExtend(TH1::kAllAxes);
77  fOscParamDrawHist->SetStats(0);
78  }
79 
80  //......................................................................
82  std::vector<cmf::Location> & vecLoc)
83  {
84  // check that the vector is the correct size
85  vecLoc.resize(fNumUniverses, nominalLoc);
86 
87  for(unsigned int u = 0; u < fNumUniverses; ++u){
88 
89  this->CreateLocationWithVariedSystematics(nominalLoc, vecLoc[u], u);
90  } // end loop over the number of universes
91 
92  // fSystDrawHist->LabelsDeflate("Y");
93  // fSystDrawHist->LabelsOption("h", "Y");
94 
95  }
96 
97  // This function will pick a random value for the specified parameter
98  // from a guassian distribution. The random value will be within +/- 1 sigma
99  // of the nominal value
100  //......................................................................
102  cmf::Location & vecLoc,
103  size_t universeNum)
104  {
105  // we will pick a random value from a unit gaussian, ie gaussian with
106  // sigma = 1. That will be the fractional change on the sigma we want
107  // to apply to the new input point
108  double sigFrac = 0.;
109 
110  // get the names of the systematic parameters in the ND and FD
111  std::map<std::string, double> systParVals;
112  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
113  for(auto const& itr : nominalLoc.DetectorLocation(detItr).Parameters()){
114  if(itr.IsOscPar() ||
115  itr.IsFixed() ||
116  systParVals.count(itr.Name()) ) continue;
117 
118  // each systematic parameter should have a different random number
119  // change in sigma so that we don't artificially make them change
120  // in a correlated way if doing more than one parameter in a job
121  sigFrac = (this->UniformDrawParameters(itr.Name())) ? fRand.Rndm() : fRand.Gaus();
122 
123  systParVals.emplace(itr.Name(), sigFrac);
124  }
125  }
126 
127  // prepare some output for the random values in this universe
128  std::string systStr("Systematic Parameters " + std::to_string(universeNum));
129 
130  for(auto const& itr : systParVals){
131 
132  systStr += "\n\t" + itr.first + ": " + std::to_string(itr.second);
133 
134  fSystDrawHist->Fill(universeNum, itr.first.c_str(), itr.second);
135 
136  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
137  vecLoc.SetParameterValue(itr.first, itr.second, detItr);
138  }
139  } // end loop over systematic parameters to change
140 
141  MF_LOG_VERBATIM("RandomUniverseUtility")
142  << systStr;
143  }
144 
145  //......................................................................
147  cmf::Location & vecLoc,
148  size_t universeNum)
149  {
150  // loop over all the parameters and vary the ones we want
151  double val;
153 
154  // set the location to the nominal location
155  vecLoc = nominalLoc;
156 
157  // prepare some output for the random values in this universe
158  std::string oscStr("Oscillation Parameters " + std::to_string(universeNum) + ": ");
159 
160  // get the names of the oscillation parameters. Only the baseline parameter
161  // varies from one detector to the next. By default we get the NOvA FD
162  std::map<std::string, double> oscParVals;
163  std::map<std::string, bool> oscParFixed;
164  for(auto const& itr : nominalLoc.OscillationParameters()){
165 
166  name = itr.Name();
167  val = itr.Value();
168 
169  MF_LOG_DEBUG("RandomUniverseUtility")
170  << name
171  << " is osc par: "
172  << itr.IsOscPar()
173  << " is fixed: "
174  << itr.IsFixed();
175 
176  // The distance the neutrino travels should never be adjusted, so
177  // don't bother putting it into oscParVals
178  if(itr.IsOscPar() && (cmf::OscParm_t)itr.Key() != cmf::kL){
179  if(!itr.IsFixed()){
180  // we have a parameter that is not fixed, so let's figure out how to
181  // vary it
182  // first check if it is in the list of nuisance parameters, if not
183  // look to see if it is constrained.
184  if(itr.IsNuisance()){
185  // treat the parameter as Gaussian distributed
186  val = fRand.Gaus(itr.CentralValue(),
187  itr.Variance());
188 
189  MF_LOG_DEBUG("RandomUniverseUtility")
190  << "nuisance parameter "
191  << name
192  << " val "
193  << val
194  << " "
195  << itr.CentralValue()
196  << " "
197  << itr.Variance()
198  << " "
199  << itr.CentralValue() - itr.Variance();
200  }
201  else if(itr.IsConstrained()){
202  val = fRand.Rndm() * (itr.UpperBound() - itr.LowerBound()) +
203  itr.LowerBound();
204 
205  MF_LOG_DEBUG("RandomUniverseUtility")
206  << "non-nuisance parameter "
207  << name
208  << " val "
209  << val
210  << " "
211  << itr.LowerBound()
212  << " "
213  << itr.UpperBound()
214  << " "
215  << itr.UpperBound() - itr.LowerBound();
216  }
217 
218  oscStr += name + ": " + std::to_string(val) + "\t";
219  }
220  oscParVals.emplace(name, val);
221  oscParFixed.emplace(name, itr.IsFixed());
222  } // end if we have an oscillation parameter
223  } // end loop over parameters
224 
225  //loop over oscParVals and set the oscillation parameter values that have
226  // changed
227  for(auto const& itr : oscParVals){
228 
229  // only add the parameter to the draw histogram if it isn't fixed
230  if(!oscParFixed.find(itr.first)->second){
231  fOscParamDrawHist->Fill(universeNum, itr.first.c_str(), itr.second);
232  }
233 
234  // do this for both fixed and not fixed parameters to ensure we have
235  // all oscillation parameters set
236  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
237  vecLoc.SetParameterValue(itr.first, itr.second, detItr);
238  }
239  } // end loop over oscillation parameters to change
240 
241  MF_LOG_DEBUG("RandomUniverseUtility")
242  << vecLoc;
243  }
244 
245  //......................................................................
247  std::vector<cmf::Location> & vecLocs)
248  {
249  // check that the vector is the correct size
250  vecLocs.resize(fNumUniverses, nominalLoc);
251 
252  for(unsigned int u = 0; u < fNumUniverses; ++u){
253 
254  this->CreateLocationWithVariedOscParameters(nominalLoc, vecLocs[u], u);
255 
256  MF_LOG_DEBUG("RandomUniverseUtility")
257  << " new point "
258  << u;
259  } // end loop over the number of universes
260 
261  // fOscParamDrawHist->LabelsDeflate("Y");
262  // fOscParamDrawHist->LabelsOption("h", "Y");
263 
264  }
265 
266 
267  // Function to fill the spectra for each meta data type
268  //......................................................................
270  cmf::Spectrum & poissonSpectrum)
271  {
272  // just copy the spectrum into the poissonSpectrum if we aren't configured for this
273  if(!fPoissonFluctuate){
274  poissonSpectrum = spectrum;
275  return;
276  }
277 
278  // first reset the poisson spectrum to be all 0's
279  poissonSpectrum.resize(spectrum.size(), 0.);
280 
281  if(fPoissonByBin) this->PoissonSpectrumBinByBin(spectrum, poissonSpectrum);
282 
283  // for(size_t b = 0; b < poissonSpectrum.size(); ++b){
284  // MF_LOG_VERBATIM("RandomUniverseUtility")
285  // << "FillPoisonSpectrum bin "
286  // << b
287  // << " detector "
288  // << cmf::cDetType_Strings[cmf::CovarianceBinUtility::Instance()->BinToDetector(b)]
289  // << " without fluctuations "
290  // << spectrum[b]
291  // << " poisson expectation "
292  // << poissonSpectrum[b];
293  // }
294 
295  }
296 
297  // poisson fluctuates each bin in the spectra
298  //......................................................................
300  cmf::Spectrum & poissonSpectrum)
301  {
302  for(size_t b = 0; b < spectrum.size(); ++b){
303  poissonSpectrum[b] = fRand.Poisson(spectrum[b]);
304 
305  MF_LOG_DEBUG("RandomUniverseUtility")
306  << "PoissonSpectrumBinByBin, bin: "
307  << b
308  << " "
310  << " "
311  << spectrum[b]
312  << "->"
313  << poissonSpectrum[b]
314  << " ratio: "
315  << 1. - poissonSpectrum[b] / spectrum[b];
316 
317  }
318  }
319 
320  //......................................................................
322  unsigned int uNumber)
323  {
325 
326  auto unumStr = std::to_string(uNumber);
327  // first off, let's make a new directory for this universe
328  art::TFileDirectory dir = tfs->mkdir("Universe_" + unumStr);
329 
330  std::string title = "Universe " + unumStr;
331  std::string name = "PoissonUniverse_" + unumStr;
332 
333  title += ";Bin;Events";
334  auto poissonHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir,
335  name,
336  title,
337  fakeU.PoissonSpectrum());
338 
339  // comment out the next lines as FakeUniverse no longer carries the raw count
340  // information. Leave these here in case we want to pass in the vector from
341  // elsewhere
342 // name = "RawCount_" + unumStr;
343 // TH1D *countHist = dir.make<TH1D>(name.data(),
344 // title.data(),
345 // fakeU.EventCounts().size() - 1,
346 // 0.,
347 // fakeU.EventCounts().size() - 1);
348 
349  name = "AsimovUniverse_" + unumStr;
350  auto asimovHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir,
351  name,
352  title,
353  fakeU.AsimovSpectrum());
354 
355  name = "RatioToAsimovUniverse_" + unumStr;
356  auto ratioHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir,
357  name,
358  title,
359  fakeU.PoissonSpectrum());
360 
361  // for(size_t b = 0; b < fakeU.PoissonSpectrum().size(); ++b){
362  // poissonHist->SetBinContent(b + 1, fakeU.PoissonSpectrum()[b]);
363  // asimovHist ->SetBinContent(b + 1, fakeU.AsimovSpectrum()[b] );
364  // //countHist ->SetBinContent(b + 1, fakeU.EventCounts()[b] );
365  // }
366 
367  ratioHist->Divide(asimovHist);
368  }
369 
370 }
const XML_Char * name
Definition: expat.h:151
enum cmf::osc_params OscParm_t
std::vector< cmf::Parameter > const & Parameters() const
Definition: Parameter.h:43
void Initialize(fhicl::ParameterSet const &pset)
cmf::Spectrum const & PoissonSpectrum() const
TH2D * fSystDrawHist
keep track of the systematic parameter values
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:743
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
void CreateLocationWithVariedOscParameters(cmf::Location const &nominalPoint, cmf::Location &vecPoint, size_t universeNum)
void CreateLocationWithVariedSystematics(cmf::Location const &nominalLoc, cmf::Location &vecLoc, size_t universeNum)
cmf::Spectrum const & AsimovSpectrum() const
T get(std::string const &key) const
Definition: ParameterSet.h:231
unsigned int fNumUniverses
number of fake universes to create
void MakeUniverseHistograms(cmf::FakeUniverse &fakeU, unsigned int uNumber)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
std::set< cmf::DetType_t > const & DetectorsToUse()
bool fPoissonByBin
use bin by bin method for picking poisson numbers
void CreateLocationsWithVariedSystematics(cmf::Location const &nominalLoc, std::vector< cmf::Location > &vecLocs)
static RandomUniverseUtility * Instance()
TH2D * fOscParamDrawHist
keep track of the oscillation parameter values
TDirectory * dir
Definition: macro.C:5
bool UniformDrawParameters(std::string const &name)
#define MF_LOG_VERBATIM(category)
bool fPoissonByHist
use histogram method for picking poisson numbers
void PoissonSpectrumBinByBin(cmf::Spectrum &spectrum, cmf::Spectrum &poissonSpectrum)
#define MF_LOG_DEBUG(id)
const hit & b
Definition: hits.cxx:21
static PlotUtilities * Instance()
cmf::ParameterSpaceLoc const & DetectorLocation(cmf::DetType_t const &det) const
Definition: Parameter.cxx:339
std::vector< Parameter > OscillationParameters(cmf::DetType_t const &det=cmf::kFARDET) const
Definition: Parameter.cxx:357
TRandom3 fRand
RNG for picking the variation in parameters.
bool fPoissonFluctuate
do we want to poisson fluctuate the FD prediction?
void FillPoissonSpectrum(cmf::Spectrum &spectrum, cmf::Spectrum &poissonSpectrum)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
void SetParameterValue(std::string const &parName, double const &value, cmf::DetType_t const &det)
Definition: Parameter.cxx:302
static std::string KeyToString(long const &key)
Definition: StaticFuncs.h:227
static CovarianceBinUtility * Instance()
static RandomUniverseUtility * gRUU
void CreateLocationsWithVariedOscParameters(cmf::Location const &nominalLoc, std::vector< cmf::Location > &vecLoc)
enum BeamMode string