Regressor function object for JShower3EZ fit using JGandalf minimiser. More...

#include <JShower3EZRegressor.hh>

Inheritance diagram for JFIT::JRegressor< JShower3EZ, JGandalf >:

Public Types
typedef JTOOLS::JSplineFunction1S_t	JFunction1D_t

typedef JTOOLS::JMapList < JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalGridMap > > > >	JPDFMaplist_t

typedef JPHYSICS::JPDFTable < JFunction1D_t, JPDFMaplist_t >	JPDF_t

typedef JTOOLS::JMapList < JTOOLS::JPolint1FunctionalMap, JTOOLS::JMapList < JTOOLS::JPolint1FunctionalMapH, JTOOLS::JMapList < JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JMapList < JTOOLS::JPolint1FunctionalGridMap > > > >	JNPEMaplist_t

typedef JPHYSICS::JNPETable < double, double, JNPEMaplist_t >	JNPE_t

typedef JGandalf< JShower3EZ >	minimiser_type

typedef JRegressor< JShower3EZ, JGandalf >	regressor_type

typedef minimiser_type::result_type	result_type

typedef JFIT_LOCAL::JTypedef_t < JShower3EZ >::parameter_type	parameter_type
	Data type of fit parameter. More...

Public Member Functions
	JRegressor (const std::string &fileDescriptor)
	Parameterized constructor. More...

result_type	operator() (const JShower3EZ &shower, const JPMTW0 &pmt) const
	Fit function. More...

JNPE_t::result_type	getH0 (const double R_Hz) const
	Get background hypothesis value for time integrated PDF. More...

JNPE_t::result_type	getH1 (const double D, const double cd, const double theta, const double phi, const double E) const
	Get signal hypothesis value for time integrated PDF. More...

result_type	operator() (const JShower3EZ &value, T __begin, T __end)
	Global fit. More...

result_type	operator() (const JFunction_t &fit, T __begin, T __end, Args...args)
	Multi-dimensional fit of multiple data sets. More...

Public Attributes
JNPE_t	npe [NUMBER_OF_PDFS]
	PDF. More...

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

std::vector< parameter_type >	parameters
	fit parameters More...

int	numberOfIterations
	number of iterations More...

double	lambda
	control parameter More...

JShower3EZ	value
	value More...

JShower3EZ	error
	error More...

JMATH::JMatrixNS	V
	Hesse matrix. More...

result_type	result

Static Public Attributes
static JTimeRange	T_ns
	Time window with respect to Cherenkov hypothesis [ns]. More...

static double	Vmax_npe = std::numeric_limits<double>::max()
	Maximal integral of PDF [npe]. More...

static const int	NUMBER_OF_PDFS = 2

static const JFIT::JPDFType_t	pdf_t [NUMBER_OF_PDFS]

static int	MAXIMUM_ITERATIONS
	maximal number of iterations More...

static double	EPSILON
	maximal distance to minimum More...

static double	LAMBDA_MIN
	minimal value control parameter More...

static double	LAMBDA_MAX
	maximal value control parameter More...

static double	LAMBDA_UP
	multiplication factor control parameter More...

static double	LAMBDA_DOWN
	multiplication factor control parameter More...

static double	PIVOT
	minimal value diagonal element of Hesse matrix More...

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

Detailed Description

template<>
struct JFIT::JRegressor< JShower3EZ, JGandalf >

Regressor function object for JShower3EZ fit using JGandalf minimiser.

Definition at line 276 of file JShower3EZRegressor.hh.

Member Typedef Documentation

typedef JTOOLS::JSplineFunction1S_t JFIT::JRegressor< JShower3EZ, JGandalf >::JFunction1D_t

Definition at line 281 of file JShower3EZRegressor.hh.

typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > > JFIT::JRegressor< JShower3EZ, JGandalf >::JPDFMaplist_t

Definition at line 285 of file JShower3EZRegressor.hh.

typedef JPHYSICS::JPDFTable<JFunction1D_t, JPDFMaplist_t> JFIT::JRegressor< JShower3EZ, JGandalf >::JPDF_t

Definition at line 286 of file JShower3EZRegressor.hh.

typedef JTOOLS::JMapList<JTOOLS::JPolint1FunctionalMap, JTOOLS::JMapList<JTOOLS::JPolint1FunctionalMapH, JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap> > > > JFIT::JRegressor< JShower3EZ, JGandalf >::JNPEMaplist_t

Definition at line 291 of file JShower3EZRegressor.hh.

typedef JPHYSICS::JNPETable<double, double, JNPEMaplist_t> JFIT::JRegressor< JShower3EZ, JGandalf >::JNPE_t

Definition at line 292 of file JShower3EZRegressor.hh.

typedef JGandalf <JShower3EZ > JFIT::JAbstractRegressor< JShower3EZ , JGandalf >::minimiser_type

inherited

Definition at line 80 of file JRegressor.hh.

typedef JRegressor<JShower3EZ , JGandalf > JFIT::JAbstractRegressor< JShower3EZ , JGandalf >::regressor_type

