Jpp  18.5.0
the software that should make you happy
 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::JRegressorStorage< JLine3Z, JGandalf > JFIT::JGandalf< JLine3Z > JFIT::JRegressorHelper< JLine3Z, JGandalf > JEEP::JMessage< T >

Public Types

typedef JRegressorStorage
< JLine3Z, JGandalf
JRegressorStorage_t
 
typedef JGandalf< JLine3Zminimiser_type
 
typedef JRegressor< JLine3Z,
JGandalf
regressor_type
 
typedef minimiser_type::result_type result_type
 
typedef JLine3Z::parameter_type parameter_type
 Data type of fit parameter. More...
 
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
 time dependent PDF More...
 
typedef JTOOLS::JMAPLIST
< JTOOLS::JPolint1FunctionalMapH,
JTOOLS::JPolint1FunctionalGridMap,
JTOOLS::JPolint1FunctionalGridMap >
::maplist 
JNPEMaplist_t
 
typedef JPHYSICS::JNPETable
< double, double,
JNPEMaplist_t
JNPE_t
 time integrated PDF More...
 
typedef JPDF_t::transformer_type transformer_type
 
typedef std::array< JPDF_t,
NUMBER_OF_PDFS
JPDFs_t
 PDFs. More...
 
typedef std::array< JNPE_t,
NUMBER_OF_PDFS
JNPEs_t
 NPEs. 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...
 
 JRegressor (const JRegressorStorage< JLine3Z, JGandalf > &storage)
 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...
 
