Regressor function object for JPoint4D fit using JSimplex minimiser. More...

#include <JPoint4DRegressor.hh>

Inheritance diagram for JFIT::JRegressor< JPoint4D, JSimplex >:

Public Types
typedef JSimplex< JPoint4D >	minimiser_type

typedef JRegressor< JPoint4D, JSimplex >	regressor_type

Public Member Functions
	JRegressor (const bool __prefit)
	Constructor. More...

	JRegressor (TH2D *hPDF2D, bool __prefit)
	Constructor for a 2D histogram as PDF. More...

template<class JHit_t >
double	operator() (const JPoint4D &vertex, const JHit_t &hit) const

double	evaluatePrefit (const double t_res, const double cosT, const double cosT_th) const

double	evaluateFit (const double t_res, const double D) const

double	operator() (const JPoint4D &value, T __begin, T __end)
	Global fit. More...

double	operator() (const JPoint4D &value, T1 __begin1, T1 __end1, T2 __begin2, T2 __end2)
	Global fit. More...

double	operator() (const JFunction_t &fit, T __begin, T __end)
	Multi-dimensional fit. More...

double	operator() (const JFunction_t &fit, T __begin, T __end, const JPoint4D &step)
	1D fit. More...

Public Attributes
TH2D *	hpdf2d [NUMBER_OF_PDFS]
	PDF 2D. More...

JLANG::JSharedPointer < JMEstimator >	estimator
	M-Estimator function. More...

const bool	prefit

JPoint4D	value

std::vector< JPoint4D >	step

int	numberOfIterations

Static Public Attributes
static const double	a = 2.0
	constant to prevent the likelihood going to infinity More...

static const double	d_ref = 10.0

static const double	p_max = 100

static const int	NUMBER_OF_PDFS = 1

static int	MAXIMUM_ITERATIONS
	maximal number of iterations More...

static double	EPSILON
	maximal distance to minimum More...

static int	debug = 0
	debug level (default is off). More...

Detailed Description

template<>
struct JFIT::JRegressor< JPoint4D, JSimplex >

Regressor function object for JPoint4D fit using JSimplex minimiser.

Definition at line 47 of file JPoint4DRegressor.hh.

Member Typedef Documentation

typedef JSimplex <JPoint4D > JFIT::JAbstractRegressor< JPoint4D , JSimplex >::minimiser_type

inherited

Definition at line 76 of file JRegressor.hh.

typedef JRegressor<JPoint4D , JSimplex > JFIT::JAbstractRegressor< JPoint4D , JSimplex >::regressor_type

inherited

Definition at line 77 of file JRegressor.hh.

Constructor & Destructor Documentation

JFIT::JRegressor< JPoint4D, JSimplex >::JRegressor ( const bool __prefit )

inline

Constructor.

Definition at line 55 of file JPoint4DRegressor.hh.

                                    :
       estimator(new JMEstimatorTukey(30.0)), 
       prefit(__prefit)
     {}

JFIT::JRegressor< JPoint4D, JSimplex >::JRegressor	(	TH2D *	hPDF2D,
		bool	__prefit
	)

inline

Constructor for a 2D histogram as PDF.

Definition at line 64 of file JPoint4DRegressor.hh.

                                            :
       estimator(new JMEstimatorTukey(30.0)),
       prefit(__prefit)
     {
       hpdf2d[0] = (TH2D*)hPDF2D->Clone();
       hpdf2d[0]->Smooth();
     }

Member Function Documentation

template<class JHit_t >

double JFIT::JRegressor< JPoint4D, JSimplex >::operator()	(	const JPoint4D &	vertex,
		const JHit_t &	hit
	)		const

inline

Definition at line 82 of file JPoint4DRegressor.hh.

     {
       using namespace std;
       using namespace JGEOMETRY3D;
       using namespace JTOOLS;
       
       const JPosition3D hit_pos(hit.getPosition());
       const double D = hit_pos.getDistance(vertex);
       const double t_res = hit.getT() - vertex.getT(hit_pos);
       const JDirection3D photonDir(hit_pos - vertex.getPosition());
       const double cosT = photonDir.getDot(hit.getDirection());
       const double cosT_th = d_ref / sqrt(d_ref*d_ref + D*D);
 
       double LogLik = 0;
       
       if(prefit){
         
         LogLik = evaluatePrefit(t_res, cosT, cosT_th);
         
       } else {
         
         LogLik = evaluateFit(t_res, D);
         
       }
 
       return LogLik;
     }

double JFIT::JRegressor< JPoint4D, JSimplex >::evaluatePrefit	(	const double	t_res,
		const double	cosT,
		const double	cosT_th
	)		const

inline

Definition at line 111 of file JPoint4DRegressor.hh.

     {
 
       double p_T = 0;
 
       if(cosT < 0){
         p_T = 0;
       } else if(cosT > cosT_th){
         p_T = p_max;
       } else {
         p_T = p_max * cosT/cosT_th;
       }
       
       double lik = 1 / (sqrt(a*a + t_res*t_res) + p_T);
       double u   = -log(lik);
       
       return estimator->getRho(u);      
     }

double JFIT::JRegressor< JPoint4D, JSimplex >::evaluateFit	(	const double	t_res,
		const double	D
	)		const

inline