inherited

Definition at line 81 of file JRegressor.hh.

typedef minimiser_type::result_type JFIT::JAbstractRegressor< JShower3EZ , JGandalf >::result_type

inherited

Definition at line 82 of file JRegressor.hh.

typedef JFIT_LOCAL::JTypedef_t<JShower3EZ >::parameter_type JFIT::JGandalf< JShower3EZ >::parameter_type

inherited

Data type of fit parameter.

Definition at line 95 of file JGandalf.hh.

Constructor & Destructor Documentation

JFIT::JRegressor< JShower3EZ, JGandalf >::JRegressor ( const std::string & fileDescriptor )

inline

Parameterized constructor.

The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which will be replaced by the corresponding PDF types listed in JRegressor<JShower3Z, JGandalf>::pdf_t.

Parameters

fileDescriptor PDF file descriptor

Definition at line 303 of file JShower3EZRegressor.hh.

                                                :
       estimator(new JMEstimatorNull())
     {
       using namespace std;
       using namespace JPP;
 
       const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
       
       for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
         
         try {
             
           JPDF_t pdf;
 
           const string file_name = getFilename(fileDescriptor, pdf_t[i]);
             
           NOTICE("loading PDF from file " << file_name << "... " << flush);
 
           pdf.load(file_name.c_str());
           
           NOTICE("OK" << endl);
 
           pdf.setExceptionHandler(supervisor);
 
           npe[i] = JNPE_t(pdf);
         }
         catch(const JException& error) {
           FATAL(error.what() << endl);
         }
       }
 
       // Add PDFs
       for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
 
         npe[ i ].add(npe[i-1]);
 
         JNPE_t buffer;
 
         npe[i-1].swap(buffer);
       }    
     }

Member Function Documentation

result_type JFIT::JRegressor< JShower3EZ, JGandalf >::operator()	(	const JShower3EZ &	shower,
		const JPMTW0 &	pmt
	)		const

inline

Fit function.

This method is used to determine the chi2 of given PMT with respect to shower hypothesis.

Parameters

shower	shower
pmt	pmt

Returns: chi2

Definition at line 353 of file JShower3EZRegressor.hh.

     {
       using namespace JPP;
       using namespace std;
 
       JPosition3D  D(pmt.getPosition());
       JDirection3D U(pmt.getDirection());
 
       D.sub(shower.getPosition());
 
       const double x  = D.getX();
       const double y  = D.getY();
       const double d  = D.getLength();  
       const double cd = D.getDot(shower.getDirection())/d;  // cosine angle between shower direction and PMT position
 
       U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
       const double theta = getPMTAngle(U.getTheta());
       const double phi   = getPMTAngle(U.getPhi());
 
       JNPE_t::result_type H0 = getH0(pmt.getR());                        // background hypothesis value for time integrated PDF. 
       JNPE_t::result_type H1 = getH1(d, cd, theta, phi, shower.getE());  // signal hypothesis value for time integrated PDF. 
 
       if (get_value(H1) >= Vmax_npe) {
         H1 *= Vmax_npe / get_value(H1);
       }
 
       double signal_npe = get_value(H1);
 
       H1 += H0; // now H1 is signal + background
 
       double expectation = get_value(H1);
 
       const bool hit = pmt.getN() != 0;
       const double u = H1.getChi2(hit);
 
       result_type result;
 
       result.chi2     = estimator->getRho(u);
 
       double energy_gradient = signal_npe/shower.getE(); // dP/dE
       if(hit) energy_gradient *= -exp(-expectation)/(1-exp(-expectation)); //dchi2/d(H1), if !hit is 1
 
       result.gradient = JShower3EZ(JPoint4D(JVector3D(0,       // d(cos_th0)/d(x)
                                                       0,       // d(cos_th0)/d(y)
                                                       0),          // d(cos_th0)/d(z) 
                                             0.0),                                 // d(cos_th0)/d(t) 
                                    JVersor3Z(x/d,                // d(cos_th0)/d(dx)
                                              y/d),               // d(cos_th0)/d(dy)  
                                    energy_gradient);                                          // d(chi2)/d(E)
      
       result.gradient.mul(estimator->getPsi(u));
       static_cast<JShower3Z&>(result.gradient).mul(H1.getDerivativeOfChi2(hit));                    // x d(chi2)/d(cos_th0)
       
       return result;
     }

JNPE_t::result_type JFIT::JRegressor< JShower3EZ, JGandalf >::getH0 ( const double R_Hz ) const

inline

Get background hypothesis value for time integrated PDF.

Parameters

R_Hz rate [Hz]

Returns: hypothesis value

Definition at line 416 of file JShower3EZRegressor.hh.

     {
       return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
     }

JNPE_t::result_type JFIT::JRegressor< JShower3EZ, JGandalf >::getH1	(	const double	D,
		const double	cd,
		const double	theta,
		const double	phi,
		const double	E
	)		const

inline

Get signal hypothesis value for time integrated PDF.