result_type operator() (const JLine3Z &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...
 
const JPDFs_tgetPDF () const
 Get PDFs. More...
 
const JNPEs_tgetNPE () const
 Get NPEs. More...
 
void transform (const transformer_type &transformer)
 Transform PDFs and NPEs. More...
 
void setRmax (const double Rmax)
 Set maximal road width of PDF. More...
 

Public Attributes

const JPDFs_tpdf
 PDF. More...
 
const JNPEs_tnpe
 PDF. More...
 
double E_GeV = 0.0
 Energy of muon at vertex [GeV]. More...
 
JLANG::JSharedPointer
< JMEstimator
estimator = new JMEstimatorNormal()
 M-Estimator function. More...
 
std::vector< parameter_typeparameters
 fit parameters More...
 
int numberOfIterations
 number of iterations More...
 
double lambda
 control parameter More...
 
JLine3Z value
 value More...
 
JLine3Z 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 double Rmin_m = 0.1
 Minimal distance of [m]. 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 Hesse matrix More...
 
static int debug = 0
 debug level (default is off). More...
 
static const int NUMBER_OF_PDFS = 6
 Number of PDFs. More...
 
static const JPDFType_t pdf_t [NUMBER_OF_PDFS]
 PDF types. More...
 

Detailed Description

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

Regressor function object for JLine3Z fit using JGandalf minimiser.

Definition at line 103 of file JLine3ZRegressor.hh.

Member Typedef Documentation

Definition at line 109 of file JLine3ZRegressor.hh.

Definition at line 80 of file JRegressor.hh.

Definition at line 81 of file JRegressor.hh.

typedef minimiser_type::result_type JFIT::JAbstractRegressor< JLine3Z , JGandalf >::result_type
inherited

Definition at line 82 of file JRegressor.hh.

Data type of fit parameter.

Definition at line 64 of file JGandalf.hh.

Definition at line 57 of file JRegressorHelper.hh.

Definition at line 60 of file JRegressorHelper.hh.

time dependent PDF

Definition at line 61 of file JRegressorHelper.hh.

Definition at line 65 of file JRegressorHelper.hh.

time integrated PDF

Definition at line 66 of file JRegressorHelper.hh.

Definition at line 68 of file JRegressorHelper.hh.

PDFs.

Definition at line 72 of file JRegressorHelper.hh.

NPEs.

Definition at line 73 of file JRegressorHelper.hh.

Constructor & Destructor Documentation

Default constructor.

Definition at line 114 of file JLine3ZRegressor.hh.

114  :
116  pdf(getPDF()),
117  npe(getNPE())
118  {}
const JNPEs_t & getNPE() const
Get NPEs.
const JPDFs_t & getPDF() const
Get PDFs.
JRegressorStorage< JLine3Z, JGandalf > JRegressorStorage_t
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::WILDCARD 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 134 of file JLine3ZRegressor.hh.

137  :
138  JRegressorStorage_t(fileDescriptor, TTS, numberOfPoints, epsilon),
139  pdf(getPDF()),
140  npe(getNPE())
141  {}
const JNPEs_t & getNPE() const
Get NPEs.
const JPDFs_t & getPDF() const
Get PDFs.
int numberOfPoints
Definition: JResultPDF.cc:22
JRegressorStorage< JLine3Z, JGandalf > JRegressorStorage_t
const double epsilon
Definition: JQuadrature.cc:21

Constructor.

Parameters
storagePDF storage

Definition at line 149 of file JLine3ZRegressor.hh.

149  :
150  pdf(storage.getPDF()),
151  npe(storage.getNPE())
152  {}
const JNPEs_t & getNPE() const
Get NPEs.
const JPDFs_t & getPDF() const
Get PDFs.

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 173 of file JLine3ZRegressor.hh.

174  {
175  using namespace JPP;
176 
177  JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
178  JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
179 
180  D.sub(track.getPosition());
181 
182  const double z = D.getDot(track.getDirection());
183  const double x = D.getX() - z * track.getDX();
184  const double y = D.getY() - z * track.getDY();
185  const double R2 = D.getLengthSquared() - z*z;
186  const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
187 
188  const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
189 
190  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
191 
192  const double theta = U.getTheta();
193  const double phi = fabs(U.getPhi());
194 
195  const double E = gWater.getE(E_GeV, z);
196  const double dt = T_ns.constrain(hit.getT() - t1);
197 
198  JPDF_t::result_type H0 = getH0(hit.getR(), dt);
199  JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
200 
201  if (H1.V >= Vmax_npe) {
202  H1 *= Vmax_npe / H1.V;
203  }
204 
205  H1 += H0; // signal + background
206 
208 
209  result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio
210 
211  const double wc = 1.0 - getTanThetaC() * z / R; // d(ct1)/d(z)
212 
213  result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R, // d(ct1)/d(x)
214  -getTanThetaC() * D.getY() / R, // d(ct1)/d(y)
215  0.0),
216  getSpeedOfLight()), // d(ct1)/d(t)
217  JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()), // d(ct1)/d(dx)
218  wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ()))); // d(ct1)/d(dy)
219 
220  result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
221  H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1)
222 
223  return result;
224  }
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.
JMATH::JVectorND x
Definition: JGandalf.hh:445
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JLine3Z.hh:234
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:36
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
static double Rmin_m
Minimal distance of [m].
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:158
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Rotation around Z-axis.
Definition: JRotation3D.hh:85
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition: JLine3Z.hh:255
minimiser_type::result_type result_type
Definition: JRegressor.hh:82
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
virtual double getE(const double E, const double dx) const override
Get energy of muon after specified distance.
Definition: JGeane.hh:257
const double getSpeedOfLight()
Get speed of light.
static double Vmax_npe
Maximal integral of PDF [npe].
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getDX() const
Get x direction.
Definition: JVersor3Z.hh:147
Data structure for fit of straight line paralel to z-axis.
Definition: JLine1Z.hh:27
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
double getDZ() const
Get z direction.
Definition: JVersor3Z.hh:169
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
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 235 of file JLine3ZRegressor.hh.