Definition at line 131 of file JPoint4DRegressor.hh.

     {
      
       int x_bin = hpdf2d[0]->GetXaxis()->FindBin(t_res);
       int y_bin = hpdf2d[0]->GetYaxis()->FindBin(D);
       
       int nbinsX = hpdf2d[0]->GetXaxis()->GetNbins();
       int nbinsY = hpdf2d[0]->GetYaxis()->GetNbins();
       
       double lik = 0;
       
       if(x_bin > 1 && x_bin < nbinsX){
         if(y_bin > 1 && y_bin < nbinsY){
           lik = hpdf2d[0]->Interpolate(t_res, D);
         }
       }
       
       if(lik == 0){
         lik = 0.000004;  // value of K40 bg                               
       }
 
       double u   = -log(lik);
       return estimator->getRho(u);
       
     }

double JFIT::JAbstractRegressor< JPoint4D , JSimplex >::operator()	(	const JPoint4D &	value,
		T	__begin,
		T	__end
	)

inlineinherited

Global fit.

Parameters

value	start value
__begin	begin of data set
__end	end of data set

Returns: chi2

Definition at line 89 of file JRegressor.hh.

     {
       static_cast<minimiser_type&>(*this).value = value;
 
       return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin, __end);
     }

double JFIT::JAbstractRegressor< JPoint4D , JSimplex >::operator()	(	const JPoint4D &	value,
		T1	__begin1,
		T1	__end1,
		T2	__begin2,
		T2	__end2
	)

inlineinherited

Global fit.

Parameters

value	start value
__begin1	begin of first data set
__end1	end of first data set
__begin2	begin of second data set
__end2	end of second data set

Returns: chi2

Definition at line 108 of file JRegressor.hh.

     {
       static_cast<minimiser_type&>(*this).value = value;
 
       return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin1, __end1, __begin2, __end2);
     }

double JFIT::JSimplex< JPoint4D >::operator()	(	const JFunction_t &	fit,
		T	__begin,
		T	__end
	)

inlineinherited

Multi-dimensional fit.

The given fit function should return the equivalent of chi2 for the current value of the given model and a given data point.

Parameters

fit	fit function
__begin	begin of data
__end	end of data

Returns: chi2

Definition at line 69 of file JSimplex.hh.

     {
       using namespace std;
       using namespace JPP;
 
       double chi2_old = evaluate(fit, __begin, __end);
 
       const int N = step.size();
 
       if (N != 0) {
 
         double chi2[N];
 
         numberOfIterations = 0;
 
         while (numberOfIterations != MAXIMUM_ITERATIONS) {
 
           DEBUG("old: " << FIXED(12,5) << chi2_old << endl);
 
           const JModel_t p0(value);
 
           for (int i = 0; i != N; ++i) {
 
             DEBUG("step: " << i << ' ' << setw(5) << numberOfIterations << endl);
 
             chi2[i] = (*this)(fit, __begin, __end, step[i]);
           }
 
           // overall step direction of last iteration
 
           JModel_t wall = value - p0;
 
           const double chi2_new = (*this)(fit, __begin, __end, wall);
 
           DEBUG("new: " << FIXED(12,5) << chi2_new << endl);
 
 
           if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
             return chi2_new;
           }
 
           // double overall step
 
           wall   = value - p0;
 
           value += wall;
 
           const double fe = evaluate(fit, __begin, __end);
 
           value -= wall;
 
 
           for (int i = N-1; i != 0; --i) {
             chi2[i] = chi2[i-1] - chi2[i];
           }
 
           chi2[0] = chi2_old - chi2[0];
 
 
           double df = 0.0;
 
           for (int i = 0; i != N; ++i) {
             if (chi2[i] > df) {
               df = chi2[i];
             }
           }
 
           const double fn = chi2_new;
           const double f0 = chi2_old;
           const double ff = f0 - fn - df;
 
           // shift step directions
 
           if (fe < f0 && 2.0*(f0 - 2.0*fn + fe)*ff*ff < (f0-fe)*(f0-fe)*df) {
 
             for (int i = 0; i != N - 1; ++i) {
               step[i] = step[i+1];
             }
 
             step[N-1] = wall;
           }
 
           chi2_old = chi2_new;
         }
       }
 
       return chi2_old;
     }

double JFIT::JSimplex< JPoint4D >::operator()	(	const JFunction_t &	fit,
		T	__begin,
		T	__end,
		const JPoint4D &	step
	)

inlineinherited

1D fit.

The given fit function should return the equivalent of chi2 for the current value of the given model and a given data point.

Parameters

fit	fit function
__begin	begin of data
__end	end of data
step	step direction

Definition at line 171 of file JSimplex.hh.

     {
       using namespace std;
       using namespace JPP;
 
       double lambda = 0.5;                                     // control parameter
       double factor = 1.0;                                     // multiplication factor
 
       double chi2_old = evaluate(fit, __begin, __end);
 
       for (int i = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations, ++i) {
 
         value += lambda * step;
 
         const double chi2_new = evaluate(fit, __begin, __end);
 
         DEBUG("step: " << setw(3) << i << ' ' << FIXED(12,5) << chi2_old << ' ' << FIXED(12,5) << chi2_new << ' ' << FIXED(5,2) << lambda << endl);
 
         if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
 
           if (chi2_new > chi2_old) {
 
             value -= lambda * step;                            // undo last step
 
             return chi2_old;
 
           } else {
 
             return chi2_new;
           }
         }
 
         if (chi2_new < chi2_old) {
 
           chi2_old = chi2_new;
 
         } else {
             
           value  -=  lambda * step;                            // step back
           lambda  = -lambda;                                   // change direction
           
           if (i != 0) {
             factor = 0.5;                                      // reduce step size
           }
         }
         
         lambda = factor * lambda;
       }
 
       return chi2_old;
     }