Jpp - 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::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
 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 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...
 

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...
 
void setRmax (const double Rmax)
 Set 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 two data sets. 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 V
 

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 const int NUMBER_OF_PDFS = 6
 
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 111 of file JLine3ZRegressor.hh.

Member Typedef Documentation

Definition at line 116 of file JLine3ZRegressor.hh.

Definition at line 119 of file JLine3ZRegressor.hh.

time dependent PDF

Definition at line 120 of file JLine3ZRegressor.hh.

Definition at line 124 of file JLine3ZRegressor.hh.

time integrated PDF

Definition at line 125 of file JLine3ZRegressor.hh.

Definition at line 78 of file JRegressor.hh.

Definition at line 79 of file JRegressor.hh.

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

Definition at line 80 of file JRegressor.hh.

Data type of fit parameter.

Definition at line 63 of file JGandalf.hh.

Constructor & Destructor Documentation

Default constructor.

Definition at line 130 of file JLine3ZRegressor.hh.

131  {};
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 147 of file JLine3ZRegressor.hh.

150  :
151  E_GeV(0.0),
153  {
154  using namespace std;
155  using namespace JPP;
156 
157  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
158 
159  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
160 
161  try {
162 
163  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
164 
165  NOTICE("loading PDF from file " << file_name << "... " << flush);
166 
167  pdf[i].load(file_name.c_str());
168 
169  NOTICE("OK" << endl);
170 
171  pdf[i].setExceptionHandler(supervisor);
172  }
173  catch(const JException& error) {
174  FATAL(error.what() << endl);
175  }
176  }
177 
178  // Add PDFs
179 
180  for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
181 
182  pdf[ i ].add(pdf[i-1]);
183 
184  JPDF_t buffer;
185 
186  pdf[i-1].swap(buffer);
187 
188  if (TTS > 0.0) {
189  pdf[i].blur(TTS, numberOfPoints, epsilon);
190  } else if (TTS < 0.0) {
191  ERROR("Illegal value of TTS [ns]: " << TTS << endl);
192  }
193 
194  npe[ i ] = JNPE_t(pdf[i]);
195  }
196  }
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:110
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
void add(const JMultiFunction_t &input)
Add function.
double E_GeV
Energy of muon at vertex [GeV].
#define NOTICE(A)
Definition: JMessage.hh:64
Normal M-estimator.
Definition: JMEstimator.hh:66
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
time integrated PDF
#define FATAL(A)
Definition: JMessage.hh:67
int numberOfPoints
Definition: JResultPDF.cc:22
std::string getFilename(const std::string &file_name)
Get file name part, i.e. part after last JEEP::PATHNAME_SEPARATOR if any.
Definition: JeepToolkit.hh:88
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
time dependent PDF
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 217 of file JLine3ZRegressor.hh.

218  {
219  using namespace JPP;
220 
221  JPosition3D D(hit.getX(), hit.getY(), hit.getZ());
222  JDirection3D U(hit.getDX(), hit.getDY(), hit.getDZ());
223 
224  D.sub(track.getPosition());
225 
226  const double z = D.getDot(track.getDirection());
227  const double x = D.getX() - z * track.getDX();
228  const double y = D.getY() - z * track.getDY();
229  const double R2 = D.getLengthSquared() - z*z;
230  const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
231 
232  const double t1 = track.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
233 
234  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
235 
236  const double theta = U.getTheta();
237  const double phi = fabs(U.getPhi());
238 
239  const double E = gWater.getE(E_GeV, z);
240  const double dt = T_ns.constrain(hit.getT() - t1);
241 
242  JPDF_t::result_type H0 = getH0(hit.getR(), dt);
243  JPDF_t::result_type H1 = getH1(E, R, theta, phi, dt);
244 
245  if (H1.V >= Vmax_npe) {
246  H1 *= Vmax_npe / H1.V;
247  }
248 
249  H1 += H0; // signal + background
250 
252 
253  result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio
254 
255  const double wc = 1.0 - getTanThetaC() * z / R; // d(ct1)/d(z)
256 
257  result.gradient = JLine3Z(JLine1Z(JPosition3D(-getTanThetaC() * D.getX() / R, // d(ct1)/d(x)
258  -getTanThetaC() * D.getY() / R, // d(ct1)/d(y)
259  0.0),
260  getSpeedOfLight()), // d(ct1)/d(t)
261  JVersor3Z(wc * (D.getX() - D.getZ()*track.getDX()/track.getDZ()), // d(ct1)/d(dx)
262  wc * (D.getY() - D.getZ()*track.getDY()/track.getDZ()))); // d(ct1)/d(dy)
263 
264  result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
265  H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1)
266 
267  return result;
268  }
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
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:33
do rm f tmp H1
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:323
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:80
return result
Definition: JPolint.hh:727
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
virtual double getE(const double E, const double dx) const override
Get energy of muon after specified distance.
Definition: JGeane.hh:255
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
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:35
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 279 of file JLine3ZRegressor.hh.

