Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Types | Public Member Functions | Public Attributes | Static Public Attributes | List of all members
JFIT::JRegressor< JLine3Z, JGandalf > Struct Template Reference

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

#include <JLine3ZRegressor.hh>

Inheritance diagram for JFIT::JRegressor< JLine3Z, JGandalf >:
JFIT::JAbstractRegressor< JLine3Z, JGandalf > JFIT::JGandalf< JLine3Z > JEEP::JMessage< T >

Public Types

typedef JTOOLS::JSplineFunction1S_t JFunction1D_t
 
typedef JTOOLS::JMAPLIST
< JTOOLS::JPolint1FunctionalMap,
JTOOLS::JPolint0FunctionalGridMap,
JTOOLS::JPolint0FunctionalGridMap >
::maplist 
JPDFMaplist_t
 
typedef JPHYSICS::JPDFTable
< JFunction1D_t, JPDFMaplist_t
JPDF_t
 
typedef JTOOLS::JMAPLIST
< JTOOLS::JPolint1FunctionalMapH,
JTOOLS::JPolint1FunctionalGridMap,
JTOOLS::JPolint1FunctionalGridMap >
::maplist 
JNPEMaplist_t
 
typedef JPHYSICS::JNPETable
< double, double,
JNPEMaplist_t
JNPE_t
 
typedef JGandalf< JLine3Zminimiser_type
 
typedef JRegressor< JLine3Z,
JGandalf
regressor_type
 
typedef JLine3Z::parameter_type parameter_type
 Data type of fit parameter. More...
 

Public Member Functions

 JRegressor ()
 Default constructor. More...
 
 JRegressor (const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
 Constructor. More...
 
template<class JHit_t >
result_type operator() (const JLine3Z &track, const JHit_t &hit) const
 Fit function. More...
 
result_type operator() (const JLine3Z &track, const JPMTW0 &pmt) const
 Fit function. More...
 
JPDF_t::result_type getH0 (const double R_Hz, const double t1) const
 Get background hypothesis value for time differentiated PDF. More...
 
JPDF_t::result_type getH1 (const double E, const double R, const double theta, const double phi, const double t1) const
 Get signal hypothesis value for time differentiated PDF. 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 E, const double R, const double theta, const double phi) const
 Get signal hypothesis value for time integrated PDF. More...
 
double getRmax () const
 Get maximal road width of PDF. More...
 
double operator() (const JLine3Z &value, T __begin, T __end)
 Global fit. More...
 
double operator() (const JLine3Z &value, T1 __begin1, T1 __end1, T2 __begin2, T2 __end2)
 Global fit. More...
 
double operator() (const JFunction_t &fit, T1 __begin1, T1 __end1, T2 __begin2, T2 __end2)
 Multi-dimensional fit of two data sets. More...
 
double operator() (const JFunction_t &fit, T __begin, T __end)
 Multi-dimensional fit of one data set. More...
 

Public Attributes

JPDF_t pdf [NUMBER_OF_PDFS]
 PDF. More...
 
JNPE_t npe [NUMBER_OF_PDFS]
 PDF. More...
 
double E_GeV
 Energy of muon at vertex [GeV]. More...
 
JLANG::JSharedPointer
< JMEstimator
estimator
 M-Estimator function. More...
 
double lambda
 
JLine3Z value
 
JLine3Z error
 
std::vector< parameter_typeparameters
 
int numberOfIterations
 
JMATH::JMatrixNS H
 

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 = 4
 
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 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 matrix More...
 
static int debug = 0
 debug level (default is off). More...
 

Detailed Description

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

Regressor function object for JLine3Z fit using JGandalf minimiser.

Definition at line 110 of file JLine3ZRegressor.hh.

Member Typedef Documentation

Definition at line 115 of file JLine3ZRegressor.hh.

Definition at line 118 of file JLine3ZRegressor.hh.

Definition at line 119 of file JLine3ZRegressor.hh.

Definition at line 123 of file JLine3ZRegressor.hh.

Definition at line 124 of file JLine3ZRegressor.hh.

Definition at line 76 of file JRegressor.hh.

Definition at line 77 of file JRegressor.hh.

Data type of fit parameter.

Definition at line 56 of file JGandalf.hh.

Constructor & Destructor Documentation

Default constructor.

Definition at line 129 of file JLine3ZRegressor.hh.

