FitAlg_ROOTFactoryFit.cxx
Go to the documentation of this file.
4 
6 
7 //ROOT includes
8 #include "TMath.h"
9 
10 namespace fnex{
11 
12  // Must be outside of the class for ROOT factory functor ctor to be happy
14  std::unique_ptr<std::vector<double>> gFFVals;
15  std::unique_ptr<std::vector<std::string>> gFFNames;
16 
18  ROOT::Math::Minimizer * gFFMin;
19 
20  //----------------------------------------------------------------------------
21  double FitFunction_ForROOT(const double * params)
22  {
23  // debug
24  LOG_VERBATIM("ROOTFactoryFit") << "Processing guess: ";
25  for(size_t i = 0; i < gFFVals->size(); ++i){
26  if((*gFFNames)[i].compare("Th23") == 0)
27  LOG_VERBATIM("ROOTFactoryFit")
28  << "\t"
29  << (*gFFNames)[i]
30  << " : "
31  << params[i]
32  << "; sin^2(Th23) : "
33  << std::pow(std::sin(params[i]), 2.);
34  else if((*gFFNames)[i].compare("dCP") == 0)
35  LOG_VERBATIM("ROOTFactoryFit")
36  << "\t"
37  << (*gFFNames)[i]
38  << " : "
39  << params[i]
40  << " = "
41  << params[i]/TMath::Pi()
42  << " (pi)";
43  else
44  LOG_VERBATIM("ROOTFactoryFit")
45  << "\t"
46  << (*gFFNames)[i]
47  << " : "
48  << params[i];
49 
50  }
51 
52  std::vector<double> guessVals = *gFFVals;
53  for(size_t i = 0; i < gFFVals->size(); ++i) guessVals[i] = params[i];
54 
55  // Now make a new InputPoint
56  InputPoint newPoint = *gPoint;
57  ConvertVectorsToPoint(guessVals, *gFFNames, newPoint);
58 
59  // Send this to the MultiExperiment, and update the experiment
60  gRFFMultiExperiment->ApplyPoint(newPoint);
61 
62  //LOG_VERBATIM("ROOTFactoryFit") << "\t" << "Applied Point = " << newPoint << std::endl;
63 
64  // Get Chisq
65  double chisq = gRFFMultiExperiment->FitValue();
66 
67  LOG_VERBATIM("ROOTFactoryFit") << "\t" << "FitValue = " << chisq << std::endl;
68 
69  // Return the fit value; that's what we're looking to minimize
70  return chisq;
71 
72  }
73 
74  //----------------------------------------------------------------------------
76  : FitAlgorithm(pset)
77  {
78  this->Reconfigure(pset);
79  gFFVals = std::make_unique<std::vector<double> >(0);
80  gFFNames = std::make_unique<std::vector<std::string> >(0);
81  }
82 
83  // Begin normal class functions
84 
85  //----------------------------------------------------------------------------
87 
88  // Look for minName and algoName (following ROOT factory nomenclature,
89  // see NumericalMinimization.C in ROOT fit documentation,
90  // root.cern.ch/root/html/tutorials/fit/
91  // NumericalMinimization.C.html
92  // Default here is set to Minuit2 && Migrad
93  fMinName = pset.get< std::string >("MinName", "Minuit2");
94  fAlgoName = pset.get< std::string >("AlgoName", "Migrad");
95  fMaxFunctionCalls = pset.get< unsigned int >("MaxFunctionCalls", 1000);
96  fMaxIterations = pset.get< unsigned int >("MaxIterations", 1000);
97  fTolerance = pset.get< double >("Tolerance", 0);
98  fPrecision = pset.get< double >("Precision", 0);
99  fStrategyCode = pset.get< int >("StrategyCode", 0);
100  fPrintLevel = pset.get< int >("MINUITPrintLevel", 3);
101 
102  // The MultiExperiment is configured by the baseclass ctor.
103 
104  // We're done here
105  fIsConfigured = true;
106  }
107 
108  //----------------------------------------------------------------------------
109  bool FitAlg_ROOTFactoryFit::Fit(EventListMap const& fEventLists,
110  InputPoint const& pInputPoint,
111  PointResult & pPointResult){
112 
113 
114  // Run through the fit sequence
115  //std::vector<double> guessVals = *gFFVals;
116  //for (int i = 0;i < gFFVals->size();i++)
117  //{
118  // LOG_VERBATIM("ROOTFactoryFit") << guessVals[i];
119  //}
120 
121  LOG_VERBATIM("ROOTFactoryFit") << "Fit::Going to fill global variables\n";
122  FillGlobalVariables(fEventLists, pInputPoint);
123 
124  LOG_VERBATIM("ROOTFactoryFit") << "Fit::Going to initialize minimizer\n";
126 
127  LOG_VERBATIM("ROOTFactoryFit") << "Fit::Going to set initial values\n";
129 
130  LOG_VERBATIM("ROOTFactoryFit") << "Fit::Going to set constraints\n";
131  SetConstraints();
132 
133  LOG_VERBATIM("ROOTFactoryFit") << "Fit::About to minimize\n";
134  bool tConverged = Minimize();
135 
136  LOG_VERBATIM("ROOTFactoryFit") << "Fit::Going to load final results\n";
137  LoadFinalResults(pPointResult);
138 
139  // Return true if the fit was successful; else false
140  return tConverged;
141 
142  }
143 
144 
145  //----------------------------------------------------------------------------
146  // -----
147  // Set up the global vectors / variables needed for the fit
148  // -----
150  InputPoint const& pInputPoint){
151 
152  // Send the event list to the multiexperiment
153  //LOG_VERBATIM("ROOTFactoryFit") << "WTF HERE 0\n";
154  fMultiExperiment.SetEventListMap(fEventLists);
155 
156  // Set the input multiexperiment to the global used by the FitValue function
157  //LOG_VERBATIM("ROOTFactoryFit") << "WTF HERE 1\n";
158  gRFFMultiExperiment = &fMultiExperiment;
159 
160  //LOG_VERBATIM("ROOTFactoryFit") << "WTF HERE 2\n";
161  gPoint = new InputPoint(pInputPoint);
162 
163  //LOG_VERBATIM("ROOTFactoryFit") << "WTF HERE 3\n";
164  ConvertPointToVectors(pInputPoint, *gFFVals, *gFFNames);
165 
166  // We're done here
167  return true;
168  }
169 
170  //----------------------------------------------------------------------------
171  // -----
172  // Set up the ROOT::Math::Factory minimizer (basic)
173  // -----
175 
176  // Make the minimizer object
177  gFFMin = ROOT::Math::Factory::CreateMinimizer(fMinName.data(), fAlgoName.data());
178 
179  // Set up the minimizer particulars
180  gFFMin->SetMaxFunctionCalls(fMaxFunctionCalls);
181  gFFMin->SetMaxIterations(fMaxIterations);
182  gFFMin->SetStrategy(fStrategyCode);
183  if(fTolerance > 0) gFFMin->SetTolerance(fTolerance);
184  if(fPrecision > 0) gFFMin->SetPrecision(fPrecision);
185  gFFMin->SetPrintLevel(fPrintLevel);
186 
187  // Vector sizes for bookkeeping
188  size_t numFloatables = gFFNames->size();
189  LOG_VERBATIM("ROOTFactoryFit")
190  << "\nNum floatables: " << numFloatables
191  << "\nMax function calls: " << fMaxFunctionCalls
192  << "\nMax iterations: " << fMaxIterations
193  << "\nTolerance: " << fTolerance;
194  LOG_VERBATIM("ROOTFactoryFit")
195  << "Initial guess: ";
196  for(size_t i = 0; i < gFFNames->size(); ++i){
197  LOG_VERBATIM("ROOTFactoryFit")
198  << "\n" << (*gFFNames)[i]
199  << " : " << (*gFFVals)[i];
200  }
201 
202  // Link in the function to be minimized
203  ROOT::Math::Functor * f = new ROOT::Math::Functor(&FitFunction_ForROOT,numFloatables);
204  gFFMin->SetFunction(*f);
205 
206  // We're done here
207  return true;
208  }
209 
210  //----------------------------------------------------------------------------
212 
213  // Set initial values for all parameterSpaceLocations
214  for(size_t ipar = 0; ipar < gFFNames->size(); ++ipar){
215  gFFMin->SetVariable(ipar,
216  (*gFFNames)[ipar],
217  (*gFFVals)[ipar],
218  0.01);
219  }
220 
221  // We're done here
222  return true;
223 
224  }
225 
226 
227  //----------------------------------------------------------------------------
228  // -----
229  // Set up the ROOT::Math::Factory minimizer variable constraints
230  // -----
232  {
233 
234  // Set constraints for all parameter space locations
237 
238  std::string key("");
239 
240  fnex::ConstraintInfo constraint;
241 
242  for(auto const & itr : constraints){
243 
244  key = itr.first;
245  constraint = itr.second;
246 
247  int ikey = -1;
248  for(size_t tkey = 0; tkey < gFFNames->size(); ++tkey){
249  if( (*gFFNames)[tkey] == key) ikey = tkey;
250  }
251  if(ikey == -1) continue;
252 
253  LOG_DEBUG("ROOTFactoryFit")
254  << " Key: " << key
255  << " IKey: " << ikey;
256 
257  if(constraint.LoValue() == constraint.HiValue()){
258  gFFMin->SetVariable(ikey,
259  (*gFFNames)[ikey],
260  constraint.LoValue(),
261  0.01);
262  gFFMin->FixVariable(ikey);
263  }else{
264  LOG_VERBATIM("ROOTFactoryFit")
265  << "Setting constraint on "
266  << std::setw(10) << key << " "
267  << std::setw(10) << constraint.LoValue()
268  << std::setw(10) << constraint.HiValue() << std::endl;
269  gFFMin->SetVariableLimits(ikey,
270  constraint.LoValue(),
271  constraint.HiValue());
272  }
273  }
274 
275  return true;
276  }
277 
278 
279 
280  //----------------------------------------------------------------------------
281  // -----
282  // Minimize!
283  // -----
285  //LOG_VERBATIM("ROOTFactoryFit") << "In Minimize. Calling ROOTFactorFit->Minimize()\n";
286  return gFFMin->Minimize();
287  }
288 
289  //----------------------------------------------------------------------------
290  // -----
291  // Fill the PointResult object with the converged-upon values
292  // -----
294 
295  // Fill the PointResult object to contain the converged-upon point
296  const double * vals = gFFMin->X();
297  std::vector<double> fitVals = *gFFVals;
298  for(size_t i = 0; i < gFFVals->size(); ++i){
299  fitVals[i] = vals[i];
300  }
301  //const double * errs = gFFMin->Errors();
302 
303  //Add a 2.*TMath::Pi() if the dcp < 0.
304  std::vector<double> guessVals = *gFFVals;
305  for(size_t i = 0; i < gFFVals->size(); ++i){
306  guessVals[i] = fitVals[i];
307 
308  if((*gFFNames)[i].compare("dCP") == 0 && guessVals[i] < 0)
309  guessVals[i] = guessVals[i] + 2.*TMath::Pi();
310  }
311 
312  // Now make a new InputPoint
313  InputPoint newPoint = *gPoint;
314  ConvertVectorsToPoint(guessVals, *gFFNames, newPoint);
315  gRFFMultiExperiment->ApplyPoint(newPoint);
316 
317  int status = gFFMin->Status();
318  double chisq = gFFMin->MinValue();
319 
320  //InputPoint newPoint = *gPoint;
321  //ConvertVectorsToPoint(fitVals, *gFFNames, newPoint);
322  pPointResult = PointResult(chisq, status, newPoint.ParameterSpaceLocation(), newPoint.SystematicPulls());
323 
324  return true;
325  }
326 
327 
328  //----------------------------------------------------------------------------
329  // Convert between InputPoint and std::vector formats
330  void ConvertPointToVectors (InputPoint const & pPoint,
331  std::vector<double> & pVals,
332  std::vector<std::string> & pNames){
333 
334  // Grab the easy-to-get vectors
335  fnex::ParameterDetMap pParameterSpaceLocations = pPoint.ParameterSpaceLocation();
336  fnex::ParameterDetMap pSystematicPulls = pPoint.SystematicPulls();
337  std::vector<fnex::ParameterDet> pFixedParameters = pPoint.FixedParameters();
338  std::vector<fnex::ParameterDet> pFixedSystematics = pPoint.FixedSystematics();
339 
340 
341  // Add all FD entry keys to the psl_key_list
342  // Add ND entry keys only if they don't shot up on FD list
343  // E.g., FD entry will override ND entry if common
344  // Don't add entries on the Fixed list
345  std::list<std::string> psl_key_list_fd;
346  for( auto psl : pPoint.ParameterSpaceLocation(novadaq::cnv::kFARDET) ){
347  ParameterDet tempPairFar = std::make_pair(psl.first, novadaq::cnv::kFARDET);
348  if(count(pFixedParameters.begin(), pFixedParameters.end(), tempPairFar) != 0) continue;
349  //LOG_VERBATIM("ROOTFactoryFit") << "FD Adding " << psl.first << std::endl;
350  psl_key_list_fd.push_back(psl.first);
351  }
352 
353  // The ND list will contain an event only if it occurs solely for the ND
354  // and not for the FD
355  std::list<std::string> psl_key_list_nd;
356  for( auto psl : pPoint.ParameterSpaceLocation(novadaq::cnv::kNEARDET) ){
357  ParameterDet tempPairNear = std::make_pair(psl.first, novadaq::cnv::kNEARDET);
358  if(count(pFixedParameters.begin(), pFixedParameters.end(), tempPairNear) != 0) continue;
359  if( count(psl_key_list_fd.begin(), psl_key_list_fd.end(), psl.first) == 0 ){
360  //LOG_VERBATIM("ROOTFactoryFit") << "ND Adding " << psl.first << std::endl;
361  psl_key_list_nd.push_back(psl.first);
362  }
363  }
364  psl_key_list_nd.sort();
365  psl_key_list_fd.sort();
366 
367  // Add all FD entry keys to the sp_key_list
368  // Add ND entry keys only if they don't shot up on FD list
369  // E.g., FD entry will override ND entry if common
370  // Don't add entries on the Fixed list
371  std::list<std::string> sp_key_list_fd;
372  for( auto sp : pPoint.SystematicPulls(novadaq::cnv::kFARDET) ){
373  ParameterDet tempPairFar = std::make_pair(sp.first, novadaq::cnv::kFARDET);
374  if(count(pFixedSystematics.begin(), pFixedSystematics.end(), tempPairFar) != 0) continue;
375  sp_key_list_fd.push_back(sp.first);
376  }
377 
378  // Again, the ND list will contain an event only if it occurs solely for the ND
379  // and not for the FD
380  std::list<std::string> sp_key_list_nd;
381  for( auto sp : pPoint.SystematicPulls(novadaq::cnv::kNEARDET) ){
382  ParameterDet tempPairNear = std::make_pair(sp.first, novadaq::cnv::kNEARDET);
383  if(count(pFixedSystematics.begin(), pFixedSystematics.end(), tempPairNear) != 0) continue;
384  if( count(sp_key_list_fd.begin(), sp_key_list_fd.end(), sp.first) == 0 ){
385  sp_key_list_nd.push_back(sp.first);
386  }
387  }
388 
389  sp_key_list_nd.sort();
390  sp_key_list_fd.sort();
391 
392  // Now fill the pVals and pNames vectors with these 'condensed' vectors
393  //
394  pVals.clear();
395  pNames.clear();
396 
397 
398  for (auto it=psl_key_list_nd.begin(); it != psl_key_list_nd.end(); ++it) {
399  pVals.push_back(pPoint.ParameterSpaceLocation(novadaq::cnv::kNEARDET).at(*it));
400  pNames.push_back(*it);
401  }
402  for (auto it=psl_key_list_fd.begin(); it != psl_key_list_fd.end(); ++it) {
403  pVals.push_back(pPoint.ParameterSpaceLocation(novadaq::cnv::kFARDET).at(*it));
404  pNames.push_back(*it);
405  }
406 
407  for (auto it=sp_key_list_nd.begin(); it != sp_key_list_nd.end(); ++it) {
408  pVals.push_back(pPoint.SystematicPulls(novadaq::cnv::kNEARDET).at(*it));
409  pNames.push_back(*it);
410  }
411  for (auto it=sp_key_list_fd.begin(); it != sp_key_list_fd.end(); ++it) {
412  pVals.push_back(pPoint.SystematicPulls(novadaq::cnv::kFARDET).at(*it));
413  pNames.push_back(*it);
414  }
415 
416  // We're done here
417  return;
418  }
419 
420 
421  //----------------------------------------------------------------------------
422  void ConvertVectorsToPoint(std::vector<double> const& pVals,
423  std::vector<std::string> const& pNames,
424  InputPoint & pPoint){
425 
426  // Grab the vectors to overwrite in the InputPoint
427  //
428  ParameterDetMap pParameterSpaceLocations = pPoint.ParameterSpaceLocation();
429  ParameterDetMap pSystematicPulls = pPoint.SystematicPulls();
430  std::vector<ParameterDet> pFixedParameters = pPoint.FixedParameters();
431  std::vector<ParameterDet> pFixedSystematics = pPoint.FixedSystematics();
432 
433 
434  // Start by looping through condensed name (key) vector
435  //
436  unsigned int keys_size = pNames.size();
437  for (unsigned int i = 0; i < keys_size; ++i) {
438  ParameterDet tempPairFar = std::make_pair(pNames.at(i), novadaq::cnv::kFARDET);
439  ParameterDet tempPairNear = std::make_pair(pNames.at(i), novadaq::cnv::kNEARDET);
440  if(pParameterSpaceLocations.count(tempPairFar) != 0){
441  pParameterSpaceLocations[tempPairFar] = pVals.at(i); }
442  if(pParameterSpaceLocations.count(tempPairNear) != 0){
443  pParameterSpaceLocations[tempPairNear] = pVals.at(i); }
444  if(pSystematicPulls.count(tempPairFar) != 0){
445  pSystematicPulls[tempPairFar] = pVals.at(i); }
446  if(pSystematicPulls.count(tempPairNear) != 0){
447  pSystematicPulls[tempPairNear] = pVals.at(i); }
448  }
449 
450  // Make the new input point
451  //
452  pPoint = InputPoint(pParameterSpaceLocations,
453  pSystematicPulls,
454  pFixedParameters,
455  pFixedSystematics );
456 
457  // We're done here
458  return;
459  }
460 
461 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void Print(std::ostream &stream=std::cout)
double const LoValue() const
int status
Definition: fabricate.py:1613
set< int >::iterator it
int fStrategyCode
Use 0 for faster fits, 1 for better errors.
constexpr T pow(T x)
Definition: pow.h:75
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
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
Create a list of fnex::Events to be used in fits.
bool LoadFinalResults(PointResult &pPointResult)
constraints
Definition: mkDefs.py:147
FitAlg_ROOTFactoryFit(fhicl::ParameterSet const &pset)
std::vector< std::string > const & FixedParameters(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:396
void Reconfigure(fhicl::ParameterSet const &pset)
fnex::MultiExperiment * gRFFMultiExperiment
InputPoint * gPoint
bool SetEventListMap(EventListMap const &pEventListMap)
std::unique_ptr< std::vector< std::string > > gFFNames
ROOT::Math::Minimizer * gFFMin
MultiExperiment fMultiExperiment
Definition: FitAlgorithm.h:50
std::vector< std::string > const & FixedSystematics(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:404
void ConvertVectorsToPoint(std::vector< double > const &pVals, std::vector< std::string > const &pNames, InputPoint &pPoint)
int fPrintLevel
3 is default, more output would be 4
Far Detector at Ash River, MN.
double FitFunction_ForROOT(const double *params)
std::map< std::string, ConstraintInfo > const & Constraints()
static ParameterConstraints * Instance()
ParameterMap const ParameterSpaceLocation(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:330
T get(std::string const &key) const
Definition: ParameterSet.h:231
Near Detector in the NuMI cavern.
void ApplyPoint(InputPoint &pInputPoint)
std::unordered_map< ParameterDet, double, ParameterDetHasher > ParameterDetMap
Definition: Constants.h:63
std::unique_ptr< std::vector< double > > gFFVals
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
bool compare(const GFluxGenerator &g1, const GFluxGenerator &g2)
double const HiValue() const
std::pair< std::string, novadaq::cnv::DetId > ParameterDet
Definition: Constants.h:38
bool Fit(EventListMap const &fEventLists, InputPoint const &pInputPoint, PointResult &pPointResult)
#define LOG_VERBATIM(category)
bool FillGlobalVariables(EventListMap const &fEventLists, InputPoint const &pInputPoint)
void ConvertPointToVectors(InputPoint const &pPoint, std::vector< double > &pVals, std::vector< std::string > &pNames)
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:364