280  {
281  using namespace JGEOMETRY3D;
282  using namespace JPHYSICS;
283 
284  JPosition3D D(pmt.getPosition());
285  JDirection3D U(pmt.getDirection());
286 
287  D.sub(track.getPosition());
288 
289  const double z = D.getDot(track.getDirection());
290  const double x = D.getX() - z * track.getDX();
291  const double y = D.getY() - z * track.getDY();
292  const double R2 = D.getLengthSquared() - z*z;
293  const double R = (R2 > Rmin_m*Rmin_m ? sqrt(R2) : Rmin_m);
294 
295  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
296 
297  const double theta = U.getTheta();
298  const double phi = fabs(U.getPhi());
299 
300  const double E = gWater.getE(E_GeV, z);
301 
302  JNPE_t::result_type H0 = getH0(pmt.getR());
303  JNPE_t::result_type H1 = getH1(E, R, theta, phi);
304 
305  if (H1.f >= Vmax_npe) {
306  H1 *= Vmax_npe / H1.f;
307  }
308 
309  H1 += H0;
310 
311  const bool hit = pmt.getN() != 0;
312  const double u = H1.getChi2(hit);
313 
315 
316  result.chi2 = estimator->getRho(u);
317 
318  result.gradient = JLine3Z(JLine1Z(JPosition3D(-D.getX() / R, // d(R)/d(x)
319  -D.getY() / R, // d(R)/d(y)
320  0.0),
321  0.0),
322  JVersor3Z(-z * D.getX() / R, // d(R)/d(dx)
323  -z * D.getY() / R)); // d(R)/d(dy)
324 
325  result.gradient.mul(estimator->getPsi(u));
326  result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(R)
327 
328  return result;
329  }
do echo Generating $dir eval D
Definition: JDrawLED.sh:50
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:33
do rm f tmp H1
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:323
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:80
return result
Definition: JPolint.hh:727
double E_GeV
Energy of muon at vertex [GeV].
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
virtual double getE(const double E, const double dx) const override
Get energy of muon after specified distance.
Definition: JGeane.hh:255
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:739
Data structure for normalised vector in positive z-direction.
Definition: JVersor3Z.hh:39
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:35
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 339 of file JLine3ZRegressor.hh.

341  {
342  return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
343  }
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 356 of file JLine3ZRegressor.hh.

361  {
362  using namespace std;
363  using namespace JPP;
364 
365  JPDF_t::result_type h1 = zero;
366 
367  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
368 
369  if (!pdf[i].empty() && R <= pdf[i].getXmax()) {
370 
371  JPDF_t::result_type y1 = pdf[i](max(R, pdf[i].getXmin()), theta, phi, t1);
372 
373  // safety measures
374 
375  if (y1.f <= 0.0) {
376  y1.f = 0.0;
377  y1.fp = 0.0;
378  }
379 
380  if (y1.v <= 0.0) {
381  y1.v = 0.0;
382  }
383 
384  // energy dependence
385 
386  if (is_deltarays(pdf_t[i])) {
387  y1 *= getDeltaRaysFromMuon(E);
388  } else if (is_bremsstrahlung(pdf_t[i])) {
389  y1 *= E;
390  }
391 
392  h1 += y1;
393  }
394  }
395 
396  return h1;
397  }
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
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:163
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:177
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:75
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:35
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 406 of file JLine3ZRegressor.hh.