236  {
237  using namespace JGEOMETRY3D;
238  using namespace JPHYSICS;
239 
240  JPosition3D D(pmt.getPosition());
241  JDirection3D U(pmt.getDirection());
242 
243  D.sub(track.getPosition());
244 
245  const double z = D.getDot(track.getDirection());
246  const double x = D.getX() - z * track.getDX();
247  const double y = D.getY() - z * track.getDY();
248  const double R2 = D.getLengthSquared() - z*z;
249  const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
250 
251  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
252 
253  const double theta = U.getTheta();
254  const double phi = fabs(U.getPhi());
255 
256  const double E = gWater.getE(E_GeV, z);
257 
258  JNPE_t::result_type H0 = getH0(pmt.getR());
259  JNPE_t::result_type H1 = getH1(E, R, theta, phi);
260 
261  if (H1.f >= Vmax_npe) {
262  H1 *= Vmax_npe / H1.f;
263  }
264 
265  H1 += H0;
266 
267  const bool hit = pmt.getN() != 0;
268  const double u = H1.getChi2(hit);
269 
271 
272  result.chi2 = estimator->getRho(u);
273 
274  result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R, // d(R)/d(x)
275  -D.getY() / R, // d(R)/d(y)
276  0.0),
277  0.0),
278  JVersor3Z(-z * D.getX() / R, // d(R)/d(dx)
279  -z * D.getY() / R)); // d(R)/d(dy)
280 
281  result.gradient.mul(estimator->getPsi(u));
282  result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(R)
283 
284  return result;
285  }
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.
JMATH::JVectorND x
Definition: JGandalf.hh:445
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:33
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
const JDirection3D & getDirection() const
Get direction.
Data structure for fit of straight line in positive z-direction.
Definition: JLine3Z.hh:36
static const JGeaneWater gWater
Function object for energy loss of muon in sea water.
Definition: JGeane.hh:381
static double Rmin_m
Minimal distance of [m].
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:158
Rotation around Z-axis.
Definition: JRotation3D.hh:85
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition: JLine3Z.hh:255
minimiser_type::result_type result_type
Definition: JRegressor.hh:82
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
virtual double getE(const double E, const double dx) const override
Get energy of muon after specified distance.
Definition: JGeane.hh:257
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:147
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:36
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
double u[N+1]
Definition: JPolint.hh:865
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
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 295 of file JLine3ZRegressor.hh.

297  {
298  return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
299  }
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 312 of file JLine3ZRegressor.hh.

317  {
318  using namespace std;
319  using namespace JPP;
320 
321  JPDF_t::result_type h1 = zero;
322 
323  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
324 
325  if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
326 
327  JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
328 
329  // safety measures
330 
331  if (y1.f <= 0.0) {
332  y1.f = 0.0;
333  y1.fp = 0.0;
334  }
335 
336  if (y1.v <= 0.0) {
337  y1.v = 0.0;
338  }
339 
340  // energy dependence
341 
342  if (is_deltarays(pdf_t[i])) {
343  y1 *= getDeltaRaysFromMuon(E);
344  } else if (is_bremsstrahlung(pdf_t[i])) {
345  y1 *= E;
346  }
347 
348  h1 += y1;
349  }
350  }
351 
352  return h1;
353  }
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:137
static const int NUMBER_OF_PDFS
Number of PDFs.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:151
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
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 362 of file JLine3ZRegressor.hh.

363  {
364  return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
365  }
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 377 of file JLine3ZRegressor.hh.

381  {
382  using namespace std;
383  using namespace JPP;
384 
385  JNPE_t::result_type h1 = JMATH::zero;
386 
387  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
388 
389  if (!npe[i].empty() && R <= npe[i].getXmax()) {
390 
391  try {
392 
393  JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
394 
395  // safety measures
396 
397  if (y1.f > 0.0) {
398 
399  // energy dependence
400 
401  if (is_deltarays(pdf_t[i])) {
402  y1 *= getDeltaRaysFromMuon(E);
403  } else if (is_bremsstrahlung(pdf_t[i])) {
404  y1 *= E;
405  }
406 
407  h1 += y1;
408  }
409  }
410  catch(JException& error) {
411  ERROR(error << endl);
412  }
413  }
414  }
415 
416  return h1;
417  }
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:137
static const int NUMBER_OF_PDFS
Number of PDFs.
#define ERROR(A)
Definition: JMessage.hh:66
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:151
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:67
double JFIT::JRegressor< JLine3Z, JGandalf >::getRmax ( ) const
inline

Get maximal road width of PDF.

Returns
road width [m]

Definition at line 425 of file JLine3ZRegressor.hh.

426  {
427  double xmax = 0.0;
428 
429  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
430  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
431  xmax = pdf[i].getXmax();
432  }
433  }
434 
435  return xmax;
436  }
const double xmax
Definition: JQuadrature.cc:24
static const int NUMBER_OF_PDFS
Number of PDFs.
result_type JFIT::JAbstractRegressor< JLine3Z , JGandalf >::operator() ( const JLine3Z value,
T  __begin,
T  __end 
)
inlineinherited

