laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2p1
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_(0),
44  NLL_(0.0)
45 {
46  TVirtualFitter::SetDefaultFitter( "Minuit" );
47  minuit_ = TVirtualFitter::Fitter( 0, maxPar_ );
48 }
49 
51 {
52 }
53 
54 void LauMinuit::initialise( LauFitObject* fitObj, const std::vector<LauParameter*>& parameters )
55 {
56  // Check whether we're going to use asymmetric errors
57  if (useAsymmFitErrors_ == kTRUE) {
58  std::cout << "INFO in LauMinuit::fit : We are going to calculate the asymmetric fit errors." << std::endl;
59  std::cout << " : This will, in general, significantly increase the CPU time required for fitting." << std::endl;
60  } else {
61  std::cout << "INFO in LauMinuit::fit : Only parabolic errors will be calculated." << std::endl;
62  }
63 
64  // Store the parameters
65  params_ = parameters;
66 
67  // Hook the external likelihood function to this LauFitter::fitter() and this class.
68  minuit_->SetFCN( logLikeFun );
69  minuit_->SetObjectFit( fitObj );
70 
71  // Clear any stored parameters etc... before using
72  minuit_->Clear();
73 
74  nParams_ = params_.size();
75  std::cout << "INFO in LauMinuit::initialise : Setting fit parameters" << std::endl;
76  std::cout << " : Total number of parameters = " << nParams_ << std::endl;
77 
78  // Define the default relative error
79  const Double_t defaultError(0.01);
80 
81  // Set-up the parameters
82  for (UInt_t i = 0; i < nParams_; ++i) {
83  TString name = params_[i]->name();
84  Double_t initVal = params_[i]->initValue();
85  Double_t initErr = params_[i]->error();
86  if ( initVal == 0.0 ) {
87  initErr = defaultError;
88  } else if ( TMath::Abs(initErr/initVal) < 1e-6 ) {
89  initErr = TMath::Abs(defaultError * initVal);
90  }
91  Double_t minVal = params_[i]->minValue();
92  Double_t maxVal = params_[i]->maxValue();
93  Bool_t secondStage = params_[i]->secondStage();
94  if (this->twoStageFit() && secondStage == kTRUE) {
95  params_[i]->fixed(kTRUE);
96  }
97  Bool_t fixVar = params_[i]->fixed();
98 
99  std::cout << " : Setting parameter " << i << " called " << name << " to have initial value " << initVal << ", error " << initErr << " and range " << minVal << " to " << maxVal << std::endl;
100  minuit_->SetParameter(i, name, initVal, initErr, minVal, maxVal);
101 
102  // Fix parameter if required
103  if (fixVar == kTRUE) {
104  std::cout << " : Fixing parameter " << i << std::endl;
105  minuit_->FixParameter(i);
106  }
107  }
108 
109  LauParamFixed pred;
110  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
111 
112  // Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
113  // for maximum likelihood fit. Very important command, otherwise all
114  // extracted errors will be too big, and pull distributions will be too narrow!
115  // TODO - The alternative to this is to make FCN = -2log(L) rather than -log(L)
116  Double_t argL[2];
117  argL[0] = 0.5;
118  fitStatus_ = minuit_->ExecuteCommand("SET ERR", argL, 1);
119 
120  //argL[0] = 0;
121  //fitStatus_ = minuit_->ExecuteCommand("SET STRATEGY", argL, 1);
122 }
123 
125 {
126  return (minuit_!=0) ? dynamic_cast<LauFitObject*>( minuit_->GetObjectFit() ) : 0;
127 }
128 
129 std::pair<Int_t,Double_t> LauMinuit::minimise()
130 {
131  Double_t arglist[2];
132  arglist[0] = 1000*nParams_; // maximum iterations
133  arglist[1] = 0.05; // tolerance -> min EDM = 0.001*tolerance (0.05)
134  fitStatus_ = minuit_->ExecuteCommand("MIGRAD", arglist, 2);
135 
136  // Dummy variables - need to feed them to the function
137  // used for getting NLL and error matrix status
138  Double_t edm, errdef;
139  Int_t nvpar, nparx;
140 
141  if (fitStatus_ != 0) {
142 
143  std::cerr << "ERROR in LauMinuit::minimise : Error in minimising loglike." << std::endl;
144 
145  } else {
146 
147  // Check that the error matrix is ok
148  NLL_ = 0.0;
149  fitStatus_ = minuit_->GetStats(NLL_, edm, errdef, nvpar, nparx);
150  std::cout << "INFO in LauMinuit::minimise : Error matrix status after MIGRAD is: " << fitStatus_ << std::endl;
151  // 0= not calculated at all
152  // 1= approximation only, not accurate
153  // 2= full matrix, but forced positive-definite
154  // 3= full accurate covariance matrix
155 
156  // Fit result was OK. Now get the more precise errors.
157  fitStatus_ = minuit_->ExecuteCommand("HESSE", arglist, 1);
158 
159  if (fitStatus_ != 0) {
160 
161  std::cerr << "ERROR in LauMinuit::minimise : Error in HESSE routine." << std::endl;
162 
163  } else {
164 
165  // Check that the error matrix is ok
166  NLL_ = 0.0;
167  fitStatus_ = minuit_->GetStats(NLL_, edm, errdef, nvpar, nparx);
168  std::cout << "INFO in LauMinuit::minimise : Error matrix status after HESSE is: " << fitStatus_ << 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_ = minuit_->ExecuteCommand("MINOS", arglist, 1);
180  fitObj->withinAsymErrorCalc( kFALSE );
181  if (fitStatus_ != 0) {
182  std::cerr << "ERROR in LauMinuit::minimise : Error in MINOS routine." << std::endl;
183  }
184  }
185  }
186  }
187 
188  // Print results
189  NLL_ = 0.0;
190  fitStatus_ = minuit_->GetStats(NLL_, edm, errdef, nvpar, nparx);
191  std::cout << "INFO in LauMinuit::minimise : Final error matrix status is: " << fitStatus_ << std::endl;
192  // 0= not calculated at all
193  // 1= approximation only, not accurate
194  // 2= full matrix, but forced positive-definite
195  // 3= full accurate covariance matrix
196 
197  minuit_->PrintResults(3, NLL_);
198 
199  // Retrieve the covariance matrix from the fitter
200  // For some reason the array returned is as if the matrix is of dimension nParams_ x nParams_
201  // but only the elements within the sub-matrix nFreeParams_ x nFreeParams_ have values,
202  // the "trailing" elements are zero, so we trim them off.
203  Double_t* covMatrix = minuit_->GetCovarianceMatrix();
204  covMatrix_.Clear();
205  covMatrix_.ResizeTo( nParams_, nParams_ );
206  covMatrix_.SetMatrixArray( covMatrix );
207  covMatrix_.ResizeTo( nFreeParams_, nFreeParams_ );
208 
209  return std::make_pair( fitStatus_, NLL_ );
210 }
211 
213 {
214  for (UInt_t i = 0; i < nParams_; ++i) {
215  Bool_t firstStage = params_[i]->firstStage();
216  if (firstStage == kTRUE) {
217  params_[i]->fixed(kTRUE);
218  minuit_->FixParameter(i);
219  }
220  }
221  LauParamFixed pred;
222  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
223 }
224 
226 {
227  for (UInt_t i = 0; i < nParams_; ++i) {
228  Bool_t firstStage = params_[i]->firstStage();
229  if (firstStage == kTRUE) {
230  params_[i]->fixed(kFALSE);
231  minuit_->ReleaseParameter(i);
232  }
233  }
234  LauParamFixed pred;
235  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
236 }
237 
239 {
240  for (UInt_t i = 0; i < nParams_; ++i) {
241  Bool_t secondStage = params_[i]->secondStage();
242  if (secondStage == kTRUE) {
243  params_[i]->fixed(kTRUE);
244  minuit_->FixParameter(i);
245  }
246  }
247  LauParamFixed pred;
248  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
249 }
250 
252 {
253  for (UInt_t i = 0; i < nParams_; ++i) {
254  Bool_t secondStage = params_[i]->secondStage();
255  if (secondStage == kTRUE) {
256  params_[i]->fixed(kFALSE);
257  minuit_->ReleaseParameter(i);
258  }
259  }
260  LauParamFixed pred;
261  nFreeParams_ = nParams_ - std::count_if(params_.begin(),params_.end(),pred);
262 }
263 
265 {
266  for (UInt_t i = 0; i < nParams_; ++i) {
267  // Get the value and errors from MINUIT
268  Double_t value = minuit_->GetParameter(i);
269  Double_t error(0.0);
270  Double_t negError(0.0);
271  Double_t posError(0.0);
272  Double_t globalcc(0.0);
273  minuit_->GetErrors(i, posError, negError, error, globalcc);
274  params_[i]->valueAndErrors(value, error, negError, posError);
275  params_[i]->globalCorrelationCoeff(globalcc);
276  }
277 }
278 
279 // Definition of the fitting function for Minuit
280 void logLikeFun(Int_t& npar, Double_t* /*first_derivatives*/, Double_t& f, Double_t* par, Int_t /*iflag*/)
281 {
282  // Routine that specifies the negative log-likelihood function for the fit.
283  // Used by the MINUIT minimising code.
284 
286 
287  // Set the internal parameters for this model using parameters from Minuit (pars):
288  theModel->setParsFromMinuit( par, npar );
289 
290  // Set the value of f to be the total negative log-likelihood for the data sample.
291  f = theModel->getTotNegLogLikelihood();
292 }
293 
UInt_t nParams_
The total number of parameters.
Definition: LauMinuit.hh:131
Double_t NLL_
The negative log-likelihood.
Definition: LauMinuit.hh:146
const UInt_t maxPar_
The maximum number of parameters.
Definition: LauMinuit.hh:125
virtual void releaseSecondStageParameters()
Release parameters marked as &quot;second stage&quot;.
Definition: LauMinuit.cc:251
virtual std::pair< Int_t, Double_t > minimise()
Perform the minimisation of the fit function.
Definition: LauMinuit.cc:129
virtual Double_t getTotNegLogLikelihood()=0
Calculate the new value of the negative log likelihood.
File containing declaration of LauParamFixed class.
virtual void updateParameters()
Update the values and errors of the parameters based on the fit minimum.
Definition: LauMinuit.cc:264
const TString & name() const
The parameter name.
Int_t fitStatus_
The status of the fit.
Definition: LauMinuit.hh:143
virtual void initialise(LauFitObject *fitObj, const std::vector< LauParameter * > &parameters)
Initialise the fitter, setting the information on the parameters.
Definition: LauMinuit.cc:54
virtual ~LauMinuit()
Destructor.
Definition: LauMinuit.cc:50
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:280
Double_t posError() const
The upper error on the parameter.
File containing LauFitter class.
virtual void releaseFirstStageParameters()
Release parameters marked as &quot;first stage&quot;.
Definition: LauMinuit.cc:225
Bool_t firstStage() const
Check whether the parameter should be floated only in the first stage of a two stage fit...
The abstract interface to the fitter.
Definition: LauAbsFitter.hh:35
virtual Bool_t twoStageFit() const
Determine whether the two-stage fit is enabled.
Definition: LauMinuit.hh:59
TMatrixD covMatrix_
The covariance matrix.
Definition: LauMinuit.hh:149
File containing declaration of LauParameter class.
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:140
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:134
File containing declaration of LauFitObject class.
TVirtualFitter * minuit_
The interface to Minuit.
Definition: LauMinuit.hh:122
virtual void fixFirstStageParameters()
Fix parameters marked as &quot;first stage&quot;.
Definition: LauMinuit.cc:212
Double_t value() const
The value of the parameter.
virtual void fixSecondStageParameters()
Fix parameters marked as &quot;second stage&quot;.
Definition: LauMinuit.cc:238
virtual void withinAsymErrorCalc(Bool_t inAsymErrCalc)=0
Mark that the fit is calculating asymmetric errors.
virtual LauFitObject * getFitObject()
Get the object that controls the calculation of the likelihood.
Definition: LauMinuit.cc:124
std::vector< LauParameter * > params_
The fit parameters.
Definition: LauMinuit.hh:128
File containing declaration of LauMinuit class.
The abstract interface for the objects that control the calculation of the likelihood.
Definition: LauFitObject.hh:26