Parameters

D	PMT distance from shower [m]
cd	cosine angle between shower direction and PMT position
theta	PMT zenith angle [deg]
phi	PMT azimuth angle [deg]
E	shower energy [GeV]

Returns: hypothesis value

Definition at line 431 of file JShower3EZRegressor.hh.

     {
       JNPE_t::result_type h1 = JMATH::zero;
         
       for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
 
         if (!npe[i].empty() && D <= npe[i].getXmax()) {
 
           try {
 
             JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi);
           
             if (get_value(y1) > 0.0) {
               h1 += y1;
             }
           
           }
           catch(JLANG::JException& error) {
             ERROR(error << std::endl);
           }
         }
       }
 
       return h1;
     }

result_type JFIT::JAbstractRegressor< JShower3EZ , JGandalf >::operator()	(	const JShower3EZ &	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 94 of file JRegressor.hh.

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

result_type JFIT::JGandalf< JShower3EZ >::operator()	(	const JFunction_t &	fit,
		T	__begin,
		T	__end,
		Args...	args
	)

inlineinherited

Multi-dimensional fit of multiple data sets.

The fit function should return the chi2 as well as the partial derivatives for the current value of the model and a given data point.

Parameters

fit	fit function
__begin	begin of data
__end	end of data
args	optional data

Returns: chi2 and gradient

Definition at line 160 of file JGandalf.hh.

     {
       using namespace std;
       using namespace JPP;
 
       // note that all model values should be assigned to the start value of the model before use
       // because the actual list of model parameters can vary from fit to fit
       // (e.g. if model consists of a container).
 
       const size_t N = parameters.size();
 
       V.resize(N);
       h.resize(N);
       x.resize(N);
 
       previous.result.chi2 = numeric_limits<double>::max();
 
       current.result.chi2     = numeric_limits<double>::max();
       current.result.gradient = value;
       current.result.gradient = zero;
 
       error     = value;
       error     = zero;
 
       lambda    = LAMBDA_MIN;
 
       for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
 
         DEBUG("step:     " << numberOfIterations << endl);
 
         reset();
 
         update(fit, __begin, __end, args...);
 
         DEBUG("lambda:   " << FIXED(12,5) << lambda              << endl);
         DEBUG("chi2:     " << FIXED(12,5) << current.result.chi2 << endl);
 
         if (current.result.chi2 < previous.result.chi2) {
 
           if (numberOfIterations != 0) {
 
             if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
 
               // normal end
 
               const result_type result = current.result;
 
               lambda  = LAMBDA_MIN;
 
               reset();
 
               update(fit, __begin, __end, args...);
 
               try {
                 V.invert();
               }
               catch (const exception& error) {
                 V.reset();
               }
 
               for (size_t i = 0; i != N; ++i) {
                 get(error, parameters[i]) = sqrt(V(i,i));
               }
 
               return result;
             }
 
             if (lambda > LAMBDA_MIN) {
               lambda /= LAMBDA_DOWN;
             }
           }
 
           // store current values
 
           previous.value  = value;
           previous.result = current.result;
 
         } else {
 
           value   = previous.value;                                 // restore value
 
           lambda *= LAMBDA_UP;
           
           if (lambda > LAMBDA_MAX) {
             break;
           }
 
           reset();
 
           update(fit, __begin, __end, args...);
         }
 
         DEBUG("Hesse matrix:" << endl << V << endl);
 
         // force definite positiveness
 
         for (size_t i = 0; i != N; ++i) {
 
           if (V(i,i) < PIVOT) {
             V(i,i) = PIVOT;
           }
 
           h[i] = 1.0 / sqrt(V(i,i));
         }
 
         // normalisation
         
         for (size_t row = 0; row != N; ++row) {
           for (size_t col = 0; col != row; ++col) {
             V(row,col) *= h[row] * h[col];
             V(col,row)  = V(row,col);
           }
         }
         
         for (size_t i = 0; i != N; ++i) {
           V(i,i) = 1.0 + lambda;
         }
 
         // solve A x = b
 
         for (size_t col = 0; col != N; ++col) {
           x[col] = h[col] * get(current.result.gradient, parameters[col]);
         }
 
         try {
           V.solve(x);
         }
         catch (const exception& error) {
 
           ERROR("JGandalf: " << error.what() << endl << V << endl);
 
           break;
         }
 
         // update value
 
         for (size_t row = 0; row != N; ++row) {
 
           DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
 
           get(value, parameters[row]) -= h[row] * x[row];
 
           DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
         }
 
         model(value);
       }
 
       // abnormal end
 
       const result_type result = previous.result;
 
       value   = previous.value;                                     // restore value
 
       lambda  = LAMBDA_MIN;
 
       reset();
 
       update(fit, __begin, __end, args...);
 
       try {
         V.invert();
       }
       catch (const exception& error) {
         V.reset();
       }
 
       for (size_t i = 0; i != N; ++i) {
         get(error, parameters[i]) = sqrt(V(i,i));
       }
 
       return result;
     }