130  {};
JFIT::JRegressor< JLine3Z, JGandalf >::JRegressor ( const std::string &  fileDescriptor,
const double  TTS,
const int  numberOfPoints = 25,
const double  epsilon = 1.0e-10 
)
inline

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<JLine3Z, JGandalf>::pdf_t.

The TTS corresponds to the additional time smearing applied to the PDFs.

Parameters
fileDescriptorPDF file descriptor
TTSTTS [ns]
numberOfPointsnumber of points for Gauss-Hermite integration of TTS
epsilonprecision for Gauss-Hermite integration of TTS

Definition at line 146 of file JLine3ZRegressor.hh.

149  :
150  E_GeV(0.0),
152  {
153  using namespace std;
154  using namespace JPHYSICS;
155 
156  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
157 
158  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
159 
160  try {
161 
162  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
163 
164  NOTICE("loading PDF from file " << file_name << "... " << flush);
165 
166  pdf[i].load(file_name.c_str());
167 
168  NOTICE("OK" << endl);
169 
170  pdf[i].setExceptionHandler(supervisor);
171 
172  if (TTS > 0.0) {
173  pdf[i].blur(TTS, numberOfPoints, epsilon);
174  } else if (TTS < 0.0) {
175  ERROR("Illegal value of TTS [ns]: " << TTS << endl);
176  }
177  }
178  catch(const JException& error) {
179  FATAL(error.what() << endl);
180  }
181  }
182 
183  // Add PDFs
184 
185  for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
186 
187  pdf[ i ].add(pdf[i-1]);
188  pdf[i-1].clear();
189 
190  npe[ i ] = JNPE_t(pdf[i]);
191  npe[ i ].transform(JNPE_t::transformer_type::getDefaultTransformer()); // get rid of R-dependent weight function
192  }
193  }
void load(const char *file_name)
Load from input file.
void setExceptionHandler(const supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
void blur(const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10, const double quantile=0.99)
Blur PDF.
Definition: JPDFTable.hh:108
virtual void transform(const transformer_type &transformer)
Application of weight function.
Definition: JNPETable.hh:194
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
void add(const JMultiFunction_t &input)
Add function.
double E_GeV
Energy of muon at vertex [GeV].
#define NOTICE(A)
Definition: JMessage.hh:62
Normal M-estimator.
Definition: JMEstimator.hh:66
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
#define FATAL(A)
Definition: JMessage.hh:65
int numberOfPoints
Definition: JResultPDF.cc:22
std::string getFilename(const std::string &file_name)
Get file name part, i.e.
Definition: JeepToolkit.hh:85
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.

Member Function Documentation

template<class JHit_t >
result_type JFIT::JRegressor< JLine3Z, JGandalf >::operator() ( const JLine3Z track,
const JHit_t &  hit 
) const
inline

Fit function.

This method is used to determine the chi2 and gradient of given hit with respect to trajectory of muon.

The template argument JHit_t refers to a data structure which should have the following member methods:

  • double getX(); // [m]
  • double getY(); // [m]
  • double getZ(); // [m]
  • double getDX(); // [u]
  • double getDY(); // [u]
  • double getDZ(); // [u]
  • double getT(); // [ns]
Parameters
tracktrack
hithit
Returns
chi2 and gradient

Definition at line 214 of file JLine3ZRegressor.hh.

215  {
216  using namespace JGEOMETRY3D;
217  using namespace JTOOLS;
218  using namespace JPHYSICS;
219 
220  JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
221  JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
222 
223  D.sub(track.getPosition());
224 
225  const double z = D.getDot(track.getDirection());
226  const double x = D.getX() - z * track.getDX();
227  const double y = D.getY() - z * track.getDY();
228  const double R = sqrt(D.getLengthSquared() - z*z);
229 
230  const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
231 
232  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
233 
234  const double theta = U.getTheta();
235  const double phi = fabs(U.getPhi());
236 
237  const double E = gWater.getE(E_GeV, z);
238  const double dt = T_ns.constrain(hit.getT() - t1);
239 
240  JPDF_t::result_type H0 = getH0(hit.getR(), dt);
241  JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
242 
243  if (H1.V >= Vmax_npe) {
244  H1 *= Vmax_npe / H1.V;
245  }
246 
247  H1 += H0; // signal + background
248 
249  result_type result;
250 
251  result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio
252 
253  const double wc = 1.0 - getTanThetaC() * z / R; // d(ct)/d(z)
254 
255  result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R, // d(ct)/d(x)
256  -getTanThetaC() * D.getY() / R, // d(ct)/d(y)
257  0.0),
258  getSpeedOfLight()), // d(ct)/d(t)
259  JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()), // d(ct)/d(dx)
260  wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ()))); // d(ct)/d(dy)
261 
262  result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
263  H0.getDerivativeOfChi2())); // x d(chi2)/d(ct)
264 
265  return result;
266  }
static const JGeane gWater(2.67e-1 *JTOOLS::DENSITY_SEA_WATER, 3.4e-4 *JTOOLS::DENSITY_SEA_WATER)
Function object for Energy loss of muon in sea water.
double getE(const double E, const double dx) const
Get energy of muon after specified distance.
Definition: JGeane.hh:81
JPDF_t::result_type getH0(const double R_Hz, const double t1) const
Get background hypothesis value for time differentiated PDF.
JPDF_t::result_type getH1(const double E, const double R, const double theta, const double phi, const double t1) const
Get signal hypothesis value for time differentiated PDF.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
const double getSpeedOfLight()
Number of bytes in a gigabyte.
Definition: JConstants.hh:89
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JLine3Z.hh:233
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:35
double getDot(const JAngle3D &angle) const
Get dot product.
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:155
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Rotation around Z-axis.
Definition: JRotation3D.hh:82
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition: JLine3Z.hh:254
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
JLine3Z & mul(const double value)
Multiplication operator.
Definition: JLine3Z.hh:160
static double Vmax_npe
Maximal integral of PDF [npe].
double getDX() const
Get x direction.
Definition: JVersor3Z.hh:144
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
double getDZ() const
Get z direction.
Definition: JVersor3Z.hh:166
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:122
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:36
result_type JFIT::JRegressor< JLine3Z, JGandalf >::operator() ( const JLine3Z track,
const JPMTW0 pmt 
) const
inline

