DeDx/src/UnbinnedLikelihoodFit.cc

0001 #include "RecoTracker/DeDx/interface/UnbinnedLikelihoodFit.h"
0002 #include <TMath.h>
0003
0004 // a class to perform a likelihood fit
0005 // Author: Christophe Delaere
0006
0007 /* Example of a Landau fit:
0008  * ------------------------
0009  * UnbinnedLikelihoodFit myfit;
0010  * double x[4] = {89,110,70,80};
0011  * myfit.setData(4,x);
0012  * TF1* myfunction = new TF1("myLandau","TMath::Landau(x,[0],[1],1)",0,255);
0013  * myfunction->SetParameters(100,10);
0014  * myfit.setFunction(myfunction);
0015  * myfit.fit();
0016  * myfit.getFunction()->Print();
0017  * double MPV = myfit.getFunction()->GetParameter(0);
0018  * double MPVerror = myfit.getFunction()->GetParError(0);
0019  */
0020
0021 // the function passed to minuit
0022 void UnbinnedLL(Int_t&, Double_t*, Double_t& val, Double_t* par, Int_t) {
0023   // retrieve the data object (it's also the fitter)
0024   // - sign to have a minimum
0025   // factor 2 to have the right errors (see for example the pdg)
0026   val = -2 * ((dynamic_cast<const UnbinnedLikelihoodFit*>((TVirtualFitter::GetFitter())->GetObjectFit()))->logL(par));
0027 }
0028
0029 // the constructor
0030 UnbinnedLikelihoodFit::UnbinnedLikelihoodFit() {
0031   nparameters_ = 0;
0032   datasize_ = 0;
0033   x_ = nullptr;
0034   min = nullptr;
0035   tolerance_ = 0.01;
0036   maxIterations_ = 1000;
0037 }
0038
0039 // the destructor
0040 UnbinnedLikelihoodFit::~UnbinnedLikelihoodFit() {}
0041
0042 // sets the data
0043 // the class is not owner of the data... it only keeps a pointer to it.
0044 void UnbinnedLikelihoodFit::setData(uint32_t n, double* x) {
0045   datasize_ = n;
0046   x_ = x;
0047 }
0048
0049 // sets the function for the fit
0050 void UnbinnedLikelihoodFit::setFunction(TF1* f) {
0051   function_ = f;
0052   nparameters_ = function_ ? function_->GetNpar() : 0;
0053 }
0054
0055 // The fit itself
0056 int32_t UnbinnedLikelihoodFit::fit(int32_t verbosity) {
0057   // creates a fitter
0058   min = TVirtualFitter::Fitter(this, nparameters_);
0059   min->SetFCN(UnbinnedLL);
0060
0061   // set print level: no output
0062   arglist_[0] = 0;
0063   min->ExecuteCommand("SET NOWarnings", arglist_, 1);
0064   arglist_[0] = verbosity;
0065   min->ExecuteCommand("SET PRINT", arglist_, 1);
0066
0067   // initial values, error, range
0068   double parmin, parmax;
0069   for (uint32_t i = 0; i < nparameters_; ++i) {
0070     function_->GetParLimits(i, parmin, parmax);
0071     min->SetParameter(i, function_->GetParName(i), function_->GetParameter(i), tolerance_, parmin, parmax);
0072   }
0073
0074   // run MIGRAD
0075   arglist_[0] = maxIterations_;  // number of function calls
0076   arglist_[1] = tolerance_;      // tolerance
0077   int32_t status = min->ExecuteCommand("MIGRAD", arglist_, 2);
0078
0079   // get fit parameters and errors
0080   for (uint32_t i = 0; i < nparameters_; ++i) {
0081     function_->SetParameter(i, min->GetParameter(i));
0082     function_->SetParError(i, min->GetParError(i));
0083   }
0084
0085   // returns the status
0086   return status;
0087 }
0088
0089 // the log-likelihood function
0090 double UnbinnedLikelihoodFit::logL(const double* parameters) const {
0091   double val = 0;
0092   if (!function_)
0093     return val;
0094   for (uint32_t i = 0; i < datasize_; ++i) {
0095     val += TMath::Log(function_->EvalPar(&(x_[i]), parameters));
0096   }
0097   return val;
0098 }