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

#include <JShower3EZRegressor.hh>

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

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::JPolint0FunctionalMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList < JTOOLS::JPolint0FunctionalGridMap > > > >	JNPEMaplist_t

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

typedef JSimplex< JShower3EZ >	minimiser_type

typedef JRegressor< JShower3EZ, JSimplex >	regressor_type

typedef minimiser_type::result_type	result_type

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

double	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 cosDelta, 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...

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 JShower3EZ &step)
	1D fit. More...

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

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

JShower3EZ	value

std::vector< JShower3EZ >	step

int	numberOfIterations

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 JPDFType_t	pdf_t [NUMBER_OF_PDFS]
	PDF types. More...

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< JShower3EZ, JSimplex >

Regressor function object for JShower3EZ fit using JSimplex minimiser.

Definition at line 52 of file JShower3EZRegressor.hh.

Member Typedef Documentation

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

Definition at line 57 of file JShower3EZRegressor.hh.

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

Definition at line 61 of file JShower3EZRegressor.hh.

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

Definition at line 62 of file JShower3EZRegressor.hh.

typedef JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList<JTOOLS::JPolint0FunctionalGridMap> > > > JFIT::JRegressor< JShower3EZ, JSimplex >::JNPEMaplist_t

Definition at line 67 of file JShower3EZRegressor.hh.

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

Definition at line 68 of file JShower3EZRegressor.hh.

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

inherited

Definition at line 78 of file JRegressor.hh.

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

inherited

Definition at line 79 of file JRegressor.hh.

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

inherited

Definition at line 80 of file JRegressor.hh.

Constructor & Destructor Documentation

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

inline

Parameterized constructor.

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

Parameters

fileDescriptor PDF file descriptor

Definition at line 79 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

double JFIT::JRegressor< JShower3EZ, JSimplex >::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 129 of file JShower3EZRegressor.hh.

     {
       using namespace JPP;
       
       JPosition3D D(pmt.getPosition());
       JDirection3D U(pmt.getDirection());
 
       D.sub(shower.getPosition());
 
       JVersor3D shower_dir(shower.getDirection());
 
       const double z  = D.getDot(shower_dir);
       const double x  = D.getX()  -  z * shower.getDX();
       const double y  = D.getY()  -  z * shower.getDY();
       const double cosDelta = z/D.getLength();  // Delta = angle between shower direction and PMT position
  
       U.rotate(JRotation3Z(-atan2(y,x)));                  // rotate PMT axis to x-z plane
 
       const double theta = U.getTheta();
       const double phi   = fabs(U.getPhi());
 
       JNPE_t::result_type H0 = getH0(pmt.getR());  // getH0 = Get background hypothesis value for time integrated PDF. 
       JNPE_t::result_type H1 = getH1(D.getLength(), cosDelta, theta, phi, shower.getE());  // getH1 = Get signal hypothesis value for time integrated PDF. 
 
       if (get_value(H1) >= Vmax_npe) {
         H1 *= Vmax_npe / get_value(H1);
       }
 
       H1 += H0; // now H1 is signal + background
 
       const bool hit = pmt.getN() != 0;
       const double u = getChi2(get_value(H1), hit); // -log(lik)
       
       return estimator->getRho(u);
     }

JNPE_t::result_type JFIT::JRegressor< JShower3EZ, JSimplex >::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 171 of file JShower3EZRegressor.hh.

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

JNPE_t::result_type JFIT::JRegressor< JShower3EZ, JSimplex >::getH1	(	const double	D,
		const double	cosDelta,
		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]
cosDelta	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 186 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()), cosDelta, theta, phi);
 
             // safety measure
 
             if(y1 < 0){
               y1 = 0;  
             }
             
             h1 += y1;
           
           }
           catch(JLANG::JException& error) {
             ERROR(error << std::endl);
           }
         }
       }
 
       return h1;
     }

result_type JFIT::JAbstractRegressor< JShower3EZ , JSimplex >::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 92 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::JSimplex< JShower3EZ >::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 71 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) {
 
         p0   = value;
         p1   = value;
         wall = value;
 
         double chi2[N];
 
         numberOfIterations = 0;
 
         while (numberOfIterations != MAXIMUM_ITERATIONS) {
 
           DEBUG("old: " << FIXED(12,5) << chi2_old << endl);
 
           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
 
           wall  = value;
           wall -= 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;
           wall  -= 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< JShower3EZ >::operator()	(	const JFunction_t &	fit,
		T	__begin,
		T	__end,
		const JShower3EZ &	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 178 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) {
 
         p1      = step;
         p1     *= lambda;
         value  += p1;
 
         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) {
 
             p1      = step;
             p1     *= lambda;
             value  -= p1;                                      // undo last step
 
             return chi2_old;
 
           } else {
 
             return chi2_new;
           }
         }
 
         if (chi2_new < chi2_old) {
 
           chi2_old = chi2_new;
 
         } else {
             
           p1      = step;
           p1     *= lambda;
           value  -=  p1;                                       // step back
           lambda  = -lambda;                                   // change direction
           
           if (i != 0) {
             factor = 0.5;                                      // reduce step size
           }
         }
         
         lambda = factor * lambda;
       }
 
       return chi2_old;
     }