Fit function.

This method is used to determine the chi2 and gradient of given PMT with respect to trajectory of muon.

Parameters
tracktrack
pmtpmt
Returns
chi2 and gradient

Definition at line 277 of file JLine3ZRegressor.hh.

278  {
279  using namespace JGEOMETRY3D;
280  using namespace JPHYSICS;
281 
282  JPosition3D D(pmt.getPosition());
283  JDirection3D U(pmt.getDirection());
284 
285  D.sub(track.getPosition());
286 
287  const double z = D.getDot(track.getDirection());
288  const double x = D.getX() - z * track.getDX();
289  const double y = D.getY() - z * track.getDY();
290  const double R = sqrt(D.getLengthSquared() - z*z);
291 
292  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
293 
294  const double theta = U.getTheta();
295  const double phi = fabs(U.getPhi());
296 
297  const double E = gWater.getE(E_GeV, z);
298 
299  JNPE_t::result_type H0 = getH0(pmt.getR());
300  JNPE_t::result_type H1 = getH1(E, R, theta, phi);
301 
302  if (H1.f >= Vmax_npe) {
303  H1 *= Vmax_npe / H1.f;
304  }
305 
306  H1 += H0;
307 
308  const bool hit = pmt.getN() != 0;
309  const double u = H1.getChi2(hit);
310 
311  result_type result;
312 
313  result.chi2 = estimator->getRho(u);
314 
315  result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R, // d(R)/d(x)
316  -D.getY() / R, // d(R)/d(y)
317  0.0),
318  0.0),
319  JVersor3Z(-z * D.getX() / R, // d(R)/d(dx)
320  -z * D.getY() / R)); // d(R)/d(dy)
321 
322  result.gradient.mul(estimator->getPsi(u));
323  result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(R)
324 
325  return result;
326  }
static const JGeane gWater(2.67e-1 *JTOOLS::DENSITY_SEA_WATER, 3.4e-4 *JTOOLS::DENSITY_SEA_WATER)
Function object for Energy loss of muon in sea water.
double getE(const double E, const double dx) const
Get energy of muon after specified distance.
Definition: JGeane.hh:81
double getR() const
Get rate.
Definition: JPMTW0.hh:56
JPDF_t::result_type getH0(const double R_Hz, const double t1) const
Get background hypothesis value for time differentiated PDF.
JPDF_t::result_type getH1(const double E, const double R, const double theta, const double phi, const double t1) const
Get signal hypothesis value for time differentiated PDF.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
const JDirection3D & getDirection() const
Get direction.
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:35
double getDot(const JAngle3D &angle) const
Get dot product.
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:155
Rotation around Z-axis.
Definition: JRotation3D.hh:82
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition: JLine3Z.hh:254
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
int getN() const
Get number of hits.
Definition: JPMTW0.hh:67
static double Vmax_npe
Maximal integral of PDF [npe].
double getDX() const
Get x direction.
Definition: JVersor3Z.hh:144
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:36
JPDF_t::result_type JFIT::JRegressor< JLine3Z, JGandalf >::getH0 ( const double  R_Hz,
const double  t1 
) const
inline

