CMFContourFromLibrary_plugin.cc
Go to the documentation of this file.
1 //
2 // Plugin to create contours from library predictions.
3 // Pass the art job a list of input library prediction files for the space
4 // along with a best fit file and a file of random universes
5 //
6 // Created by Brian Rebel on 2/24/20.
7 //
8 
9 #include <string>
10 #include <set>
11 #include <map>
12 #include <memory>
13 #include <sstream>
14 
15 //#include "TStopwatch.h"
16 #include "TMath.h"
17 #include "TGraph2D.h"
18 
19 // Framework includes
25 #include "art_root_io/TFileService.h"
26 #include "art_root_io/TFileDirectory.h"
27 #include "fhiclcpp/ParameterSet.h"
29 #include "cetlib_except/exception.h"
30 
31 // NOvA includes
40 
41 namespace cmf{
42 
44  public:
45  explicit ContourFromLibrary(fhicl::ParameterSet const& pset);
46 
47  ~ContourFromLibrary() override = default;
48 
49  void readResults(art::Results const& r) override;
50 
51  void writeResults(art::Results& r) override;
52 
53  void reconfigure(fhicl::ParameterSet const& p);
54 
55  void clear() override;
56  void endJob() override;
57 
58  private:
59 
60  void ContourForUniverse(std::string const& baseName,
61  cmf::Spectrum const& dataSpectrum);
62 
63  void CreateChiSqrMaps(cmf::ChiSqrInfo & chiSqrInfo);
64 
65  void CheckMinimum1DChiSqr(cmf::ParamChiSqr & paramChiSqr,
66  double paramVal,
67  double chiSqr);
68 
69  double FindChiSqrForPoint(cmf::ChiSqrInfo & chiSqrInfo,
70  cmf::PointSpectra const& pointSpectra,
71  cmf::OscParamPoint const& oscPoint);
72 
73  void MakeAndStorePlots(std::string const& baseName,
74  cmf::ChiSqrInfo & chiSqrInfo);
75 
76  void MakeMedianContour();
77 
78  float fDeltaChiSqrToSkip; ///< how many units of chi^2 for a grid point before we ignore it when making the contours
79  float fFourFlavord24Val; ///< value of delta_24 to use
80  float fNDFHCExposureScale; ///< how much the ND FHC bins in the prediction libraries should be scaled by
81  float fNDRHCExposureScale; ///< how much the ND RHC bins in the prediction libraries should be scaled by
82  float fFDFHCExposureScale; ///< how much the FD FHC bins in the prediction libraries should be scaled by
83  float fFDRHCExposureScale; ///< how much the FD RHC bins in the prediction libraries should be scaled by
84  std::string fPredictionModule; ///< label for module containing the prediction point library
85  std::string fInputUniverseLabel; ///< label for module containing the prediction point library
86  cmf::ChiSqrCalc_t fChiSqrToUse; ///< which calculator to use
87  cmf::OscParm_t fXParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
88  cmf::OscParm_t fYParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
89  fhicl::ParameterSet fPlotConfig; ///< configuration for the plot utilities
90  bool fLogXAxis; ///< are parameters on a log scale on the X axis
91  bool fLogYAxis; ///< are parameters on a log scale on the Y axis
92  bool fScaleExposure; ///< should the exposure of the prediction libraries be scaled?
93  size_t fNumUniverses; ///< number of universes to use when making contours
94  size_t fStartUniverse; ///< which universe to start with in the fake universe file
95  cmf::FakeUniColl fFakeUniverses; ///< the input fake universes
96  std::vector<cmf::PointSpectra> fPointSpectraVec; ///< the point spectra for all the oscillation points
97  std::vector<cmf::ChiSqrInfo> fUniverseChiSqrInfo; ///< maps for the different universes
98 
99  std::map<cmf::DetBeamPair, float> fDBScaleMap; ///< how do we want to scale our exposure?
100 
101  };
102 
103  //......................................................................
107  , fPlotConfig(pset.get<fhicl::ParameterSet>("PlotUtilities"))
108  {
109  this->reconfigure(pset);
110  }
111 
112  //......................................................................
114  {
117  cmf::ParameterUtility::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse" ));
118  cmf::ChiSqrCalculator::Instance() ->InitializeCovarianceMatrix(pset.get<fhicl::ParameterSet>("ChiSqrCalculationParameters" ));
119 
120  fChiSqrToUse = cmf::StringToChiSqrType(pset.get<std::string>("ChiSqrCalcToUse"));
121 
122  fDeltaChiSqrToSkip = pset.get< float >("DeltaChiSqrToSkip", 25. );
123  fFourFlavord24Val = pset.get< float >("FourFlavord24Value", 0.0 );
124  fNDFHCExposureScale = pset.get< float >("NDFHCExposureScale", 1.0 );
125  fFDFHCExposureScale = pset.get< float >("FDFHCExposureScale", 1.0 );
126  fNDRHCExposureScale = pset.get< float >("NDRHCExposureScale", 1.0 );
127  fFDRHCExposureScale = pset.get< float >("FDRHCExposureScale", 1.0 );
128  fInputUniverseLabel = pset.get<std::string >("InputUniverseModule" );
129  fLogXAxis = pset.get< bool >("LogXAxis", false);
130  fLogYAxis = pset.get< bool >("LogYAxis", false);
131  fScaleExposure = pset.get< bool >("ScaleExposure", false);
132  fNumUniverses = pset.get<unsigned int>("NumberOfUniverses", 0 );
133  fPredictionModule = pset.get<std::string >("PredictionLibraryModule" );
134  fStartUniverse = pset.get<unsigned int>("StartUniverse", 0 );
135 
136  fDBScaleMap.emplace(std::make_pair(cmf::kNEARDET, cmf::kFHC), fNDFHCExposureScale);
137  fDBScaleMap.emplace(std::make_pair(cmf::kFARDET, cmf::kFHC), fFDFHCExposureScale);
138  fDBScaleMap.emplace(std::make_pair(cmf::kNEARDET, cmf::kRHC), fNDRHCExposureScale);
139  fDBScaleMap.emplace(std::make_pair(cmf::kFARDET, cmf::kRHC), fFDRHCExposureScale);
140 
141 
144  }
145 
146  //......................................................................
147  // Method to clear out the collections of data products after the
148  // writeResults method is called.
150  {}
151 
152  //......................................................................
153  // We need to read in several types of files - the best fit,
154  // prediction library and random universe files. The latter is what we
155  // want if we are making median sensitivities and also has the Asimov
156  // predictions for the fake data.
158  {
159 
161  if( r.getByLabel(fInputUniverseLabel, universes) &&
162  universes.isValid() ){
163 
164  fFakeUniverses.reserve(universes->size() + fFakeUniverses.size());
165  for(auto const& itr : *universes) fFakeUniverses.emplace_back(itr);
166  MF_LOG_VERBATIM("FitFeldmanCousinsPoint")
167  << "Valid handle to std::vector<cmf::FakeUniverse> "
168  << "objects found in "
170  << " vector size "
171  << fFakeUniverses.size();
172  r.removeCachedProduct(universes);
173  }
174 
176  if( r.getByLabel(fPredictionModule, predLib) &&
177  predLib.isValid() ){
178  MF_LOG_VERBATIM("CovarianceMatrixFitter")
179  << "Valid handle to std::vector<cmf::OscParamPoint> "
180  << "objects found in "
182  << " "
183  << " vector size "
184  << predLib->size();
185 
186  auto useTrig = fPlotConfig.get<bool>("UseTrig");
187  auto psvItr = fPointSpectraVec.begin();
188  float xVal;
189  float yVal;
190 
191  // maybe we should first figure out how many grid points we need in the space
192  // and how many entries we need for the hidden parameter combinations at each
193  // point and reserve that space
194  auto const& dbBinsVec = cmf::CovarianceBinUtility::Instance()->DetectorBeamBins();
195 
196  for (auto const &itrUnscaled : *predLib){
197 
199  itrUnscaled.ScaleSpectra(dbBinsVec, fDBScaleMap) :
200  itrUnscaled;
201 
202  MF_LOG_DEBUG("CMFContourFromLibrary")
203  << itr;
204 
205  auto tmpMap = itr.OscPoint();
206  if(tmpMap.count(cmf::kd24) > 0){
207  if(tmpMap.find(cmf::kd24)->second != fFourFlavord24Val){
208  MF_LOG_DEBUG("CMFContourFromLibrary")
209  << "don't use this point";
210 
211  continue;
212  }
213  }
214 
215  xVal = (useTrig && cmf::IsAngleParameter(fXParamEnum)) ? std::pow(std::sin(itr.OscPoint().find(fXParamEnum)->second), 2.) : itr.OscPoint().find(fXParamEnum)->second;
216  yVal = (useTrig && cmf::IsAngleParameter(fYParamEnum)) ? std::pow(std::sin(itr.OscPoint().find(fYParamEnum)->second), 2.) : itr.OscPoint().find(fYParamEnum)->second;
217 
218  auto gridPoint = cmf::GridPoint(xVal,
219  yVal,
220  fXParamEnum,
221  fYParamEnum);
222  psvItr = std::find(fPointSpectraVec.begin(),
223  fPointSpectraVec.end(),
224  PointSpectra(gridPoint));
225  if (psvItr == fPointSpectraVec.end()){
226  fPointSpectraVec.emplace_back(gridPoint);
227  fPointSpectraVec.back().oscPoints.emplace_back(itr);
228  }
229  else
230  psvItr->oscPoints.emplace_back(itr);
231  }
232 
233  // release space needed to read in the prediction library
234  r.removeCachedProduct(predLib);
235  }
236 
237  // for(auto const& itr : fPointSpectraVec)
238  // MF_LOG_VERBATIM("ContourFromLibrary")
239  // << itr.point
240  // << " has "
241  // << itr.oscPoints.size()
242  // << " spectra associated with it";
243 
244  }
245 
246  //......................................................................
248  cmf::PointSpectra const& pointSpectra,
249  cmf::OscParamPoint const& oscPoint)
250  {
251  // set the MC spectrum to the input Asimov spectrum to find the chi^2 for this universe
252  // compared to the known input values
254  false);
256  oscPoint.OscPoint());
257 
258  if(chiSqr < chiSqrInfo.minChiSqr){
259  chiSqrInfo.minChiSqr = chiSqr;
260  chiSqrInfo.bestFitSpectrum = cmf::Spectrum(oscPoint.PredictedSpectrum());
261  chiSqrInfo.bestFit = cmf::GridPoint(pointSpectra);
262  chiSqrInfo.bestFitVals = cmf::OscillationParameterMap(oscPoint.OscPoint());
263  }
264 
265  auto twoDItr = chiSqrInfo.twoDChiSqr.find(pointSpectra);
266  if(twoDItr == chiSqrInfo.twoDChiSqr.end())
267  chiSqrInfo.twoDChiSqr.emplace(pointSpectra, chiSqr);
268  else{
269  if(chiSqr < twoDItr->second) twoDItr->second = chiSqr;
270  }
271 
272  this->CheckMinimum1DChiSqr(chiSqrInfo.param1ChiSqr, pointSpectra.X(), chiSqr);
273  this->CheckMinimum1DChiSqr(chiSqrInfo.param2ChiSqr, pointSpectra.Y(), chiSqr);
274 
275  return chiSqr;
276  }
277 
278  //......................................................................
280  {
281 
283 
284  // we can save ourselves some time by first getting a basic idea
285  // of what the chi^2 map for the space looks like. We will loop
286  // over all the entries in fPointSpectraVec and get the chi^2
287  // for just the first, middle and last set of parameters for each
288  // grid point. Then we will loop over it again and ignore all points
289  // with a Delta chi^2 between the minimum from the basic stan and
290  // the current point of more than fDeltaChiSqrToSkip - ie very far
291  // outside the 3sigma confidence interval
292  for(auto & itr : fPointSpectraVec){
293  this->FindChiSqrForPoint(chiSqrInfo, itr, itr.oscPoints.front());
294  this->FindChiSqrForPoint(chiSqrInfo, itr, itr.oscPoints[itr.oscPoints.size() / 2]);
295  this->FindChiSqrForPoint(chiSqrInfo, itr, itr.oscPoints.back());
296  }
297 
298  double chiSqr;
299  double minChiSqrForPoint;
300 
301  // the vector of PointSpectra was sorted before this method
302  // was called. start by filling the maps with just the chi^2 values,
303  // we'll subtract the minimum off later.
304  for(auto & itr : fPointSpectraVec){
305 
306  // ignore points that are more than fDeltaChiSqrToSkip units of chi^2 away from the initial
307  // minimum as they are not likely to be near any of the confidence intervals we are interested in
308  if(chiSqrInfo.twoDChiSqr.find(itr)->second > chiSqrInfo.minChiSqr + fDeltaChiSqrToSkip)
309  continue;
310 
311  minChiSqrForPoint = cmf::kGarbageDouble;
312 
313  MF_LOG_DEBUG("ContourFromLibrary")
314  << " grid point: "
315  << itr;
316 
317  for(auto const& pitr : itr.oscPoints){
318 
319  chiSqr = this->FindChiSqrForPoint(chiSqrInfo, itr, pitr);
320 
321  // if(itr == *(fPointSpectraVec.begin()))
322  // MF_LOG_VERBATIM("ContourFromLibrary")
323  // << "check using first point in space: "
324  // << pitr
325  // << " chi^2 = "
326  // << chiSqr;
327 
328  MF_LOG_VERBATIM("ContourFromLibrary")
329  << "check point in space: "
330  << pitr
331  << " chi^2 = "
332  << chiSqr;
333 
334  if(chiSqr < minChiSqrForPoint){
335  minChiSqrForPoint = chiSqr;
336  itr.bestFitForPoint = cmf::OscillationParameterMap(pitr.OscPoint());
337  }
338  } // end loop over the OscParamPoints for this location in the space
339 
340  MF_LOG_DEBUG("ContourFromLibrary")
341  << "Point: "
342  << itr
343  << " has min chi^2 "
344  << minChiSqrForPoint
345  << " min chi^2 for space: "
346  << chiSqrInfo.minChiSqr;
347 
348  } // end loop over the point spectra
349 
350  // now subtract off the minimum chi^2 value from the map entries
351  for(auto & itr : chiSqrInfo.param1ChiSqr) itr.second -= chiSqrInfo.minChiSqr;
352  for(auto & itr : chiSqrInfo.param2ChiSqr) itr.second -= chiSqrInfo.minChiSqr;
353  for(auto & itr : chiSqrInfo.twoDChiSqr ) itr.second -= chiSqrInfo.minChiSqr;
354  }
355 
356  //----------------------------------------------------------------------------
358  double paramVal,
359  double chiSqr)
360  {
361  if(paramChiSqr.count(paramVal)){
362  paramChiSqr[paramVal] = std::min(chiSqr, paramChiSqr[paramVal]);
363  }
364  else{
365  paramChiSqr.emplace(paramVal, chiSqr);
366  }
367  }
368 
369  //----------------------------------------------------------------------------
371  cmf::Spectrum const& dataSpectrum)
372  {
373  // make the chiSqr maps
374  cmf::ChiSqrInfo chiSqrInfo;
375 
376  chiSqrInfo.dataSpectrum = dataSpectrum;
377 
378  // here is where we would pick out which universe we want to
379  // make a contour for or loop over multiple universes
380  this->CreateChiSqrMaps(chiSqrInfo);
381 
382  std::string bestFitOutput(baseName + " best fit found at");
383  for(auto const& itr : chiSqrInfo.bestFitVals)
384  bestFitOutput += " " + cmf::cOscParams_Strings[itr.first] + " " + std::to_string(itr.second);
385 
386  bestFitOutput += " with chi^{2} = " + std::to_string(chiSqrInfo.minChiSqr);
387 
388  MF_LOG_VERBATIM("ContourFromLibrary")
389  << bestFitOutput;
390 
391  this->MakeAndStorePlots(baseName,
392  chiSqrInfo);
393 
394  // store the ChiSqrInfo for fake universes, but not the Asimov universe
395  if(baseName.find("Asimov") == std::string::npos) fUniverseChiSqrInfo.emplace_back(chiSqrInfo);
396  }
397 
398  //----------------------------------------------------------------------------
400  {
401  // make sure we aren't taking up too much space
402  fPointSpectraVec.shrink_to_fit();
403  fFakeUniverses.shrink_to_fit();
404 
405  // do the same for the osc points in the vector
406  for(auto & itr : fPointSpectraVec) itr.oscPoints.shrink_to_fit();
407 
408  MF_LOG_VERBATIM("ContourFromLibrary")
409  << " There are "
410  << fPointSpectraVec.size()
411  << " library points in the space with "
412  << fPointSpectraVec.front().oscPoints.size()
413  << " predictions at each point\n There are "
414  << fFakeUniverses.size()
415  << " fake universes";
416 
417  // check that the number of bins in the predicted spectra and fake
418  // universes match the expectation for the job
419  size_t totConfigBins = cmf::CovarianceBinUtility::Instance()->TotalBins();
420  //if(fPointSpectraVec.front().oscPoints.front().PredictedSpectrum().size() != totConfigBins)
421  // throw cet::exception("ContourFromLibrary")
422  // << "The predicted spectra have "
423  // << fPointSpectraVec.front().oscPoints.front().PredictedSpectrum().size()
424  // << " bins while the configuration expects "
425  // << totConfigBins;
426 
427  //if(fFakeUniverses.front().AsimovSpectrum().size() != totConfigBins)
428  // throw cet::exception("ContourFromLibrary")
429  // << "The predicted spectra have "
430  // << fFakeUniverses.front().AsimovSpectrum().size()
431  // << " bins while the configuration expects "
432  // << totConfigBins;
433 
434  // if we ever wanted to downsample the number of spectra we use
435  // from the library, the call to do so should come here
436 
437  // now sort the universe point results in the library points vector
438  std::sort(fPointSpectraVec.begin(), fPointSpectraVec.end());
439 
440  // find the extrema in the space
441  std::pair<float, float> parXExtrema = std::make_pair(fPointSpectraVec.front().X(),
442  fPointSpectraVec.back() .X());
443  std::pair<float, float> parYExtrema = std::make_pair(fPointSpectraVec.front().Y(),
444  fPointSpectraVec.back() .Y());
445 
446  MF_LOG_DEBUG("ContourFromLibrary")
448  << " range is "
449  << parXExtrema.first
450  << " - "
451  << parXExtrema.second
452  << " "
454  << " range is "
455  << parYExtrema.first
456  << " - "
457  << parYExtrema.second;
458 
459  // figure out how many different points we have for each parameter.
460  std::set<float> xPoints;
461  std::set<float> yPoints;
462  for(auto & itr : fPointSpectraVec){
463 
464  xPoints.insert(itr.X());
465  yPoints.insert(itr.Y());
466 
467  // sort the vector of OscParamPoints
468  std::sort(itr.oscPoints.begin(), itr.oscPoints.end());
469 
470  MF_LOG_DEBUG("ContourFromLibrary")
471  << itr;
472  }
473 
474  // for(auto const& oitr : fPointSpectraVec.front().oscPoints)
475  // MF_LOG_VERBATIM("ContourFromLibrary")
476  // << "check osc point ordering: "
477  // << oitr;
478 
479  MF_LOG_DEBUG("ContourFromLibrary")
480  << "There are "
481  << xPoints.size()
482  << " points in "
484  << " and "
485  << yPoints.size()
486  << " points in "
488 
490  fXParamEnum,
491  fYParamEnum,
492  xPoints.size() - 1,
493  yPoints.size() - 1,
494  parXExtrema,
495  parYExtrema,
496  fLogXAxis,
497  fLogYAxis);
498 
499  this->ContourForUniverse("Asimov", fFakeUniverses.front().AsimovSpectrum());
500 
501  // make sure we have a clean slate for storing the ChiSqrInfo for the fake universes
502  // in ContourForUniverse
503  fUniverseChiSqrInfo.clear();
504 
505  // if NumUniverses is still set to 0, then the user did not specify a number so
506  // loop over all the fake universes
507  if(fNumUniverses < 1) fNumUniverses = fFakeUniverses.size();
508 
509  MF_LOG_VERBATIM("ContourFromLibrary")
510  << "Find contours for "
511  << fNumUniverses
512  << " universes starting with number "
513  << fStartUniverse;
514 
515  for(size_t u = fStartUniverse; u < fStartUniverse + fNumUniverses; ++u){
516  this->ContourForUniverse("Universe_" + std::to_string(u), fFakeUniverses[u].PoissonSpectrum());
517  }
518 
519  // finally make the median contours for the fake universes
520  this->MakeMedianContour();
521  }
522 
523  //----------------------------------------------------------------------------
525  cmf::ChiSqrInfo & chiSqrInfo)
526  {
528 
529  // make a TFileDirectory to store the histograms for this calculation
530  art::TFileDirectory curDir = tfs->mkdir(baseName);
531 
532  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param1ChiSqr, true);
533  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param2ChiSqr, false);
534  cmf::PlotUtilities::Instance()->Make2DContours(curDir, baseName, chiSqrInfo.twoDChiSqr, chiSqrInfo.bestFit, true);
535 
536  if(baseName.find("Median") != std::string::npos){
537  // let's make some heat map plots.
538  cmf::PlotUtilities::Instance()->MakeCLHeatMaps(curDir, "RandomUniverses", fUniverseChiSqrInfo);
539 
540  // we are done here if looking at plots for median contours as best fits
541  // don't mean anything in that context
542  return;
543  }
544  // 2D histograms for the best fit hidden parameter values at each point in space
545  std::map<cmf::OscParm_t, cmf::PointMap> pointMaps;
546 
547  for(auto const& itr : fPointSpectraVec){
548  for(auto const& opItr : itr.bestFitForPoint){
549  if(opItr.first == fXParamEnum ||
550  opItr.second == fYParamEnum ) continue;
551  if(cmf::IsAngleParameter(opItr.first))
552  pointMaps[opItr.first][itr] = std::pow(std::sin(opItr.second), 2.);
553  else
554  pointMaps[opItr.first][itr] = opItr.second * cmf::cOscParams_Scales[opItr.first];
555  }
556  }
557 
558  for(auto const& itr : pointMaps)
560  (baseName + cmf::cOscParams_Strings[itr.first]).c_str(),
561  itr.second);
562 
563  auto *fitResultGr = curDir.makeAndRegister<TGraph2D>((baseName + "FitResult").c_str(),
564  (baseName + "FitResult").c_str(),
565  1);
566 
567  fitResultGr->SetPoint(0,
568  chiSqrInfo.bestFit.X() * cmf::cOscParams_Scales[chiSqrInfo.bestFit.XPar()],
569  chiSqrInfo.bestFit.Y() * cmf::cOscParams_Scales[chiSqrInfo.bestFit.YPar()],
570  chiSqrInfo.minChiSqr);
571  fitResultGr->SetMarkerColor(kBlue);
572  //fitResultGr->GetHistogram()->SetXTitle(cmf::cOscParams_Strings_Latex[bestFit.XPar()].c_str());
573  //fitResultGr->GetHistogram()->SetYTitle(cmf::cOscParams_Strings_Latex[bestFit.YPar()].c_str());
574  //fitResultGr->GetHistogram()->SetZTitle("Best Fit #chi^{2}");
575 
577  baseName + "DataSpectrum",
578  ";Bin;Events",
579  chiSqrInfo.dataSpectrum);
580 
582  baseName + "BestFitSpectrum",
583  ";Bin;Events",
584  chiSqrInfo.bestFitSpectrum);
585 
586  }
587 
588  //----------------------------------------------------------------------------
590  {
591  // first off we need to sort the the Delta chi^2 values for each point in the space
592  std::map<cmf::GridPoint, std::vector<double>> pointDeltaChiSqrs;
593  for(auto const& csItr : fUniverseChiSqrInfo){
594  for(auto const& twoDItr : csItr.twoDChiSqr){
595  pointDeltaChiSqrs[twoDItr.first].emplace_back(twoDItr.second);
596  }
597  }
598 
599  cmf::ChiSqrInfo medianChiSqrInfo;
600  double chiSqr;
601  // now sort the Delta chi^2 values
602  for(auto & itr : pointDeltaChiSqrs){
603  std::sort(itr.second.begin(), itr.second.end());
604 
605  chiSqr = itr.second[itr.second.size() / 2];
606 
607  this->CheckMinimum1DChiSqr(medianChiSqrInfo.param1ChiSqr, itr.first.X(), chiSqr);
608  this->CheckMinimum1DChiSqr(medianChiSqrInfo.param2ChiSqr, itr.first.Y(), chiSqr);
609  medianChiSqrInfo.twoDChiSqr.emplace(itr.first, chiSqr);
610  } // end loop to create the medianChiSqrInfo
611 
612  // now make the contours
613  this->MakeAndStorePlots("Median",
614  medianChiSqrInfo);
615 
616  }
617 
618  //----------------------------------------------------------------------------
620  {
621  }
622 
624 
625 } // end namespace
cmf::OscillationParameterMap OscPoint() const
enum cmf::osc_params OscParm_t
void InitializeCovarianceMatrix(fhicl::ParameterSet const &pset)
std::vector< cmf::FakeUniverse > FakeUniColl
Definition: FakeUniverse.h:14
cmf::PointMap twoDChiSqr
double FindChiSqrForPoint(cmf::ChiSqrInfo &chiSqrInfo, cmf::PointSpectra const &pointSpectra, cmf::OscParamPoint const &oscPoint)
cmf::ParamChiSqr param2ChiSqr
cmf::OscillationParameterMap bestFitVals
float fFDRHCExposureScale
how much the FD RHC bins in the prediction libraries should be scaled by
void Initialize(fhicl::ParameterSet const &pset)
~ContourFromLibrary() override=default
cmf::OscParm_t fXParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
std::vector< cmf::PointSpectra > fPointSpectraVec
the point spectra for all the oscillation points
const char * p
Definition: xmltok.h:285
static const double kGarbageDouble
Definition: Constants.h:18
constexpr T pow(T x)
Definition: pow.h:72
std::map< double, double > ParamChiSqr
Definition: Constants.h:747
void CreateChiSqrMaps(cmf::ChiSqrInfo &chiSqrInfo)
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:746
cmf::GridPoint bestFit
std::vector< cmf::DetBeamBins > const & DetectorBeamBins()
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
static bool IsAngleParameter(cmf::OscParm_t const &par)
Definition: StaticFuncs.h:354
cmf::Spectrum dataSpectrum
void Make1DPlot(art::TFileDirectory &tfd, std::string const &chiSqrName, std::map< double, double > const &paramChiSqr, bool xPar)
cmf::OscParm_t fYParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
cmf::ParamChiSqr param1ChiSqr
void Make2DContours(art::TFileDirectory &tfd, std::string const &chiSqrName, cmf::PointMap const &twoDChiSqr, cmf::GridPoint const &bestFit, bool useBestFit=true)
cmf::OscParamPoint ScaleSpectra(std::vector< cmf::DetBeamBins > const &dbBinsVec, std::map< cmf::DetBeamPair, float > const &dBScaleMap) const
void readResults(art::Results const &r) override
void CheckMinimum1DChiSqr(cmf::ParamChiSqr &paramChiSqr, double paramVal, double chiSqr)
bool isValid() const
Definition: Handle.h:183
float fNDFHCExposureScale
how much the ND FHC bins in the prediction libraries should be scaled by
static ParameterUtility * Instance()
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
float fFDFHCExposureScale
how much the FD FHC bins in the prediction libraries should be scaled by
void reconfigure(fhicl::ParameterSet const &p)
std::map< cmf::OscParm_t, float > OscillationParameterMap
Definition: Constants.h:748
std::string fPredictionModule
label for module containing the prediction point library
cmf::FakeUniColl fFakeUniverses
the input fake universes
#define DEFINE_ART_RESULTS_PLUGIN(klass)
void ContourForUniverse(std::string const &baseName, cmf::Spectrum const &dataSpectrum)
ContourFromLibrary(fhicl::ParameterSet const &pset)
T get(std::string const &key) const
Definition: ParameterSet.h:231
cmf::Spectrum const & PredictedSpectrum() const
std::map< cmf::DetBeamPair, float > fDBScaleMap
how do we want to scale our exposure?
cmf::ChiSqrCalc_t fChiSqrToUse
which calculator to use
void Initialize(fhicl::ParameterSet const &pset)
std::vector< cmf::ChiSqrInfo > fUniverseChiSqrInfo
maps for the different universes
void Initialize(fhicl::ParameterSet const &plottingPars)
bool removeCachedProduct(Handle< PROD > &) const
Definition: DataViewImpl.h:971
float const & X() const
static cmf::ChiSqrCalculator * Instance()
double ChiSqr(cmf::ChiSqrCalc_t const &chiSqrCalc, std::vector< std::string > const &pars, const double *vals)
size_t fStartUniverse
which universe to start with in the fake universe file
float fDeltaChiSqrToSkip
how many units of chi^2 for a grid point before we ignore it when making the contours ...
size_t TotalBins(bool allSels=false)
void MakeAndStorePlots(std::string const &baseName, cmf::ChiSqrInfo &chiSqrInfo)
static cmf::ChiSqrCalc_t StringToChiSqrType(std::string const &str)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
size_t fNumUniverses
number of universes to use when making contours
T sin(T number)
Definition: d0nt_math.hpp:132
bool fLogYAxis
are parameters on a log scale on the Y axis
cmf::OscParm_t const & XPar() const
cmf::Spectrum bestFitSpectrum
#define MF_LOG_VERBATIM(category)
bool fScaleExposure
should the exposure of the prediction libraries be scaled?
fhicl::ParameterSet fPlotConfig
configuration for the plot utilities
#define MF_LOG_DEBUG(id)
float const & Y() const
void writeResults(art::Results &r) override
static PlotUtilities * Instance()
float fNDRHCExposureScale
how much the ND RHC bins in the prediction libraries should be scaled by
TRandom3 r(0)
const std::vector< double > cOscParams_Scales({1., 1., 1.e5, 1.e3, 1., 1., 1., 1./3.14159, 1., 1., 1., 1., 1./3.14159, 1., 1., 1., 1., 1., 1., 1., 1., 1.})
void MakeCLHeatMaps(art::TFileDirectory &tfd, std::string const &baseName, std::vector< cmf::ChiSqrInfo > const &chiSqrInfoVec)
T min(const caf::Proxy< T > &a, T b)
void Initialize(fhicl::ParameterSet const &pset)
enum BeamMode kBlue
cmf::OscParm_t const & YPar() const
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
static cmf::OscParm_t StringToOscParmType(std::string const &str)
float fFourFlavord24Val
value of delta_24 to use
void Make2DHiddenParHistogram(art::TFileDirectory &tfd, std::string const &baseName, cmf::PointMap const &parVals)
std::string fInputUniverseLabel
label for module containing the prediction point library
enum cmf::chisqr_type ChiSqrCalc_t
void SetSpectrum(cmf::Spectrum const &spectrum, bool isData)
static CovarianceBinUtility * Instance()
bool fLogXAxis
are parameters on a log scale on the X axis
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