Global fit.

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

Definition at line 94 of file JRegressor.hh.

95  {
96  static_cast<minimiser_type&>(*this).value = value;
97 
98  return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin, __end);
99  }
JModel_t value
value
Definition: JGandalf.hh:303
JRegressor< JLine3Z, JGandalf > regressor_type
Definition: JRegressor.hh:81
result_type JFIT::JGandalf< JLine3Z >::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
fitfit function
__beginbegin of data
__endend of data
argsoptional data
Returns
chi2 and gradient

Definition at line 129 of file JGandalf.hh.

130  {
131  using namespace std;
132  using namespace JPP;
133 
134  // note that all model values should be assigned to the start value of the model before use
135  // because the actual list of model parameters can vary from fit to fit
136  // (e.g. if model consists of a container).
137 
138  const size_t N = parameters.size();
139 
140  V.resize(N);
141  h.resize(N);
142  x.resize(N);
143 
144  previous.result.chi2 = numeric_limits<double>::max();
145 
146  current.result.chi2 = numeric_limits<double>::max();
147  current.result.gradient = value;
148  current.result.gradient = zero;
149 
150  error = value;
151  error = zero;
152 
153  lambda = LAMBDA_MIN;
154 
155 
157 
158  DEBUG("step: " << numberOfIterations << endl);
159 
160  reset();
161 
162  update(fit, __begin, __end, args...);
163 
164  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
165  DEBUG("chi2: " << FIXED(12,5) << current.result.chi2 << endl);
166 
167  if (current.result.chi2 < previous.result.chi2) {
168 
169  if (numberOfIterations != 0) {
170 
171  if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
172 
173  const result_type result = current.result;
174 
175  reset();
176 
177  update(fit, __begin, __end, args...);
178 
179  try {
180  V.invert();
181  }
182  catch (const exception& error) {
183  V.reset();
184  }
185 
186  return result; // normal end
187  }
188 
189  if (lambda > LAMBDA_MIN) {
190  lambda /= LAMBDA_DOWN;
191  }
192 
193  // store current values
194 
195  previous.value = value;
196  previous.error = error;
197  previous.result = current.result;
198  }
199 
200  } else {
201 
202  value = previous.value; // restore value
203 
204  lambda *= LAMBDA_UP;
205 
206  if (lambda > LAMBDA_MAX) {
207  break;
208  }
209 
210  reset();
211 
212  update(fit, __begin, __end, args...);
213  }
214 
215  DEBUG("Hesse matrix:" << endl << V << endl);
216 
217  // force definite positiveness
218 
219  for (size_t i = 0; i != N; ++i) {
220 
221  if (V(i,i) < PIVOT) {
222  V(i,i) = PIVOT;
223  }
224 
225  h[i] = 1.0 / sqrt(V(i,i));
226  }
227 
228  // normalisation
229 
230  for (size_t row = 0; row != N; ++row) {
231  for (size_t col = 0; col != row; ++col) {
232  V(row,col) *= h[row] * h[col];
233  V(col,row) = V(row,col);
234  }
235  }
236 
237  for (size_t i = 0; i != N; ++i) {
238  V(i,i) = 1.0 + lambda;
239  }
240 
241  // solve A x = b
242 
243  for (size_t col = 0; col != N; ++col) {
244  x[col] = h[col] * get(current.result.gradient, parameters[col]);
245  }
246 
247  try {
248  V.solve(x);
249  }
250  catch (const exception& error) {
251 
252  ERROR("JGandalf: " << error.what() << endl << V << endl);
253 
254  break;
255  }
256 
257  // update value and error
258 
259  for (size_t row = 0; row != N; ++row) {
260 
261  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
262 
263  get(value, parameters[row]) -= h[row] * x[row];
264  get(error, parameters[row]) = h[row];
265 
266  DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
267  }
268  }
269 
270  // abnormal end
271 
272  const result_type result = previous.result;
273 
274  value = previous.value; // restore value
275  error = previous.error; // restore error
276 
277  reset();
278 
279  update(fit, __begin, __end, args...);
280 
281  try {
282  V.invert();
283  }
284  catch (const exception& error) {
285  V.reset();
286  }
287 
288  return result;
289  }
double lambda
control parameter
Definition: JGandalf.hh:302
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:293
JMATH::JVectorND x
Definition: JGandalf.hh:445
static double PIVOT
minimal value diagonal element of Hesse matrix
Definition: JGandalf.hh:298
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:305
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:456
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:300
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:443
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:294
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:297
void reset()
Reset current parameters.
Definition: JGandalf.hh:311
#define ERROR(A)
Definition: JMessage.hh:66
void update(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method to update current parameters.
Definition: JGandalf.hh:331
struct JFIT::JGandalf::@10 current
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:296
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:75
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JMatrixNS.hh:308
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:292
int numberOfIterations
number of iterations
Definition: JGandalf.hh:301
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:295
struct JFIT::JGandalf::@11 previous
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:444
const JPDFs_t& JFIT::JRegressorStorage< JLine3Z, JGandalf >::getPDF ( ) const
inlineinherited

Get PDFs.

Returns
PDFs

Definition at line 172 of file JRegressorHelper.hh.

173  {
174  return _pdf;
175  }
const JNPEs_t& JFIT::JRegressorStorage< JLine3Z, JGandalf >::getNPE ( ) const
inlineinherited

Get NPEs.

Returns
NPEs

Definition at line 183 of file JRegressorHelper.hh.

184  {
185  return _npe;
186  }
void JFIT::JRegressorStorage< JLine3Z, JGandalf >::transform ( const transformer_type transformer)
inlineinherited

Transform PDFs and NPEs.

Parameters
transformertransformer

Definition at line 194 of file JRegressorHelper.hh.

195  {
196  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
197  _pdf[i].transform(transformer);
198  _npe[i].transform(transformer);
199  }
200  }
static const int NUMBER_OF_PDFS
Number of PDFs.
void JFIT::JRegressorStorage< JLine3Z, JGandalf >::setRmax ( const double  Rmax)
inlineinherited