407  {
408  return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
409  }
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 421 of file JLine3ZRegressor.hh.

425  {
426  using namespace std;
427  using namespace JPP;
428 
429  JNPE_t::result_type h1 = JMATH::zero;
430 
431  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
432 
433  if (!npe[i].empty() && R <= npe[i].getXmax()) {
434 
435  try {
436 
437  JNPE_t::result_type y1 = npe[i](max(R, npe[i].getXmin()), theta, phi);
438 
439  // safety measures
440 
441  if (y1.f > 0.0) {
442 
443  // energy dependence
444 
445  if (is_deltarays(pdf_t[i])) {
446  y1 *= getDeltaRaysFromMuon(E);
447  } else if (is_bremsstrahlung(pdf_t[i])) {
448  y1 *= E;
449  }
450 
451  h1 += y1;
452  }
453  }
454  catch(JException& error) {
455  ERROR(error << endl);
456  }
457  }
458  }
459 
460  return h1;
461  }
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
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:163
then echo Variable JPP_DIR undefined exit fi source $JPP_DIR setenv sh $JPP_DIR set_variable NORTH set_variable EAST set_variable SOUTH set_variable WEST set_variable WORKDIR tmp set_variable R set_variable CT set_variable YMAX set_variable YMIN if do_usage *then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:35
bool is_deltarays(const int pdf)
Test if given PDF type corresponds to Cherenkov light from delta-rays.
Definition: JPDFTypes.hh:177
double getDeltaRaysFromMuon(const double E)
Equivalent EM-shower energy due to delta-rays per unit muon track length.
Definition: JPDFToolkit.hh:75
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:35
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 469 of file JLine3ZRegressor.hh.

470  {
471  double xmax = 0.0;
472 
473  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
474  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
475  xmax = pdf[i].getXmax();
476  }
477  }
478 
479  return xmax;
480  }
void JFIT::JRegressor< JLine3Z, JGandalf >::setRmax ( const double  Rmax)
inline

Set maximal road width of PDF.

Parameters
Rmaxroad width [m]

Definition at line 488 of file JLine3ZRegressor.hh.

489  {
490  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
491  if (!pdf[i].empty()) {
493  }
494  }
495  }
void compress(const JRange< typename function_type::abscissa_type > &range)
Compresses PDF to given abscissa range.
Definition: JPDFTable.hh:195
Range of values.
Definition: JRange.hh:38
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 92 of file JRegressor.hh.

93  {
94  static_cast<minimiser_type&>(*this).value = value;
95 
96  return static_cast<minimiser_type&>(*this)(static_cast<regressor_type&>(*this), __begin, __end);
97  }
JModel_t value
Definition: JGandalf.hh:257
JRegressor< JLine3Z, JGandalf > regressor_type
Definition: JRegressor.hh:79
result_type JFIT::JGandalf< JLine3Z >::operator() ( const JFunction_t &  fit,
T  __begin,
T  __end,
Args...  args 
)
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
__beginbegin of data
__endend of data
argsoptional data
Returns
chi2

Definition at line 128 of file JGandalf.hh.

