CMFContourFromGrid_plugin.cc
Go to the documentation of this file.
1 //
2 // Plugin to create contours from a set of PointBestFits.
3 // Pass the art job a list of input PointBestFit files for the space
4 //
5 // Created by Brian Rebel on 3/23/21.
6 //
7 
8 #include <string>
9 #include <set>
10 #include <map>
11 #include <memory>
12 #include <sstream>
13 
14 //#include "TStopwatch.h"
15 #include "TMath.h"
16 #include "TGraph2D.h"
17 
18 // Framework includes
24 #include "art_root_io/TFileService.h"
25 #include "art_root_io/TFileDirectory.h"
26 #include "fhiclcpp/ParameterSet.h"
28 #include "cetlib_except/exception.h"
29 
30 // NOvA includes
39 
40 namespace cmf{
41 
43  public:
44  explicit ContourFromGrid(fhicl::ParameterSet const& pset);
45 
46  ~ContourFromGrid() override = default;
47 
48  void readResults(art::Results const& r) override;
49 
50  void writeResults(art::Results& r) override;
51 
52  void reconfigure(fhicl::ParameterSet const& p);
53 
54  void clear() override;
55  void endJob() override;
56 
57  private:
58 
59  void ContourForUniverse(std::string const& baseName,
60  std::vector<cmf::PointBestFit> const& randomColl,
61  std::vector<cmf::PointBestFit> const& asimovColl,
62  cmf::ChiSqrInfo & randomChiSqr,
63  cmf::ChiSqrInfo & asimovChiSqr);
64 
65  void CreateChiSqrMaps(std::string const& baseName,
66  std::vector<cmf::PointBestFit> const& pbfColl,
67  cmf::ChiSqrInfo & chiSqrInfo);
68 
69  void CheckMinimum1DChiSqr(cmf::ParamChiSqr & paramChiSqr,
70  double paramVal,
71  double chiSqr);
72 
73  void MakeAndStorePlots(std::string const& baseName,
74  std::vector<cmf::PointBestFit> const& pbfColl,
75  cmf::ChiSqrInfo & chiSqrInfo);
76 
77  void MakeAndStoreMedianPlots(std::string const& baseName,
78  cmf::ChiSqrInfo & chiSqrInfo);
79  void MakeMedianContour();
80 
81  std::string fPBFModule; ///< label for module containing the best fits for the grid points
82  cmf::OscParm_t fXParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
83  cmf::OscParm_t fYParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
84  fhicl::ParameterSet fPlotConfig; ///< configuration for the plot utilities
85  bool fLogXAxis; ///< are parameters on a log scale on the X axis
86  bool fLogYAxis; ///< are parameters on a log scale on the Y axis
87  std::vector<cmf::ChiSqrInfo> fRandomUniChiSqrInfo; ///< maps for the different universes
88  std::vector<cmf::ChiSqrInfo> fAsimovUniChiSqrInfo; ///< maps for the different universes
89  std::vector<cmf::UniversePointBestFits> fUniverseFits; ///< collection of best fit info for universes
90 
91  };
92 
93  //......................................................................
97  , fPlotConfig(pset.get<fhicl::ParameterSet>("PlotUtilities"))
98  {
99  this->reconfigure(pset);
100  }
101 
102  //......................................................................
104  {
107  cmf::ParameterUtility::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse" ));
108 
109  fLogXAxis = pset.get< bool >("LogXAxis", false);
110  fLogYAxis = pset.get< bool >("LogYAxis", false);
111  fPBFModule = pset.get<std::string >("PBFModule" );
112 
115  }
116 
117  //......................................................................
118  // Method to clear out the collections of data products after the
119  // writeResults method is called.
121  {}
122 
123  //......................................................................
124  // We need to read in several types of files - the best fit,
125  // prediction library and random universe files. The latter is what we
126  // want if we are making median sensitivities and also has the Asimov
127  // predictions for the fake data.
129  {
130 
132  if( r.getByLabel(fPBFModule, universes) &&
133  universes.isValid() ){
134 
135  if(fUniverseFits.size() < universes->size())
136  fUniverseFits.resize(universes->size(), cmf::UniversePointBestFits());
137 
138  for(size_t u = 0; u < universes->size(); ++u){
139  fUniverseFits[u].AddPointsToUniverse((*universes)[u].RandomUniversePBF(),
140  (*universes)[u].AsimovUniversePBF());
141  }
142 
143  MF_LOG_VERBATIM("FitFeldmanCousinsPoint")
144  << "Valid handle to std::vector<cmf::UniversePointBestFits> "
145  << "objects found in "
146  << fPBFModule
147  << " vector size "
148  << fUniverseFits.size();
149  r.removeCachedProduct(universes);
150  }
151 
152  }
153 
154  //......................................................................
156  std::vector<cmf::PointBestFit> const& pbfColl,
157  cmf::ChiSqrInfo & chiSqrInfo)
158  {
159  // start by filling the maps with just the chi^2 values,
160  // we'll subtract the minimum off later.
161  chiSqrInfo.minChiSqr = cmf::kGarbageDouble;
162 
163  for(auto const& itr : pbfColl){
164 
165  chiSqrInfo.twoDChiSqr.emplace(itr, itr.MinChiSqr());
166  this->CheckMinimum1DChiSqr(chiSqrInfo.param1ChiSqr,
167  itr.X(),
168  itr.MinChiSqr());
169  this->CheckMinimum1DChiSqr(chiSqrInfo.param2ChiSqr,
170  itr.Y(),
171  itr.MinChiSqr());
172 
173  if(itr.MinChiSqr() < chiSqrInfo.minChiSqr){
174  chiSqrInfo.minChiSqr = itr.MinChiSqr();
175  chiSqrInfo.bestFitVals = cmf::OscillationParameterMap(itr.BestFitOscPars());
176  chiSqrInfo.bestFitSpectrum = itr.BestFitSpectrum();
177  }
178 
179  MF_LOG_VERBATIM("ContourFromGrid")
180  << "Point: "
181  << itr
182  << " has min chi^2 "
183  << itr.MinChiSqr()
184  << " min chi^2 for space: "
185  << chiSqrInfo.minChiSqr;
186 
187  } // end loop over the point best fits
188 
189  std::string bestFitOutput(baseName + " universe best fit found at");
190  for(auto const& itr : chiSqrInfo.bestFitVals)
191  bestFitOutput += " " + cmf::cOscParams_Strings[itr.first] + " " + std::to_string(itr.second);
192 
193  bestFitOutput += " with chi^{2} = " + std::to_string(chiSqrInfo.minChiSqr);
194 
195  MF_LOG_VERBATIM("ContourFromGrid")
196  << bestFitOutput;
197 
198  // now subtract off the minimum chi^2 value from the map entries
199  for(auto & itr : chiSqrInfo.param1ChiSqr) itr.second -= chiSqrInfo.minChiSqr;
200  for(auto & itr : chiSqrInfo.param2ChiSqr) itr.second -= chiSqrInfo.minChiSqr;
201  for(auto & itr : chiSqrInfo.twoDChiSqr ) itr.second -= chiSqrInfo.minChiSqr;
202  }
203 
204  //----------------------------------------------------------------------------
206  double paramVal,
207  double chiSqr)
208  {
209  if(paramChiSqr.count(paramVal)){
210  paramChiSqr[paramVal] = std::min(chiSqr, paramChiSqr[paramVal]);
211  }
212  else{
213  paramChiSqr.emplace(paramVal, chiSqr);
214  }
215  }
216 
217  //----------------------------------------------------------------------------
219  std::vector<cmf::PointBestFit> const& randomColl,
220  std::vector<cmf::PointBestFit> const& asimovColl,
221  cmf::ChiSqrInfo & randomChiSqr,
222  cmf::ChiSqrInfo & asimovChiSqr)
223  {
224 
225  randomChiSqr.dataSpectrum = randomColl.front().DataSpectrum();
226  std::string base = baseName + "Random";
227  this->CreateChiSqrMaps(base,
228  randomColl,
229  randomChiSqr);
230  this->MakeAndStorePlots(base,
231  randomColl,
232  randomChiSqr);
233 
234  asimovChiSqr.dataSpectrum = asimovColl.front().DataSpectrum();
235  base = baseName + "Asimov";
236  this->CreateChiSqrMaps(base,
237  asimovColl,
238  asimovChiSqr);
239  this->MakeAndStorePlots(base,
240  asimovColl,
241  asimovChiSqr);
242 
243  }
244 
245  //----------------------------------------------------------------------------
247  {
248  // make sure we aren't taking up too much space
249  fUniverseFits.shrink_to_fit();
250 
251  MF_LOG_VERBATIM("ContourFromGrid")
252  << " There are "
253  << fUniverseFits.size()
254  << " universes with "
255  << fUniverseFits.front().RandomUniversePBF().size()
256  << " points in the space";
257 
258  // now sort the points for each universe
259  for(auto & itr : fUniverseFits) itr.SortPoints();
260 
261  // find the extrema in the space
262  std::pair<float, float> parXExtrema = std::make_pair(fUniverseFits.front().RandomUniversePBF().front().X(),
263  fUniverseFits.front().RandomUniversePBF().back() .X());
264  std::pair<float, float> parYExtrema = std::make_pair(fUniverseFits.front().RandomUniversePBF().front().Y(),
265  fUniverseFits.front().RandomUniversePBF().back() .Y());
266 
267  MF_LOG_DEBUG("ContourFromGrid")
269  << " range is "
270  << parXExtrema.first
271  << " - "
272  << parXExtrema.second
273  << " "
275  << " range is "
276  << parYExtrema.first
277  << " - "
278  << parYExtrema.second;
279 
280  // figure out how many different points we have for each parameter.
281  std::set<float> xPoints;
282  std::set<float> yPoints;
283 
284  for(auto & itr : fUniverseFits.front().RandomUniversePBF()){
285  xPoints.insert(itr.X());
286  yPoints.insert(itr.Y());
287  }
288 
289  // for(auto const& oitr : fPointSpectraVec.front().oscPoints)
290  // MF_LOG_VERBATIM("ContourFromGrid")
291  // << "check osc point ordering: "
292  // << oitr;
293 
294  MF_LOG_DEBUG("ContourFromGrid")
295  << "There are "
296  << xPoints.size()
297  << " points in "
299  << " and "
300  << yPoints.size()
301  << " points in "
303 
305  fXParamEnum,
306  fYParamEnum,
307  xPoints.size() - 1,
308  yPoints.size() - 1,
309  parXExtrema,
310  parYExtrema,
311  fLogXAxis,
312  fLogYAxis);
313 
314  // make sure we have a clean slate for storing the ChiSqrInfo for the fake universes
315  // in ContourForUniverse
316  fRandomUniChiSqrInfo.resize(fUniverseFits.size());
317  fAsimovUniChiSqrInfo.resize(fUniverseFits.size());
318 
319  for(size_t u = 0; u < fUniverseFits.size(); ++u){
320  this->ContourForUniverse("Universe_" + std::to_string(u),
321  fUniverseFits[u].RandomUniversePBF(),
322  fUniverseFits[u].AsimovUniversePBF(),
325  }
326 
327  // finally make the median contours for the fake universes
328  this->MakeMedianContour();
329  }
330 
331  //----------------------------------------------------------------------------
333  std::vector<cmf::PointBestFit> const& pbfColl,
334  cmf::ChiSqrInfo & chiSqrInfo)
335  {
337 
338  // make a TFileDirectory to store the histograms for this calculation
339  art::TFileDirectory curDir = tfs->mkdir(baseName);
340 
341  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param1ChiSqr, true);
342  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param2ChiSqr, false);
343  cmf::PlotUtilities::Instance()->Make2DContours(curDir, baseName, chiSqrInfo.twoDChiSqr, chiSqrInfo.bestFit, true);
344 
345  if(baseName.find("Median") != std::string::npos){
346  // let's make some heat map plots.
348 
349  // we are done here if looking at plots for median contours as best fits
350  // don't mean anything in that context
351  return;
352  }
353 
354  // 2D histograms for the best fit hidden parameter values at each point in space
355  std::map<cmf::OscParm_t, cmf::PointMap> pointMaps;
356 
357  for(auto const& itr : pbfColl){
358  for(auto const& opItr : itr.BestFitOscPars()){
359  if(opItr.first == fXParamEnum ||
360  opItr.second == fYParamEnum ) continue;
361  if(cmf::IsAngleParameter(opItr.first))
362  pointMaps[opItr.first][itr] = std::pow(std::sin(opItr.second), 2.);
363  else
364  pointMaps[opItr.first][itr] = opItr.second * cmf::cOscParams_Scales[opItr.first];
365  }
366  }
367 
368  for(auto const& itr : pointMaps)
370  (baseName + cmf::cOscParams_Strings[itr.first]).c_str(),
371  itr.second);
372 
373  auto *fitResultGr = curDir.makeAndRegister<TGraph2D>((baseName + "FitResult").c_str(),
374  (baseName + "FitResult").c_str(),
375  1);
376 
377  fitResultGr->SetPoint(0,
378  chiSqrInfo.bestFit.X() * cmf::cOscParams_Scales[chiSqrInfo.bestFit.XPar()],
379  chiSqrInfo.bestFit.Y() * cmf::cOscParams_Scales[chiSqrInfo.bestFit.YPar()],
380  chiSqrInfo.minChiSqr);
381  fitResultGr->SetMarkerColor(kBlue);
382  //fitResultGr->GetHistogram()->SetXTitle(cmf::cOscParams_Strings_Latex[bestFit.XPar()].c_str());
383  //fitResultGr->GetHistogram()->SetYTitle(cmf::cOscParams_Strings_Latex[bestFit.YPar()].c_str());
384  //fitResultGr->GetHistogram()->SetZTitle("Best Fit #chi^{2}");
385 
387  baseName + "DataSpectrum",
388  ";Bin;Events",
389  chiSqrInfo.dataSpectrum);
390 
392  baseName + "BestFitSpectrum",
393  ";Bin;Events",
394  chiSqrInfo.bestFitSpectrum);
395 
396  }
397 
398  //----------------------------------------------------------------------------
400  cmf::ChiSqrInfo & chiSqrInfo)
401  {
403 
404  // make a TFileDirectory to store the histograms for this calculation
405  art::TFileDirectory curDir = tfs->mkdir(baseName);
406 
407  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param1ChiSqr, true);
408  cmf::PlotUtilities::Instance()->Make1DPlot (curDir, baseName, chiSqrInfo.param2ChiSqr, false);
409  cmf::PlotUtilities::Instance()->Make2DContours(curDir, baseName, chiSqrInfo.twoDChiSqr, chiSqrInfo.bestFit, true);
410 
411  // let's make some heat map plots.
413 
414  }
415 
416  //----------------------------------------------------------------------------
418  {
419  // first off we need to sort the the Delta chi^2 values for each point in the space
420  std::map<cmf::GridPoint, std::vector<double>> pointDeltaChiSqrs;
421  for(auto const& csItr : fRandomUniChiSqrInfo){
422  for(auto const& twoDItr : csItr.twoDChiSqr){
423  pointDeltaChiSqrs[twoDItr.first].emplace_back(twoDItr.second);
424  }
425  }
426 
427  cmf::ChiSqrInfo medianChiSqrInfo;
428  double chiSqr;
429  // now sort the Delta chi^2 values
430  for(auto & itr : pointDeltaChiSqrs){
431  std::sort(itr.second.begin(), itr.second.end());
432 
433  chiSqr = itr.second[itr.second.size() / 2];
434 
435  this->CheckMinimum1DChiSqr(medianChiSqrInfo.param1ChiSqr, itr.first.X(), chiSqr);
436  this->CheckMinimum1DChiSqr(medianChiSqrInfo.param2ChiSqr, itr.first.Y(), chiSqr);
437  medianChiSqrInfo.twoDChiSqr.emplace(itr.first, chiSqr);
438  } // end loop to create the medianChiSqrInfo
439 
440  this->MakeAndStoreMedianPlots("Median",
441  medianChiSqrInfo);
442 
443  }
444 
445  //----------------------------------------------------------------------------
447  {
448  }
449 
451 
452 } // end namespace
void CheckMinimum1DChiSqr(cmf::ParamChiSqr &paramChiSqr, double paramVal, double chiSqr)
enum cmf::osc_params OscParm_t
cmf::PointMap twoDChiSqr
~ContourFromGrid() override=default
cmf::ParamChiSqr param2ChiSqr
cmf::OscillationParameterMap bestFitVals
void Initialize(fhicl::ParameterSet const &pset)
std::vector< cmf::UniversePointBestFits > fUniverseFits
collection of best fit info for universes
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
static SelectionUtility * Instance()
cmf::GridPoint bestFit
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
ContourFromGrid(fhicl::ParameterSet const &pset)
std::string fPBFModule
label for module containing the best fits for the grid points
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
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
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)
std::vector< cmf::ChiSqrInfo > fAsimovUniChiSqrInfo
maps for the different universes
cmf::ParamChiSqr param1ChiSqr
void readResults(art::Results const &r) override
void Make2DContours(art::TFileDirectory &tfd, std::string const &chiSqrName, cmf::PointMap const &twoDChiSqr, cmf::GridPoint const &bestFit, bool useBestFit=true)
bool fLogXAxis
are parameters on a log scale on the X axis
bool isValid() const
Definition: Handle.h:183
static ParameterUtility * Instance()
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
bool fLogYAxis
are parameters on a log scale on the Y axis
std::map< cmf::OscParm_t, float > OscillationParameterMap
Definition: Constants.h:748
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
void Initialize(fhicl::ParameterSet const &pset)
cmf::OscParm_t fXParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
void CreateChiSqrMaps(std::string const &baseName, std::vector< cmf::PointBestFit > const &pbfColl, cmf::ChiSqrInfo &chiSqrInfo)
void Initialize(fhicl::ParameterSet const &plottingPars)
void writeResults(art::Results &r) override
bool removeCachedProduct(Handle< PROD > &) const
Definition: DataViewImpl.h:971
float const & X() const
cmf::OscParm_t fYParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
T sin(T number)
Definition: d0nt_math.hpp:132
cmf::OscParm_t const & XPar() const
cmf::Spectrum bestFitSpectrum
#define MF_LOG_VERBATIM(category)
void MakeAndStorePlots(std::string const &baseName, std::vector< cmf::PointBestFit > const &pbfColl, cmf::ChiSqrInfo &chiSqrInfo)
fhicl::ParameterSet fPlotConfig
configuration for the plot utilities
#define MF_LOG_DEBUG(id)
float const & Y() const
static PlotUtilities * Instance()
std::vector< cmf::ChiSqrInfo > fRandomUniChiSqrInfo
maps for the different universes
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)
void ContourForUniverse(std::string const &baseName, std::vector< cmf::PointBestFit > const &randomColl, std::vector< cmf::PointBestFit > const &asimovColl, cmf::ChiSqrInfo &randomChiSqr, cmf::ChiSqrInfo &asimovChiSqr)
void reconfigure(fhicl::ParameterSet const &p)
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)
void MakeAndStoreMedianPlots(std::string const &baseName, cmf::ChiSqrInfo &chiSqrInfo)
void Make2DHiddenParHistogram(art::TFileDirectory &tfd, std::string const &baseName, cmf::PointMap const &parVals)
static CovarianceBinUtility * Instance()
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