Set maximal road width of PDF.

Parameters
Rmaxroad width [m]

Definition at line 208 of file JRegressorHelper.hh.

209  {
211 
212  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
213  if (!_pdf[i].empty()) {
214  _pdf[i].compress(range);
215  }
216  }
217  }
static const int NUMBER_OF_PDFS
Number of PDFs.
Range of values.
Definition: JRange.hh:38
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $

Member Data Documentation

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

Time window with respect to Cherenkov hypothesis [ns].

Default values.

Definition at line 439 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 440 of file JLine3ZRegressor.hh.

double JFIT::JRegressor< JLine3Z, JGandalf >::Rmin_m = 0.1
static

Minimal distance of [m].

Definition at line 441 of file JLine3ZRegressor.hh.

PDF.

Definition at line 443 of file JLine3ZRegressor.hh.

PDF.

Definition at line 444 of file JLine3ZRegressor.hh.

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

Energy of muon at vertex [GeV].

Definition at line 446 of file JLine3ZRegressor.hh.

M-Estimator function.

Definition at line 448 of file JLine3ZRegressor.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 292 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 293 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 294 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 295 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 296 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 297 of file JGandalf.hh.

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

minimal value diagonal element of Hesse matrix

minimal value diagonal element of matrix

Definition at line 298 of file JGandalf.hh.

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

fit parameters

Definition at line 300 of file JGandalf.hh.

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

number of iterations

Definition at line 301 of file JGandalf.hh.

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

control parameter

Definition at line 302 of file JGandalf.hh.

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

value

Definition at line 303 of file JGandalf.hh.

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

error

Definition at line 304 of file JGandalf.hh.

Hesse matrix.

Definition at line 305 of file JGandalf.hh.

result_type JFIT::JGandalf< JLine3Z >::result
inherited

Definition at line 448 of file JGandalf.hh.

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

debug level (default is off).

Definition at line 45 of file JMessage.hh.

const int JFIT::JRegressorHelper< JLine3Z, JGandalf >::NUMBER_OF_PDFS = 6
staticinherited

Number of PDFs.

Definition at line 70 of file JRegressorHelper.hh.

const JPDFType_t JFIT::JRegressorHelper< JLine3Z, JGandalf >::pdf_t
staticinherited

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