Get background hypothesis value for time differentiated PDF.

Parameters
R_Hzrate [Hz]
t1time [ns]
Returns
hypothesis value

Definition at line 336 of file JLine3ZRegressor.hh.

338  {
339  return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
340  }
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
JPDF_t::result_type JFIT::JRegressor< JLine3Z, JGandalf >::getH1 ( const double  E,
const double  R,
const double  theta,
const double  phi,
const double  t1 
) const
inline

Get signal hypothesis value for time differentiated PDF.

Parameters
Emuon energy at minimum distance of approach [GeV]
Rminimum distance of approach [m]
thetaPMT zenith angle [rad]
phiPMT azimuth angle [rad]
t1arrival time relative to Cherenkov hypothesis [ns]
Returns
hypothesis value

Definition at line 353 of file JLine3ZRegressor.hh.

358  {
359  using namespace JPHYSICS;
360 
361  JPDF_t::result_type h1 = JMATH::zero;
362 
363  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
364 
365  if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
366 
367  JPDF_t::result_type y1 = pdf[i](std::max(R, pdf[i].getXmin()), theta, phi, t1);
368 
369  // safety measures
370 
371  if (y1.f <= 0.0) {
372  y1.f = 0.0;
373  y1.fp = 0.0;
374  }
375 
376  if (y1.v <= 0.0) {
377  y1.v = 0.0;
378  }
379 
380  // energy dependence
381 
382  if (is_deltarays(pdf_t[i])) {
383  y1 *= getDeltaRaysFromMuon(E);
384  } else if (is_bremsstrahlung(pdf_t[i])) {
385  y1 *= E;
386  }
387 
388  h1 += y1;
389  }
390  }
391 
392  return h1;
393  }
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:157
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:171
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:81
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.
JNPE_t::result_type JFIT::JRegressor< JLine3Z, JGandalf >::getH0 ( const double  R_Hz) const
inline

Get background hypothesis value for time integrated PDF.

Parameters
R_Hzrate [Hz]
Returns
hypothesis value

Definition at line 402 of file JLine3ZRegressor.hh.

403  {
404  return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
405  }
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
JNPE_t::result_type JFIT::JRegressor< JLine3Z, JGandalf >::getH1 ( const double  E,
const double  R,
const double  theta,
const double  phi 
) const
inline

Get signal hypothesis value for time integrated PDF.

Parameters
Emuon energy at minimum distance of approach [GeV]
Rminimum distance of approach [m]
thetaPMT zenith angle [rad]
phiPMT azimuth angle [rad]
Returns
hypothesis value

Definition at line 417 of file JLine3ZRegressor.hh.

421  {
422  using namespace JPHYSICS;
423 
424  JNPE_t::result_type h1 = JMATH::zero;
425 
426  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
427 
428  if (!npe[i].empty() && R <= npe[i].getXmax()) {
429 
430  try {
431 
432  JNPE_t::result_type y1 = npe[i](std::max(R, npe[i].getXmin()), theta, phi);
433 
434  // safety measures
435 
436  if (y1.f > 0.0) {
437 
438  // energy dependence
439 
440  if (is_deltarays(pdf_t[i])) {
441  y1 *= getDeltaRaysFromMuon(E);
442  } else if (is_bremsstrahlung(pdf_t[i])) {
443  y1 *= E;
444  }
445 
446  h1 += y1;
447  }
448  }
449  catch(JLANG::JException& error) {
450  ERROR(error << std::endl);
451  }
452  }
453  }
454 
455  return h1;
456  }
General exception.
Definition: JException.hh:40
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:157
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:171
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:81
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.
double JFIT::JRegressor< JLine3Z, JGandalf >::getRmax ( ) const
inline

Get maximal road width of PDF.

Returns
road width [m]

Definition at line 464 of file JLine3ZRegressor.hh.

