laura is hosted by Hepforge, IPPP Durham
Laura++  v3r3
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauMinuit.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2013 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 #include <algorithm>
17 
18 #include "TMatrixD.h"
19 #include "TVirtualFitter.h"
20 
21 #include "LauFitObject.hh"
22 #include "LauFitter.hh"
23 #include "LauMinuit.hh"
24 #include "LauParameter.hh"
25 #include "LauParamFixed.hh"
26 
27 // It's necessary to define an external function that specifies the address of the function
28 // that Minuit needs to minimise. Minuit doesn't know about any classes - therefore
29 // use gMinuit->SetFCN(external_function), gMinuit->SetObjectFit(this).
30 // Here, we use TVirtualFitter* fitter instead of gMinuit, defined below.
31 // Then, within the external function, invoke an object from this class (LauAllModel),
32 // and use the member functions to access the parameters/variables.
33 extern void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
34 
35 
36 LauMinuit::LauMinuit( Int_t maxPar ) : LauAbsFitter(),
37  minuit_(0),
38  maxPar_(maxPar),
39  nParams_(0),
40  nFreeParams_(0),
41  twoStageFit_(kFALSE),
42  useAsymmFitErrors_(kFALSE),
43  fitStatus_({-1,0.0,0.0})
44 {
45  TVirtualFitter::SetDefaultFitter( "Minuit" );
46  minuit_ = TVirtualFitter::Fitter( 0, maxPar_ );
47 }
48 
50 {
51 }
52 
53 void LauMinuit::initialise( LauFitObject* fitObj, const std::vector<LauParameter*>& parameters )
54 {
55  // Check whether we're going to use asymmetric errors
56  if (useAsymmFitErrors_ == kTRUE) {
57  std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl;
58  std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl;
59  } else {
60  std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl;
61  }
62 
63  // Store the parameters
64  params_ = parameters;
65 
66  // Hook the external likelihood function to this LauFitter::fitter() and this class.
67  minuit_->SetFCN( logLikeFun );
68  minuit_->SetObjectFit( fitObj );
69 
70  // Clear any stored parameters etc... before using
71  minuit_->Clear();
72 
73  nParams_ = params_.size();
74  std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl;
75  std::cout << " : Total number of parameters = " << nParams_ << std::endl;
76 
77  // Define the default relative error
78  const Double_t defaultError(0.01);
79 
80  // Set-up the parameters
81  for (UInt_t i = 0; i < nParams_; ++i) {
82  TString name = params_[i]->name();
83  Double_t initVal = params_[i]->initValue();
84  Double_t initErr = params_[i]->error();
85  // If we do not have a supplied estimate of the error, we should make a reasonable guess
86  if ( initErr == 0.0 ) {
87  if ( initVal == 0.0 ) {
88  initErr = defaultError;
89  } else if ( TMath::Abs(initErr/initVal) < 1e-6 ) {
90  initErr = TMath::Abs(defaultError * initVal);
91  }
92  }
93  Double_t minVal = params_[i]->minValue();
94  Double_t maxVal = params_[i]->maxValue();
95  Bool_t secondStage = params_[i]->secondStage();
96  if (this->twoStageFit() && secondStage == kTRUE) {
97  params_[i]->fixed(kTRUE);
98  }
99  Bool_t fixVar = params_[i]->fixed();
100 
101  std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl;
102  minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal);
103 
104  // Fix parameter if required
105  if (fixVar == kTRUE) {
106  std::cout << " : Fixing parameter " << i << std::endl;
107  minuit_->FixParameter(i);
108  }
109  }
110 
111  LauParamFixed pred;
112  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
113 
114  // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
115  // for maximum likelihood fit. Very important command, otherwise all
116  // extracted errors will be too big, and pull distributions will be too narrow!
117  // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L)
118  Double_t argL[2];
119  argL[0] = 0.5;
120  fitStatus_.status = minuit_->ExecuteCommand("SET ERR", argL, 1);
121 
122  //argL[0] = 0;
123  //fitStatus_.status = minuit_->ExecuteCommand("SET STRATEGY", argL, 1);
124 }
125 
127 {
128  return (minuit_!=0) ? dynamic_cast<LauFitObject*>( minuit_->GetObjectFit() ) : 0;
129 }
130 
132 {
133  Double_t arglist[2];
134  arglist[0] = 1000*nParams_; // maximum iterations
135  arglist[1] = 0.05; // tolerance -> min EDM = 0.001*tolerance (0.05)
136  fitStatus_.status = minuit_->ExecuteCommand("MIGRAD", arglist, 2);
137 
138  // Dummy variables - need to feed them to the function
139  // used for getting NLL, EDM and error matrix status
140  Double_t errdef;
141  Int_t nvpar, nparx;
142 
143  if (fitStatus_.status != 0) {
144 
145  std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl;
146 
147  } else {
148 
149  // Check that the error matrix is ok
150  fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx);
151  std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_.status << std::endl;
152  // 0= not calculated at all
153  // 1= approximation only, not accurate
154  // 2= full matrix, but forced positive-definite
155  // 3= full accurate covariance matrix
156 
157  // Fit result was OK. Now get the more precise errors.
158  fitStatus_.status = minuit_->ExecuteCommand("HESSE", arglist, 1);
159 
160  if (fitStatus_.status != 0) {
161 
162  std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl;
163 
164  } else {
165 
166  // Check that the error matrix is ok
167  fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx);
168  std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_.status << std::endl;
169  // 0= not calculated at all
170  // 1= approximation only, not accurate
171  // 2= full matrix, but forced positive-definite
172  // 3= full accurate covariance matrix
173 
174  // Symmetric errors and eror matrix were OK.
175  // Get asymmetric errors if asked for.
176  if (useAsymmFitErrors_ == kTRUE) {
177  LauFitObject* fitObj = this->getFitObject();
178  fitObj->withinAsymErrorCalc( kTRUE );
179  fitStatus_.status = minuit_->ExecuteCommand("MINOS", arglist, 1);
180  fitObj->withinAsymErrorCalc( kFALSE );
181  if (fitStatus_.status != 0) {
182  std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl;
183  }
184  }
185  }
186  }
187 
188  // Print results
189  fitStatus_.status = minuit_->GetStats(fitStatus_.NLL, fitStatus_.EDM, errdef, nvpar, nparx);
190  std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_.status << std::endl;
191  // 0= not calculated at all
192  // 1= approximation only, not accurate
193  // 2= full matrix, but forced positive-definite
194  // 3= full accurate covariance matrix
195 
196  minuit_->PrintResults(3, fitStatus_.NLL);
197 
198  // Retrieve the covariance matrix from the fitter
199  // For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_
200  // but only the elements within the sub-matrix nFreeParams_ x nFreeParams_ have values,
201  // the "trailing" elements are zero, so we trim them off.
202  Double_t* covMatrix = minuit_->GetCovarianceMatrix();
203  covMatrix_.Clear();
204  covMatrix_.ResizeTo( nParams_, nParams_ );
205  covMatrix_.SetMatrixArray( covMatrix );
206  covMatrix_.ResizeTo( nFreeParams_, nFreeParams_ );
207 
208  return fitStatus_;
209 }
210 
212 {
213  for (UInt_t i = 0; i < nParams_; ++i) {
214  Bool_t secondStage = params_[i]->secondStage();
215  if (secondStage == kTRUE) {
216  params_[i]->fixed(kTRUE);
217  minuit_->FixParameter(i);
218  }
219  }
220  LauParamFixed pred;
221  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
222 }
223 
225 {
226  for (UInt_t i = 0; i < nParams_; ++i) {
227  Bool_t secondStage = params_[i]->secondStage();
228  if (secondStage == kTRUE) {
229  params_[i]->fixed(kFALSE);
230  minuit_->ReleaseParameter(i);
231  }
232  }
233  LauParamFixed pred;
234  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
235 }
236 
238 {
239  for (UInt_t i = 0; i < nParams_; ++i) {
240  // Get the value and errors from MINUIT
241  Double_t value = minuit_->GetParameter(i);
242  Double_t error(0.0);
243  Double_t negError(0.0);
244  Double_t posError(0.0);
245  Double_t globalcc(0.0);
246  minuit_->GetErrors(i, posError, negError, error, globalcc);
247  params_[i]->valueAndErrors(value, error, negError, posError);
248  params_[i]->globalCorrelationCoeff(globalcc);
249  }
250 }
251 
252 // Definition of the fitting function for Minuit
253 void logLikeFun(Int_t& npar, Double_t* /*first_derivatives*/, Double_t& f, Double_t* par, Int_t /*iflag*/)
254 {
255  // Routine that specifies the negative log-likelihood function for the fit.
256  // Used by the MINUIT minimising code.
257 
259 
260  // Set the internal parameters for this model using parameters from Minuit (pars):
261  theModel->setParsFromMinuit( par, npar );
262 
263  // Set the value of f to be the total negative log-likelihood for the data sample.
264  f = theModel->getTotNegLogLikelihood();
265 }
266 
UInt_t nParams_
The total number of parameters.
Definition: LauMinuit.hh:125
Double_t EDM
The estimated distance to the minimum.
Definition: LauAbsFitter.hh:47
virtual void releaseSecondStageParameters()
Release parameters marked as &quot;second stage&quot;.
Definition: LauMinuit.cc:224
virtual Double_t getTotNegLogLikelihood()=0
Calculate the new value of the negative log likelihood.
File containing declaration of LauParamFixed class.
FitStatus fitStatus_
The status of the fit.
Definition: LauMinuit.hh:137
virtual void updateParameters()
Update the values and errors of the parameters based on the fit minimum.
Definition: LauMinuit.cc:237
const TString & name() const
The parameter name.
virtual void initialise(LauFitObject *fitObj, const std::vector< LauParameter * > &parameters)
Initialise the fitter, setting the information on the parameters.
Definition: LauMinuit.cc:53
virtual ~LauMinuit()
Destructor.
Definition: LauMinuit.cc:49
LauMinuit(Int_t maxPar=100)
Constructor.
Definition: LauMinuit.cc:36
virtual void setParsFromMinuit(Double_t *par, Int_t npar)=0
This function sets the parameter values from Minuit.
Double_t negError() const
The lower error on the parameter.
void logLikeFun(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: LauMinuit.cc:253
Double_t posError() const
The upper error on the parameter.
File containing declaration of LauFitter class.
The abstract interface to the fitter.
Definition: LauAbsFitter.hh:34
minuit_
Definition: LauMinuit.cc:46
Double_t NLL
The negative log-likelihood.
Definition: LauAbsFitter.hh:45
virtual Bool_t twoStageFit() const
Determine whether the two-stage fit is enabled.
Definition: LauMinuit.hh:59
Int_t status
The status code of the fit.
Definition: LauAbsFitter.hh:43
TMatrixD covMatrix_
The covariance matrix.
Definition: LauMinuit.hh:140
File containing declaration of LauParameter class.
Struct to store fit status information.
Definition: LauAbsFitter.hh:41
Predicate to allow counting of the number of fixed parameters.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
Double_t error() const
The error on the parameter.
Bool_t useAsymmFitErrors_
Option to use asymmetric errors.
Definition: LauMinuit.hh:134
virtual void withinAsymErrorCalc(const Bool_t inAsymErrCalc)
Mark that the fit is calculating asymmetric errors.
Definition: LauFitObject.hh:68
static LauAbsFitter * fitter()
Method that provides access to the singleton fitter.
Definition: LauFitter.cc:34
virtual LauFitObject * getFitObject()=0
Get the object that controls the calculation of the likelihood.
UInt_t nFreeParams_
The number of free parameters.
Definition: LauMinuit.hh:128
File containing declaration of LauFitObject class.
virtual const FitStatus & minimise()
Perform the minimisation of the fit function.
Definition: LauMinuit.cc:131
TVirtualFitter * minuit_
The interface to Minuit.
Definition: LauMinuit.hh:116
Double_t value() const
The value of the parameter.
virtual void fixSecondStageParameters()
Fix parameters marked as &quot;second stage&quot;.
Definition: LauMinuit.cc:211
virtual LauFitObject * getFitObject()
Get the object that controls the calculation of the likelihood.
Definition: LauMinuit.cc:126
std::vector< LauParameter * > params_
The fit parameters.
Definition: LauMinuit.hh:122
File containing declaration of LauMinuit class.
The abstract interface for the objects that control the calculation of the likelihood.
Definition: LauFitObject.hh:30