SpectrumPredictionMaker.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 4/21/20.
3 //
4 
6 #include "art_root_io/TFileService.h"
8 
16 
17 #include "TH1.h"
18 #include "TStopwatch.h"
19 
20 namespace cmf{
21  static SpectrumPredictionMaker *gLPU = nullptr;
22 
23  //----------------------------------------------------------------------------
25  {
26  if(gLPU == nullptr) gLPU = new SpectrumPredictionMaker();
27 
28  return gLPU;
29  }
30 
31  //----------------------------------------------------------------------------
33  {
34  }
35 
36  //----------------------------------------------------------------------------
38  {
39  }
40 
41  //----------------------------------------------------------------------------
43  {
45 
46  MF_LOG_DEBUG("SpectrumPredictionMaker")
47  << "initial location is\n"
48  << loc;
49 
51  pset.get<fhicl::ParameterSet>("ShifterAndWeighterParameters"));
52 
53  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
54  fManipulatorPars.emplace(detItr, pset.get<fhicl::ParameterSet>(cmf::cDetType_Strings[detItr] + "EventListManipulator" ));
55  }
56 
57 
58  // get the grid point parameters
59  auto gridPointPars = pset.get<fhicl::ParameterSet>("GridPointPars");
60 
61  auto gridPars = gridPointPars.get< std::vector<std::string> >("GridPars");
62  for(size_t p = 0; p < gridPars.size(); ++p)
63  fOscParams.emplace_back(cmf::StringToOscParmType(gridPars[p]));
64 
65  fOscParamVals = gridPointPars.get< std::vector<double> >("GridParVals");
66 
67  if(fOscParamVals.size() != fOscParams.size())
68  throw cet::exception("SpectrumPredictionMaker") << "grid parameter vector size"
69  << " does't match grid parameter val vector size";
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  MF_LOG_VERBATIM("SpectrumPredictionMaker")
93  << hiddenParsOutput.str();
94 
95  // find the points in space and then set the parameters for those locations
96  this->FindOscillationPoints(loc);
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 
115  if(hpcItr.oscParm != cmf::kd24){
116  for(int d = 0; d < hpcItr.divs; ++d){
117  hiddenParVals[ctr].emplace_back((d * hpcItr.stepSize) + hpcItr.lowerVal);
118 
119  MF_LOG_DEBUG("SpectrumPredictionMaker")
120  << "\t"
121  << d
122  << " "
123  << hiddenParVals[ctr].back();
124  }
125  }
126  else{
127  // d_24 comes in to the probability as sin(d_24)
128  // that means 0, pi and 2pi are degenerate with each other
129  // so are values between 0 and pi/2 and values between pi/2 and pi,
130  // as are values between pi and 3pi/2 with values between 3pi/2 and 2pi
131  hiddenParVals[ctr] = std::vector<float>({0., 0.785398, 1.570796, 3.926991, 4.712389});
132  }
133 
134  numOscPoints *= hiddenParVals[ctr].size();
135 
136  MF_LOG_DEBUG("SpectrumPredictionMaker")
137  << ctr
138  << " hidden param "
139  << cmf::cOscParams_Strings[hpcItr.oscParm]
140  << " "
141  << hpcItr.divs
142  << " "
143  << hpcItr.stepSize
144  << " "
145  << hpcItr.lowerVal
146  << " "
147  << hpcItr.upperVal
148  << " has "
149  << hiddenParVals[ctr].size()
150  << " values stored for this parameter";
151 
152  ++ctr;
153  }
154 
155  MF_LOG_DEBUG("SpectrumPredictionMaker")
156  << "begin creating "
157  << numOscPoints
158  << " points in oscillation space";
159 
160  oscillationPoints.resize(numOscPoints);
161 
162  for(size_t p = 0; p < numOscPoints; ++p){
163  for(size_t op = 0; op < fOscParams.size(); ++op)
164  oscillationPoints[p].emplace(fOscParams[op], fOscParamVals[op]);
165  }
166 
167  // For now have different loops for the number of hidden points, in the future
168  // try to figure out a smarter solution
169  ctr = 0;
170  if(fHiddenParameters.size() == 1){
171  for(size_t p = 0; p < hiddenParVals[0].size(); ++p)
172  oscillationPoints[p].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][p]);
173  }
174  else if(fHiddenParameters.size() == 2){
175  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
176  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
177  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
178  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
179 
180  MF_LOG_DEBUG("CMFPredictionLibrary")
181  << "hidden parameter location "
182  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
183  << " "
184  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second;
185 
186  ++ctr;
187  }
188  }
189  } // end if 2 hidden parameters
190  else if(fHiddenParameters.size() == 3){
191  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
192  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
193  for(size_t c = 0; c < hiddenParVals[2].size(); ++c){
194  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
195  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
196  oscillationPoints[ctr].emplace(fHiddenParameters[2].oscParm, hiddenParVals[2][c]);
197  MF_LOG_DEBUG("CMFPredictionLibrary")
198  << "hidden point location "
199  << ctr
200  << " "
201  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
202  << " "
203  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second
204  << " "
205  << oscillationPoints[ctr].find(fHiddenParameters[2].oscParm)->second;
206 
207  ++ctr;
208  }
209  }
210  }
211  } // end if 3 hidden parameters
212  else if(fHiddenParameters.size() == 4){
213  for(size_t a = 0; a < hiddenParVals[0].size(); ++a){
214  for(size_t b = 0; b < hiddenParVals[1].size(); ++b){
215  for(size_t c = 0; c < hiddenParVals[2].size(); ++c){
216  for(size_t d = 0; d < hiddenParVals[3].size(); ++d){
217  oscillationPoints[ctr].emplace(fHiddenParameters[0].oscParm, hiddenParVals[0][a]);
218  oscillationPoints[ctr].emplace(fHiddenParameters[1].oscParm, hiddenParVals[1][b]);
219  oscillationPoints[ctr].emplace(fHiddenParameters[2].oscParm, hiddenParVals[2][c]);
220  oscillationPoints[ctr].emplace(fHiddenParameters[3].oscParm, hiddenParVals[3][d]);
221 
222  MF_LOG_DEBUG("CMFPredictionLibrary")
223  << "hidden point location "
224  << oscillationPoints[ctr].find(fHiddenParameters[0].oscParm)->second
225  << " "
226  << oscillationPoints[ctr].find(fHiddenParameters[1].oscParm)->second
227  << " "
228  << oscillationPoints[ctr].find(fHiddenParameters[2].oscParm)->second
229  << " "
230  << oscillationPoints[ctr].find(fHiddenParameters[3].oscParm)->second;
231 
232  ++ctr;
233  }
234  }
235  }
236  }
237  } // end if 4 hidden parameters
238 
239  MF_LOG_VERBATIM("SpectrumPredictionMaker")
240  << "Oscillation points are: ";
241 
242  std::string point;
243  for(auto const& itr : oscillationPoints){
244  fOscillationPoints.emplace_back(defaultLoc);
245 
246  // use cmf::kALLDET when setting these oscillation parameters because
247  // only the baseline changes with detector
248  point.clear();
249  for(auto const& parItr : itr){
250  fOscillationPoints.back().SetParameterValue(cmf::cOscParams_Strings[parItr.first],
251  parItr.second,
252  cmf::kALLDET);
253 
254  point += cmf::cOscParams_Strings[parItr.first] + " " + std::to_string(parItr.second) + " ";
255  }
256 
257  MF_LOG_VERBATIM("SpectrumPredictionMaker")
258  << point;
259  }
260 
261  } // end FindHiddenPointLocations
262 
263 
264  //......................................................................
265  std::vector<float> SpectrumPredictionMaker::HiddenParStep() const
266  {
267  std::vector<float> steps;
268  for(auto const& itr : fHiddenParameters) steps.emplace_back(itr.stepSize);
269  return std::move(steps);
270  }
271 
272  //......................................................................
273  void SpectrumPredictionMaker::MakePredictions(std::vector<cmf::OscParamPoint> & predVec)
274  {
275  // containers we will need
276  cmf::EventListColl mcLists;
277  std::vector<cmf::Spectrum> pointSpectrum(fOscillationPoints.size());
279 
280  std::set<cmf::DetType_t> dets;
281  for(auto const& mpItr : fManipulatorPars){
282 
283  mcLists.clear();
284 
285  dets.clear();
286  dets.insert(mpItr.first);
287 
288  MF_LOG_VERBATIM("PredictionLibrary")
289  << cmf::cDetType_Strings[mpItr.first]
290  << "\n"
291  << mpItr.second.to_indented_string();
292 
293  //TStopwatch sw;
294 
295  //sw.Start();
296  // get the manipulator to deserialize the data and MC
297  // we'll keep the MC lists around, but will drop the data lists as
298  // soon as we leave this method
299  cmf::EventListManipulator elm(mpItr.second);
300  elm.Deserialize(mcLists, cmf::kMC, dets);
301 
302  auto const& exposureMap = elm.ExposureMap();
303 
304  //sw.Stop();
305 
306  // MF_LOG_VERBATIM("PredictionLibrary")
307  // << "deserializing takes ";
308  // sw.Print();
309  // sw.Reset();
310 
311  // print out the exposures we are requesting
312 // for(auto const& itr : exposureMap)
313 // MF_LOG_VERBATIM("PredictionLibrary")
314 // << "Exposure for "
315 // << itr.first.ToString()
316 // << " is "
317 // << itr.second;
318 
319  for(size_t p = 0; p < fOscillationPoints.size(); ++p){
320 
321  // sw.Start();
323  // sw.Stop();
324 
325  // MF_LOG_VERBATIM("PredictionLibrary")
326  // << "setting current values takes ";
327  // sw.Print();
328  // sw.Reset();
329 
330  // sw.Start();
331  cmf::FillSpectrum(mcLists,
332  exposureMap,
333  pointSpectrum[p],
334  count);
335 
336  // sw.Stop();
337  // MF_LOG_VERBATIM("PredictionLibrary")
338  // << "filling spectrum takes ";
339  // sw.Print();
340  // sw.Reset();
341 
342  // if(mpItr.first == cmf::kFARDET){
343  // for(size_t b = 0; b < pointSpectrum[p].size(); ++b)
344  // MF_LOG_VERBATIM("PredictionLibrary")
345  // << "bin "
346  // << b
347  // << " has "
348  // << pointSpectrum[p][b];
349  // }
350 
351  MF_LOG_VERBATIM("SpectrumPredictionMaker")
352  << " Finished point "
353  << cmf::cDetType_Strings[mpItr.first]
354  << " "
355  << p;
356  } // end loop to make predictions
357  } // end loop to load one detector at a time
358 
359  MF_LOG_DEBUG("SpectrumPredictionMaker")
360  << fOscillationPoints.front();
361 
362  predVec.clear();
363  for(size_t p = 0; p < fOscillationPoints.size(); ++p){
364  predVec.emplace_back(fOscillationPoints[p].AsOscillationParameterMap(),
365  pointSpectrum[p]);
366 
367  MF_LOG_DEBUG("SpectrumPredictionMaker")
368  << predVec.back();
369 
370  this->MakeHiddenParameterHistograms(fOscillationPoints[p].AsOscillationParameterMap(),
371  pointSpectrum[p]);
372  }
373  }
374 
375  //......................................................................
377  cmf::Spectrum const& spectrum)
378  {
379  // store the spectrum in a histogram that is written out to a file
380  // loop over the fHiddenPointLocs vectors and create a
381  // histogram for each location
383 
384  // first off, let's make a new directory for this universe
385  std::string dirName("Predictions_");
386 
387  for(size_t op = 0; op < fOscParams.size(); ++op)
388  dirName += cmf::cOscParams_Strings[fOscParams[op]] + "_" + std::to_string(fOscParamVals[op]);
389 
390  art::TFileDirectory dir = tfs->mkdir(dirName);
391 
393 
394  histName = "prediction";
395  for(auto const& hpitr : fHiddenParameters){
396  auto const& itr = parMap.find(hpitr.oscParm);
397  histName += "_" + cmf::cOscParams_Strings[itr->first] + "_" +
398  std::to_string(itr->second);
399  }
400 
401  MF_LOG_DEBUG("CMFPredictionLibrary")
402  << histName;
403 
405  histName,
406  ";Logical Bin;Events",
407  spectrum);
408 
409  // previous way of filling the histogram, leave it here for the
410  // debugging output
411  // for(size_t b = 0; b < spectrum.size(); ++b){
412  // hist->Fill(b, spectrum[b]);
413 
414  // MF_LOG_DEBUG("SpectrumPredictionMaker")
415  // << "filling "
416  // << hist->GetName()
417  // << b
418  // << " "
419  // << spectrum[b]
420  // << " "
421  // << hist->GetBinContent(hist->FindBin(b));
422  // }
423 
424  }
425 
426 
427 }
void MakePredictions(std::vector< cmf::OscParamPoint > &predVec)
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, cmf::Spectrum &spectrum, cmf::Spectrum &count, cmf::InteractionType_t intType, bool applyExposureNorm)
std::vector< cmf::OscParm_t > fOscParams
oscillation parameters
std::vector< double > fOscParamVals
oscillation parameter values
const char * p
Definition: xmltok.h:285
const std::vector< std::string > cDetType_Strings({"UnknownDet","NearDet","FarDet","MINOSNear","MINOSFar","AllDetectors"})
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:746
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
static SpectrumPredictionMaker * gLPU
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}))
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:748
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
static ShifterAndWeighter * Instance()
static SpectrumPredictionMaker * Instance()
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::vector< cmf::EventList > EventListColl
Definition: Event.h:150
void MakeHiddenParameterHistograms(cmf::OscillationParameterMap const &parMap, cmf::Spectrum const &spectrum)
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
std::string dirName
Definition: PlotSpectra.h:47
std::vector< cmf::HiddenParameterConfig > fHiddenParameters
marginalized oscillation parameters
TDirectory * dir
Definition: macro.C:5
#define MF_LOG_VERBATIM(category)
#define MF_LOG_DEBUG(id)
const hit & b
Definition: hits.cxx:21
static PlotUtilities * Instance()
cmf::Location ParametersAsLocation()
std::vector< float > HiddenParStep() const
std::vector< cmf::Location > fOscillationPoints
points in oscillation space where we want predictions
void SetCurrentVals(cmf::Location const &loc)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
void FindOscillationPoints(cmf::Location const &defaultLoc)
static cmf::OscParm_t StringToOscParmType(std::string const &str)
const std::vector< std::string > cOscParams_Strings({"L","Rho","Dmsq21","Dmsq32","Th12","Th13","Th23","dCP","Th14","Th24","Th34","Dmsq41","d24","Eps_ee","Eps_emu","Eps_etau","Eps_mumu","Eps_mutau","Eps_tautau","Delta_emu","Delta_etau","Delta_mutau"})
enum BeamMode string