465  {
466  double xmax = 0.0;
467 
468  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
469  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
470  xmax = pdf[i].getXmax();
471  }
472  }
473 
474  return xmax;
475  }
double JFIT::JAbstractRegressor< JLine3Z , JGandalf >::operator() ( const JLine3Z value,
__begin,
__end 
)
inlineinherited

Global fit.

Parameters
valuestart value
__beginbegin of data set
__endend of data set
Returns
chi2

Definition at line 89 of file JRegressor.hh.

90  {
91  static_cast<minimiser_type&>(*this).value = value;
92 
93  return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin, __end);
94  }
JModel_t value
Definition: JGandalf.hh:265
JRegressor< JLine3Z, JGandalf > regressor_type
Definition: JRegressor.hh:77
double JFIT::JAbstractRegressor< JLine3Z , JGandalf >::operator() ( const JLine3Z value,
T1  __begin1,
T1  __end1,
T2  __begin2,
T2  __end2 
)
inlineinherited

Global fit.

Parameters
valuestart value
__begin1begin of first data set
__end1end of first data set
__begin2begin of second data set
__end2end of second data set
Returns
chi2

Definition at line 108 of file JRegressor.hh.

111  {
112  static_cast<minimiser_type&>(*this).value = value;
113 
114  return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin1, __end1, __begin2, __end2);
115  }
JModel_t value
Definition: JGandalf.hh:265
JRegressor< JLine3Z, JGandalf > regressor_type
Definition: JRegressor.hh:77
double JFIT::JGandalf< JLine3Z >::operator() ( const JFunction_t &  fit,
T1  __begin1,
T1  __end1,
T2  __begin2,
T2  __end2 
)
inlineinherited

Multi-dimensional fit of two data sets.

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

Parameters
fitfit function
__begin1begin of first data set
__end1end of first data set
__begin2begin of second data set
__end2end of second data set
Returns
chi2

Definition at line 122 of file JGandalf.hh.

125  {
126  using namespace std;
127 
128  const int N = parameters.size();
129 
130  H.resize(N);
131  h.resize(N);
132 
133  previous = value;
134  error = JModel_t();
135  lambda = LAMBDA_MIN;
136 
137  double chi2_old = numeric_limits<double>::max();
138 
140 
141  DEBUG("step: " << numberOfIterations << endl);
142 
143  reset();
144 
145  evaluate(fit, __begin1, __end1);
146  evaluate(fit, __begin2, __end2);
147 
148  DEBUG("lambda: " << SCIENTIFIC(12,5) << lambda << endl);
149  DEBUG("chi2: " << FIXED (12,5) << chi2 << endl);
150 
151  if (chi2 < chi2_old) {
152 
153  if (numberOfIterations != 0) {
154 
155  if (fabs(chi2_old - chi2) < EPSILON*fabs(chi2_old)) {
156  return chi2;
157  }
158 
159  lambda /= LAMBDA_DOWN;
160  }
161 
162  chi2_old = chi2;
163  previous = value;
164 
165  } else {
166 
167  value = previous;
168  lambda *= LAMBDA_UP;
169 
170  if (lambda > LAMBDA_MAX) {
171  return chi2_old; // no improvement found
172  }
173 
174  reset();
175 
176  evaluate(fit, __begin1, __end1);
177  evaluate(fit, __begin2, __end2);
178  }
179 
180  DEBUG("Hesse matrix:" << endl);
181  DEBUG(H << endl);
182 
183 
184  // force definite positiveness
185 
186  for (int i = 0; i != N; ++i) {
187 
188  if (H(i,i) < PIVOT) {
189  H(i,i) = PIVOT;
190  }
191 
192  h[i] = 1.0 / sqrt(H(i,i));
193  }
194 
195 
196  // normalisation
197 
198  for (int i = 0; i != N; ++i) {
199  for (int j = 0; j != i; ++j) {
200  H(j,i) *= h[i] * h[j];
201  H(i,j) = H(j,i);
202  }
203  }
204 
205  for (int i = 0; i != N; ++i) {
206  H(i,i) = 1.0 + lambda;
207  }
208 
209 
210  try {
211  H.invert();
212  }
213  catch (const JException& error) {
214  ERROR("JGandalf: " << error.what() << endl);
215  return chi2_old;
216  }
217 
218 
219  for (int i = 0; i != N; ++i) {
220 
221  DEBUG("u[" << noshowpos << i << "] = " << showpos << FIXED(15,5) << value.*parameters[i]);
222 
223  for (int j = 0; j != N; ++j) {
224  value.*parameters[i] -= H(i,j) * gradient.*parameters[j] * h[i] * h[j];
225  }
226 
227  DEBUG(' ' << FIXED(15,5) << value.*parameters[i] << noshowpos << endl);
228 
229  error.*parameters[i] = h[i];
230  }
231  }
232 
233  return chi2_old;
234  }
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:256
static double PIVOT
minimal value diagonal element of matrix
Definition: JGandalf.hh:261
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
std::vector< parameter_type > parameters
Definition: JGandalf.hh:267
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:257
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:260
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:98
void evaluate(const JFunction_t &fit, T __begin, T __end)
Evaluate fit for given data set.
Definition: JGandalf.hh:292
#define ERROR(A)
Definition: JMessage.hh:64
JMATH::JMatrixNS H
Definition: JGandalf.hh:269
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:259
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:255
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:258
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:498
void invert()
Invert matrix.
Definition: JMatrixNS.hh:61
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
std::vector< double > h
Definition: JGandalf.hh:312
double JFIT::JGandalf< JLine3Z >::operator() ( const JFunction_t &  fit,
__begin,
__end 
)
inlineinherited