129  {
130  using namespace std;
131  using namespace JPP;
132 
133  const int N = parameters.size();
134 
135  V.resize(N);
136  h.resize(N);
137 
138  previous = value;
139 
140  successor.gradient = value;
141  successor.gradient = zero;
142 
143  error = value;
144  error = zero;
145 
146  lambda = LAMBDA_MIN;
147 
148  result_type precessor = result_type(numeric_limits<double>::max(), value);
149 
151 
152  DEBUG("step: " << numberOfIterations << endl);
153 
154  reset();
155 
156  __evaluate__(fit, __begin, __end, args...);
157 
158  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
159  DEBUG("chi2: " << FIXED(12,5) << successor.chi2 << endl);
160 
161  if (successor.chi2 < precessor.chi2) {
162 
163  if (numberOfIterations != 0) {
164 
165  if (fabs(precessor.chi2 - successor.chi2) < EPSILON*fabs(precessor.chi2)) {
166  return successor;
167  }
168 
169  if (lambda > LAMBDA_MIN) {
170  lambda /= LAMBDA_DOWN;
171  }
172  }
173 
174  precessor = successor;
175  previous = value;
176 
177  } else {
178 
179  value = previous;
180  lambda *= LAMBDA_UP;
181 
182  if (lambda > LAMBDA_MAX) {
183  return precessor; // no improvement found
184  }
185 
186  reset();
187 
188  __evaluate__(fit, __begin, __end, args...);
189  }
190 
191  DEBUG("Hesse matrix:" << endl);
192  DEBUG(V << endl);
193 
194 
195  // force definite positiveness
196 
197  for (int i = 0; i != N; ++i) {
198 
199  if (V(i,i) < PIVOT) {
200  V(i,i) = PIVOT;
201  }
202 
203  h[i] = 1.0 / sqrt(V(i,i));
204  }
205 
206 
207  // normalisation
208 
209  for (int i = 0; i != N; ++i) {
210  for (int j = 0; j != i; ++j) {
211  V(j,i) *= h[i] * h[j];
212  V(i,j) = V(j,i);
213  }
214  }
215 
216  for (int i = 0; i != N; ++i) {
217  V(i,i) = 1.0 + lambda;
218  }
219 
220 
221  try {
222  V.invert();
223  }
224  catch (const JException& error) {
225  ERROR("JGandalf: " << error.what() << endl);
226  return precessor;
227  }
228 
229 
230  for (int i = 0; i != N; ++i) {
231 
232  DEBUG("u[" << noshowpos << i << "] = " << showpos << FIXED(15,5) << alias(value, parameters[i]));
233 
234  for (int j = 0; j != N; ++j) {
235  alias(value, parameters[i]) -= V(i,j) * alias(successor.gradient, parameters[j]) * h[i] * h[j];
236  }
237 
238  DEBUG(' ' << FIXED(15,5) << alias(value, parameters[i]) << noshowpos << endl);
239 
240  alias(error, parameters[i]) = h[i];
241  }
242  }
243 
244  return precessor;
245  }
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:249
static double PIVOT
minimal value diagonal element of matrix
Definition: JGandalf.hh:254
static double alias(const JLine3Z &model, double JLine3Z::*parameter)
Read/write access to parameter value by data member.
Definition: JGandalf.hh:327
JMATH::JMatrixNS V
Definition: JGandalf.hh:261
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
std::vector< parameter_type > parameters
Definition: JGandalf.hh:259
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:362
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:250
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:253
#define ERROR(A)
Definition: JMessage.hh:66
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:252
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:80
void __evaluate__(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method for evaluation of fit.
Definition: JGandalf.hh:288
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:248
int j
Definition: JPolint.hh:666
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:251
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:35
std::vector< double > h
Definition: JGandalf.hh:400

Member Data Documentation

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

Time window with respect to Cherenkov hypothesis [ns].

Default values.

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

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

Minimal distance of [m].

Definition at line 500 of file JLine3ZRegressor.hh.

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

Definition at line 502 of file JLine3ZRegressor.hh.

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

PDF.

Definition at line 506 of file JLine3ZRegressor.hh.

PDF.

Definition at line 507 of file JLine3ZRegressor.hh.

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

Energy of muon at vertex [GeV].

Definition at line 508 of file JLine3ZRegressor.hh.

M-Estimator function.

Definition at line 510 of file JLine3ZRegressor.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 248 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 249 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 250 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 251 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 252 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 253 of file JGandalf.hh.

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

minimal value diagonal element of matrix

Definition at line 254 of file JGandalf.hh.

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

Definition at line 256 of file JGandalf.hh.

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

Definition at line 257 of file JGandalf.hh.

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

Definition at line 258 of file JGandalf.hh.

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

Definition at line 259 of file JGandalf.hh.

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

Definition at line 260 of file JGandalf.hh.

Definition at line 261 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.


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