SpectrumPredictionMaker.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 4/21/20.
3 //
4 
8 
15 
16 #include "TH1.h"
17 
18 namespace cmf{
19  static SpectrumPredictionMaker *gLPU = nullptr;
20 
21  //----------------------------------------------------------------------------
23  {
24  if(gLPU == nullptr) gLPU = new SpectrumPredictionMaker();
25 
26  return gLPU;
27  }
28 
29  //----------------------------------------------------------------------------
31  : fParamZ("NoZParam")
32  , fZVal(std::numeric_limits<double>::max())
33  {
34  }
35 
36  //----------------------------------------------------------------------------
38  {
39  }
40 
41  //----------------------------------------------------------------------------
43  {
45 
46  LOG_DEBUG("SpectrumPredictionMaker")
47  << "initial location is\n"
48  << loc;
49 
51  pset.get<fhicl::ParameterSet>("ShifterAndWeighterParameters"));
52 
53  fManipulatorPars.emplace(cmf::kFARDET, pset.get<fhicl::ParameterSet>("FarDetEventListManipulator" ));
54  fManipulatorPars.emplace(cmf::kNEARDET, pset.get<fhicl::ParameterSet>("NearDetEventListManipulator"));
55 
56 
57  // get the grid point parameters
58  auto gridPointPars = pset.get<fhicl::ParameterSet>("GridPointPars");
59 
60  fXVal = gridPointPars.get<float >("GridXVal");
61  fYVal = gridPointPars.get<float >("GridYVal");
62  fParamX = gridPointPars.get<std::string>("ParamX" );
63  fParamY = gridPointPars.get<std::string>("ParamY" );
64 
65  // check to see if we have the z parameter defined
66  if(gridPointPars.has_key("GridZVal")){
67  fZVal = gridPointPars.get<float >("GridZVal");
68  fParamZ = gridPointPars.get<std::string>("ParamZ" );
69  }
70 
71  auto const& hpPS = gridPointPars.get<fhicl::ParameterSet>("HiddenParameters");
72  for(auto const& itr : hpPS.get_names()){
73  auto const& hpConfig = hpPS.get<fhicl::ParameterSet>(itr);
74  fHiddenParameters.emplace_back(hpConfig.get<int >("Divs" ),
75  hpConfig.get<float>("LowerVal"),
76  hpConfig.get<float>("UpperVal"),
77  cmf::StringToOscParmType(hpConfig.get<std::string>("ParName")));
78  }
79 
80  std::stringstream hiddenParsOutput;
81  hiddenParsOutput << "HiddenParameters are:";
82  for(auto const& itr : fHiddenParameters)
83  hiddenParsOutput << "\n\t"
84  << cmf::cOscParams_Strings[itr.oscParm]
85  << " with "
86  << itr.divs
87  << " points between/including "
88  << itr.lowerVal
89  << " to "
90  << itr.upperVal;
91 
92  LOG_VERBATIM("SpectrumPredictionMaker")
93  << hiddenParsOutput.str();
94 
95  // find the points in space and then set the parameters for those locations
96  this->FindOscillationPoints();
97  }
98 
99  // Method to determine the point locations in space.
100  //......................................................................
102  {
103  // find the values for the grid in space defining the hidden parameters
104  std::vector< std::vector<float> > hiddenParVals(fHiddenParameters.size());
105  std::vector< cmf::OscillationParameterMap > oscillationPoints;
106  size_t ctr = 0;
107 
108  // we know we have 1 point for the X and Y parameters in oscillation space
109  size_t numOscPoints = 1;
110  for(auto const& hpcItr : fHiddenParameters){
111 
112  // each of the hidden parameters increases the number of points by a
113  // factor of the number of divisions
114  numOscPoints *= hpcItr.divs;
115 
116  LOG_DEBUG("SpectrumPredictionMaker")
117  << ctr
118  << " hidden param "
119  << cmf::cOscParams_Strings[hpcItr.oscParm]
120  << " "
121  << hpcItr.divs
122  << " "
123  << hpcItr.stepSize
124  << " "
125  << hpcItr.lowerVal
126  << " "
127  << hpcItr.upperVal;
128 
129  for(int d = 0; d < hpcItr.divs; ++d){
130  hiddenParVals[ctr].emplace_back((d * hpcItr.stepSize) + hpcItr.lowerVal);
131 
132  LOG_DEBUG("SpectrumPredictionMaker")
133  << "\t"
134  << d
135  << " "
136  << hiddenParVals[ctr].back();
137  }
138 
139  LOG_DEBUG("SpectrumPredictionMaker")
140  << " there are "
141  << hiddenParVals[ctr].size()
142  << " values stored for this parameter";
143 
144  ++ctr;
145  }
146 
147  LOG_DEBUG("SpectrumPredictionMaker")
148  << "begin creating "
149  << numOscPoints
150  << " points in oscillation space";
151 
152  oscillationPoints.resize(numOscPoints);
153 
154  auto xOscPar = cmf::StringToOscParmType(fParamX);
155  auto yOscPar = cmf::StringToOscParmType(fParamY);
156  for(size_t p = 0; p < numOscPoints; ++p){
157  oscillationPoints[p].emplace(xOscPar, fXVal);
158  oscillationPoints[p].emplace(yOscPar, fYVal);
160  oscillationPoints[p].emplace(cmf::StringToOscParmType(fParamZ), fZVal);
161  }
162 
163  // For now have different loops for the number of hidden points, in the future
164  // try to figure out a smarter solution
165  ctr = 0;
166  if(fHiddenParameters.size() == 1){
167  for(size_t p = 0; p < hiddenParVals[0].size(); ++p)
168  oscillationPoints[p].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][p]);
169  }
170  else if(fHiddenParameters.size() == 2){
171  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
172  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
173  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
174  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
175 
176  LOG_DEBUG("CMFPredictionLibrary")
177  << "hidden parameter location "
178  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
179  << " "
180  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second;
181 
182  ++ctr;
183  }
184  }
185  } // end if 2 hidden parameters
186  else if(fHiddenParameters.size() == 3){
187  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
188  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
189  for(size_t c = 0; c < hiddenParVals[2].size(); ++c){
190  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
191  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
192  oscillationPoints[ctr].emplace(fHiddenParameters[2].oscParm, hiddenParVals[2][c]);
193  LOG_DEBUG("CMFPredictionLibrary")
194  << "hidden point location "
195  << ctr
196  << " "
197  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
198  << " "
199  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second
200  << " "
201  << oscillationPoints[ctr].find(fHiddenParameters[2].oscParm)->second;
202 
203  ++ctr;
204  }
205  }
206  }
207  } // end if 3 hidden parameters
208  else if(fHiddenParameters.size() == 4){
209  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
210  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
211  for(size_t c = 0; c < hiddenParVals[2].size(); ++c){
212  for(size_t d = 0; d < hiddenParVals[3].size(); ++d){
213  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
214  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
215  oscillationPoints[ctr].emplace(fHiddenParameters[2].oscParm, hiddenParVals[2][c]);
216  oscillationPoints[ctr].emplace(fHiddenParameters[3].oscParm, hiddenParVals[3][d]);
217 
218  LOG_DEBUG("CMFPredictionLibrary")
219  << "hidden point location "
220  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
221  << " "
222  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second
223  << " "
224  << oscillationPoints[ctr].find(fHiddenParameters[2].oscParm)->second
225  << " "
226  << oscillationPoints[ctr].find(fHiddenParameters[3].oscParm)->second;
227 
228  ++ctr;
229  }
230  }
231  }
232  }
233  } // end if 4 hidden parameters
234 
235  LOG_VERBATIM("SpectrumPredictionMaker")
236  << "Hidden points are: ";
237 
238  std::string point;
239  for(auto const& itr : oscillationPoints){
240  fOscillationPoints.emplace_back(itr);
241  point.clear();
242  for(auto const& pos : itr){
243  point += cmf::cOscParams_Strings[pos.first] + " " + std::to_string(pos.second) + " ";
244  }
245  LOG_VERBATIM("SpectrumPredictionMaker")
246  << point;
247  }
248 
249  } // end FindHiddenPointLocations
250 
251 
252  //......................................................................
253  std::vector<float> SpectrumPredictionMaker::HiddenParStep() const
254  {
255  std::vector<float> steps;
256  for(auto const& itr : fHiddenParameters) steps.emplace_back(itr.stepSize);
257  return std::move(steps);
258  }
259 
260  //......................................................................
261  void SpectrumPredictionMaker::MakePredictions(std::vector<cmf::OscParamPoint> & predVec)
262  {
263  // containers we will need
264  cmf::EventListColl mcLists;
265  std::vector<cmf::Spectrum> pointSpectrum(fOscillationPoints.size());
266  std::vector<float> count;
267 
268  std::set<cmf::DetType_t> dets;
269  for(auto const& mpItr : fManipulatorPars){
270 
271  mcLists.clear();
272 
273  dets.clear();
274  dets.insert(mpItr.first);
275 
276  LOG_VERBATIM("PredictionLibrary")
277  << cmf::cDetType_Strings[mpItr.first]
278  << "\n"
279  << mpItr.second.to_indented_string();
280 
281  // get the manipulator to deserialize the data and MC
282  // we'll keep the MC lists around, but will drop the data lists as
283  // soon as we leave this method
284  cmf::EventListManipulator elm(mpItr.second);
285  elm.Deserialize(mcLists, cmf::kMC, dets);
286 
287  auto const& exposureMap = elm.ExposureMap();
288 
289  // print out the exposures we are requesting
290 // for(auto const& itr : exposureMap)
291 // LOG_VERBATIM("PredictionLibrary")
292 // << "Exposure for "
293 // << itr.first.ToString()
294 // << " is "
295 // << itr.second;
296 
297  for(size_t p = 0; p < fOscillationPoints.size(); ++p){
298 
300 
301  cmf::FillSpectrum(mcLists,
302  exposureMap,
303  pointSpectrum[p],
304  count);
305 
306  // if(mpItr.first == cmf::kFARDET){
307  // for(size_t b = 0; b < pointSpectrum[p].size(); ++b)
308  // LOG_VERBATIM("PredictionLibrary")
309  // << "bin "
310  // << b
311  // << " has "
312  // << pointSpectrum[p][b];
313  // }
314 
315  LOG_VERBATIM("SpectrumPredictionMaker")
316  << " Finished point "
317  << cmf::cDetType_Strings[mpItr.first]
318  << " "
319  << p;
320  } // end loop to make predictions
321  } // end loop to load one detector at a time
322 
323  predVec.clear();
324  for(size_t p = 0; p < fOscillationPoints.size(); ++p){
325  predVec.emplace_back(fOscillationPoints[p].AsOscillationParameterMap(),
326  pointSpectrum[p]);
327  this->MakeHiddenParameterHistograms(fOscillationPoints[p].AsOscillationParameterMap(),
328  pointSpectrum[p]);
329  }
330  }
331 
332  //......................................................................
334  cmf::Spectrum const& spectrum)
335  {
336  // store the spectrum in a histogram that is written out to a file
337  // loop over the fHiddenPointLocs vectors and create a
338  // histogram for each location
340 
341  // first off, let's make a new directory for this universe
342  art::TFileDirectory dir = tfs->mkdir("Spectra");
343 
345 
346  histName = "prediction";
347  for(auto const& hpitr : fHiddenParameters){
348  auto const& itr = parMap.find(hpitr.oscParm);
349  histName += "_" + cmf::cOscParams_Strings[itr->first] + "_" +
350  std::to_string(itr->second);
351  }
352 
353  LOG_DEBUG("CMFPredictionLibrary")
354  << histName;
355 
356  TH1D* hist = dir.make<TH1D>(histName.c_str(),
357  ";Logical Bin;Events",
359  0,
361 
362  for(size_t b = 0; b < spectrum.size(); ++b){
363  hist->Fill(b, spectrum[b]);
364 
365  LOG_DEBUG("SpectrumPredictionMaker")
366  << "filling "
367  << hist->GetName()
368  << b
369  << " "
370  << spectrum[b]
371  << " "
372  << hist->GetBinContent(hist->FindBin(b));
373  }
374 
375  }
376 
377 
378 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void MakePredictions(std::vector< cmf::OscParamPoint > &predVec)
const char * p
Definition: xmltok.h:285
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
static SpectrumPredictionMaker * gLPU
const std::string cOscParams_Strings[kNumOscParams]
Definition: Constants.h:253
static ParameterUtility * Instance()
std::map< cmf::DetType_t, fhicl::ParameterSet > fManipulatorPars
parameters to make object for deserializing events
void Initialize(fhicl::ParameterSet const &pset)
std::map< cmf::OscParm_t, float > OscillationParameterMap
Definition: Constants.h:572
void Deserialize(cmf::EventListColl &eventLists, cmf::DataMC_t dataMC=cmf::kBoth, std::set< cmf::DetType_t > const &detectors=std::set< cmf::DetType_t >({cmf::kNEARDET, cmf::kFARDET}))
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
Float_t d
Definition: plot.C:236
cmf::ExposureMap const & ExposureMap() const
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, std::vector< float > &spectrum, std::vector< float > &count, cmf::InteractionType_t intType, bool applyExposureNorm)
static ShifterAndWeighter * Instance()
std::vector< float > Spectrum
Definition: Constants.h:570
static SpectrumPredictionMaker * Instance()
size_t TotalBins(bool allSels=false)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
void MakeHiddenParameterHistograms(cmf::OscillationParameterMap const &parMap, cmf::Spectrum const &spectrum)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
T * make(ARGS...args) const
std::vector< cmf::HiddenParameterConfig > fHiddenParameters
marginalized oscillation parameters
TDirectory * dir
Definition: macro.C:5
const hit & b
Definition: hits.cxx:21
cmf::Location ParametersAsLocation()
std::vector< float > HiddenParStep() const
#define LOG_VERBATIM(category)
std::map< cmf::MetaData, cmf::EventList > EventListColl
Definition: Event.h:147
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const std::string cDetType_Strings[5]
Definition: Constants.h:552
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::vector< cmf::Location > fOscillationPoints
points in oscillation space where we want predictions
void SetCurrentVals(cmf::Location const &loc)
static cmf::OscParm_t StringToOscParmType(std::string const &str)
static CovarianceBinUtility * Instance()