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