Multi-dimensional fit of one data set.

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

Parameters
fitfit function
__beginbegin of data
__endend of data
Returns
chi2

Definition at line 249 of file JGandalf.hh.

250  {
251  return (*this)(fit, __begin, __end, __end, __end);
252  }

Member Data Documentation

JTimeRange JFIT::JRegressor< JLine3Z, JGandalf >::T_ns
static

Time window with respect to Cherenkov hypothesis [ns].

Default values.

Definition at line 478 of file JLine3ZRegressor.hh.

double JFIT::JRegressor< JLine3Z, JGandalf >::Vmax_npe = std::numeric_limits<double>::max()
static

Maximal integral of PDF [npe].

Definition at line 479 of file JLine3ZRegressor.hh.

const int JFIT::JRegressor< JLine3Z, JGandalf >::NUMBER_OF_PDFS = 4
static

Definition at line 481 of file JLine3ZRegressor.hh.

const JPDFType_t JFIT::JRegressor< JLine3Z, JGandalf >::pdf_t
static

PDF.

Definition at line 485 of file JLine3ZRegressor.hh.

PDF.

Definition at line 486 of file JLine3ZRegressor.hh.

double JFIT::JRegressor< JLine3Z, JGandalf >::E_GeV

Energy of muon at vertex [GeV].

Definition at line 487 of file JLine3ZRegressor.hh.

M-Estimator function.

Definition at line 489 of file JLine3ZRegressor.hh.

int JFIT::JGandalf< JLine3Z >::MAXIMUM_ITERATIONS
staticinherited

maximal number of iterations

maximal number of iterations.

Definition at line 255 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::EPSILON
staticinherited

maximal distance to minimum

maximal distance to minimum.

Definition at line 256 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::LAMBDA_MIN
staticinherited

minimal value control parameter

Definition at line 257 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::LAMBDA_MAX
staticinherited

maximal value control parameter

Definition at line 258 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::LAMBDA_UP
staticinherited

multiplication factor control parameter

Definition at line 259 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::LAMBDA_DOWN
staticinherited

multiplication factor control parameter

Definition at line 260 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::PIVOT
staticinherited

minimal value diagonal element of matrix

Definition at line 261 of file JGandalf.hh.

double JFIT::JGandalf< JLine3Z >::lambda
inherited

Definition at line 264 of file JGandalf.hh.

JLine3Z JFIT::JGandalf< JLine3Z >::value
inherited

Definition at line 265 of file JGandalf.hh.

JLine3Z JFIT::JGandalf< JLine3Z >::error
inherited

Definition at line 266 of file JGandalf.hh.

std::vector<parameter_type> JFIT::JGandalf< JLine3Z >::parameters
inherited

Definition at line 267 of file JGandalf.hh.

int JFIT::JGandalf< JLine3Z >::numberOfIterations
inherited

Definition at line 268 of file JGandalf.hh.

Definition at line 269 of file JGandalf.hh.

template<class T>
int JEEP::JMessage< T >::debug = 0
staticinherited

debug level (default is off).

Definition at line 43 of file JMessage.hh.


The documentation for this struct was generated from the following file: