UnbinnedLikelihoodFit.cc

Go to the documentation of this file.

#include "RecoTracker/DeDx/interface/UnbinnedLikelihoodFit.h"
#include <TMath.h>

// a class to perform a likelihood fit
// Author: Christophe Delaere

/* Example of a Landau fit:
 * ------------------------
 * UnbinnedLikelihoodFit myfit;
 * double x[4] = {89,110,70,80};
 * myfit.setData(4,x);
 * TF1* myfunction = new TF1("myLandau","TMath::Landau(x,[0],[1],1)",0,255);
 * myfunction->SetParameters(100,10);
 * myfit.setFunction(myfunction);
 * myfit.fit();
 * myfit.getFunction()->Print();
 * double MPV = myfit.getFunction()->GetParameter(0);
 * double MPVerror = myfit.getFunction()->GetParError(0);
 */

// the function passed to minuit
void UnbinnedLL(Int_t&, Double_t*, Double_t &val, Double_t *par, Int_t) {
  // retrieve the data object (it's also the fitter)
  // - sign to have a minimum
  // factor 2 to have the right errors (see for example the pdg)
  val = -2*((dynamic_cast<const UnbinnedLikelihoodFit*>((TVirtualFitter::GetFitter())->GetObjectFit()))->logL(par));
}

// the constructor
UnbinnedLikelihoodFit::UnbinnedLikelihoodFit() {
  nparameters_ = 0;
  datasize_ = 0;
  x_ = nullptr;
  min = nullptr;
  tolerance_ = 0.01;
  maxIterations_ = 1000;
}

// the destructor
UnbinnedLikelihoodFit::~UnbinnedLikelihoodFit() {
}

// sets the data
// the class is not owner of the data... it only keeps a pointer to it.
void UnbinnedLikelihoodFit::setData(uint32_t n, double* x) {
  datasize_ = n;
  x_ = x;
}

// sets the function for the fit
void UnbinnedLikelihoodFit::setFunction(TF1* f) {
  function_ = f;
  nparameters_ = function_ ? function_->GetNpar() : 0;
}

// The fit itself
int32_t UnbinnedLikelihoodFit::fit(int32_t verbosity) {
  // creates a fitter 
  min = TVirtualFitter::Fitter(this,nparameters_);
  min->SetFCN(UnbinnedLL);
  
  // set print level: no output
  arglist_[0] = 0;
  min->ExecuteCommand("SET NOWarnings",arglist_,1);
  arglist_[0] = verbosity;
  min->ExecuteCommand("SET PRINT",arglist_,1);
  
  // initial values, error, range
  double parmin,parmax;
  for(uint32_t i=0;i<nparameters_;++i) {
    function_->GetParLimits(i, parmin, parmax);
    min->SetParameter(i,
                      function_->GetParName(i),
                      function_->GetParameter(i),
                      tolerance_,
                      parmin, parmax);
  }

  // run MIGRAD
  arglist_[0] = maxIterations_; // number of function calls
  arglist_[1] = tolerance_;     // tolerance
  int32_t status = min->ExecuteCommand("MIGRAD",arglist_,2);

  // get fit parameters and errors
  for(uint32_t i=0;i<nparameters_;++i) {
    function_->SetParameter(i, min->GetParameter(i));
    function_->SetParError( i, min->GetParError(i) );
  }

  // returns the status
  return status;
}

// the log-likelihood function
double UnbinnedLikelihoodFit::logL(const double* parameters) const {
  double val=0;
  if(!function_) return val;
  for (uint32_t i=0;i<datasize_;++i){
    val += TMath::Log(function_->EvalPar(&(x_[i]),parameters));
  }
  return val; 
}