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< JPoint4D, JGandalf > Struct Template Reference

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

#include <JShowerBrightPointRegressor.hh>

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

Public Types

typedef JTOOLS::JSplineFunction1S_t JFunction1D_t
 
typedef JTOOLS::JMAPLIST
< JTOOLS::JPolint2FunctionalMap,
JTOOLS::JPolint1FunctionalGridMap >
::maplist 
JPDFMapList_t
 
typedef JPHYSICS::JPDFTable
< JFunction1D_t, JPDFMapList_t
JPDF_t
 
typedef JGandalf< JPoint4Dminimiser_type
 
typedef JRegressor< JPoint4D,
JGandalf
regressor_type
 
typedef minimiser_type::result_type result_type
 
typedef JPoint4D::parameter_type parameter_type
 Data type of fit parameter. More...
 

Public Member Functions

 JRegressor (const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
 Parameterized constructor. More...
 
template<class JHit_t >
result_type operator() (const JPoint4D &vx, const JHit_t &hit) 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 D, const double ct, const double t, const double E) const
 Get signal hypothesis value for bright point emission PDF. More...
 
double getRmax () const
 Get maximal road width of PDF. More...
 
result_type operator() (const JPoint4D &value, T __begin, T __end)
 Global fit. More...
 
result_type operator() (const JFunction_t &fit, T __begin, T __end, Args...args)
 Multi-dimensional fit of multiple data sets. More...
 

Public Attributes

JPDF_t pdf [NUMBER_OF_PDFS]
 PDF. More...
 
double E_GeV
 Energy of the shower [GeV]. More...
 
std::vector< parameter_typeparameters
 fit parameters More...
 
int numberOfIterations
 number of iterations More...
 
double lambda
 control parameter More...
 
JPoint4D value
 value More...
 
JPoint4D error
 error More...
 
JMATH::JMatrixNS V
 Hesse matrix. More...
 
result_type result
 

Static Public Attributes

static JTimeRange T_ns
 Time window with respect to Cherenkov hypothesis [ns]. More...
 
static double Vmax_npe = std::numeric_limits<double>::max()
 Maximal integral of PDF [npe]. More...
 
static const int NUMBER_OF_PDFS = 2
 
static const 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 Hesse matrix More...
 
static int debug = 0
 debug level (default is off). More...
 

Detailed Description

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

Regressor function object for JPoint4D fit using JGandalf minimiser.

Definition at line 43 of file JShowerBrightPointRegressor.hh.

Member Typedef Documentation

Definition at line 48 of file JShowerBrightPointRegressor.hh.

Definition at line 50 of file JShowerBrightPointRegressor.hh.

Definition at line 51 of file JShowerBrightPointRegressor.hh.

Definition at line 80 of file JRegressor.hh.

Definition at line 81 of file JRegressor.hh.

typedef minimiser_type::result_type JFIT::JAbstractRegressor< JPoint4D , 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.

Constructor & Destructor Documentation

JFIT::JRegressor< JPoint4D, JGandalf >::JRegressor ( const std::string fileDescriptor,
const double  TTS,
const int  numberOfPoints = 25,
const double  epsilon = 1.0e-10 
)
inline

Parameterized constructor.

The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which will be replaced by the corresponding PDF types.

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 65 of file JShowerBrightPointRegressor.hh.

68  :
69  E_GeV(0.0)
70  {
71  using namespace std;
72  using namespace JPP;
73 
74  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
75 
76  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
77 
78  try {
79 
80  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
81 
82  NOTICE("loading PDF from file " << file_name << "... " << flush);
83 
84  pdf[i].load(file_name.c_str());
85 
86  NOTICE("OK" << endl);
87 
88  pdf[i].setExceptionHandler(supervisor);
89  }
90  catch(const JException& error) {
91  FATAL(error.what() << endl);
92  }
93  }
94 
95  // Add PDFs
96  for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
97 
98  pdf[ i ].add(pdf[i-1]);
99 
100  JPDF_t buffer;
101 
102  pdf[i-1].swap(buffer);
103 
104  if (TTS > 0.0) {
105  pdf[i].blur(TTS, numberOfPoints, epsilon);
106  } else if (TTS < 0.0) {
107  ERROR("Illegal value of TTS [ns]: " << TTS << endl);
108  }
109  }
110  }
void load(const char *file_name)
Load from input file.
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMapList_t > JPDF_t
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:111
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
void add(const JMultiFunction_t &input)
Add function.
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
static const JPDFType_t pdf_t[NUMBER_OF_PDFS]
PDF types.
#define FATAL(A)
Definition: JMessage.hh:67
void setExceptionHandler(const typename function_type::supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
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:128
JPoint4D error
error
Definition: JGandalf.hh:304
const double epsilon
Definition: JQuadrature.cc:21

Member Function Documentation

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

Fit function.

This method is used to determine the chi2 and gradient of given hit with respect a bright point emitting isotropically

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
vxshower vertex
hithit
Returns
chi2 and gradient

Definition at line 130 of file JShowerBrightPointRegressor.hh.

131  {
132  using namespace JPP;
133 
134  JPosition3D D(hit.getPosition());
135  JDirection3D U(hit.getDirection());
136 
137  D.sub(vx.getPosition());
138 
139  double ct = U.getDot(D) / D.getLength();
140 
141  if (ct > +1.0) { ct = +1.0; }
142  if (ct < -1.0) { ct = -1.0; }
143 
144  const double t = vx.getT() + (D.getLength() * getIndexOfRefraction() * getInverseSpeedOfLight());
145 
146  const double dt = T_ns.constrain(hit.getT() - t);
147 
148  JPDF_t::result_type H0 = getH0(hit.getR(), dt); // getH0 = Get background hypothesis value
149  JPDF_t::result_type H1 = getH1(D.getLength(), ct, dt, E_GeV); // getH1 = Get signal hypothesis value
150 
151  if (get_value(H1) >= Vmax_npe) {
152  H1 *= Vmax_npe / get_value(H1);
153  }
154 
155  H1 += H0; // now H1 is signal + background
156 
158 
159  result.chi2 = H1.getChi2() - H0.getChi2(); // Likelihood ratio
160 
161  /*
162  * Here it is evaluated: d(chi2)/d(ct) * d(ct)/d(x0,y0,z0,t0)
163  */
164  result.gradient = JPoint4D(JVector3D(-getIndexOfRefraction() * D.getX() / D.getLength(), // d(ct)/d(x0)
165  -getIndexOfRefraction() * D.getY() / D.getLength(), // d(ct)/d(y0)
166  -getIndexOfRefraction() * D.getZ() / D.getLength()), // d(ct)/d(z0)
167  getSpeedOfLight()); // d(ct)/d(t0)
168 
169  result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
170  H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1)
171 
172  return result;
173 
174  }
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Data structure for vertex fit.
Definition: JPoint4D.hh:22
minimiser_type::result_type result_type
Definition: JRegressor.hh:82
JPDF_t::result_type getH0(const double R_Hz, const double t1) const
Get background hypothesis value for time differentiated PDF.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
static double Vmax_npe
Maximal integral of PDF [npe].
const double getSpeedOfLight()
Get speed of light.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:147
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:998
JPDF_t::result_type getH1(const double D, const double ct, const double t, const double E) const
Get signal hypothesis value for bright point emission PDF.
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
JPDF_t::result_type JFIT::JRegressor< JPoint4D, 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 183 of file JShowerBrightPointRegressor.hh.

185  {
186  using namespace JPP;
187 
188  return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
189  }
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
JPDF_t::result_type JFIT::JRegressor< JPoint4D, JGandalf >::getH1 ( const double  D,
const double  ct,
const double  t,
const double  E 
) const
inline

Get signal hypothesis value for bright point emission PDF.

Parameters
Dhit distance from shower vertex [m]
ctcosine of the HIT angle
tarrival time of the light
Eshower energy [GeV]
Returns
hypothesis value

Definition at line 200 of file JShowerBrightPointRegressor.hh.

204  {
205  using namespace JPP;
206 
207  JPDF_t::result_type h1 = JMATH::zero;
208 
209  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
210 
211  if (!pdf[i].empty() && D <= pdf[i].getXmax()) {
212 
213  try {
214 
215  JPDF_t::result_type y1 = E * pdf[i](std::max(D, pdf[i].getXmin()), ct, t);
216 
217  // safety measures
218 
219  if (y1.f <= 0.0) {
220  y1.f = 0.0;
221  y1.fp = 0.0;
222  }
223 
224  if (y1.v <= 0.0) {
225  y1.v = 0.0;
226  }
227 
228  h1 += y1;
229 
230  }
231  catch(JLANG::JException& error) {
232  ERROR(error << std::endl);
233  }
234  }
235  }
236 
237  return h1;
238  }
General exception.
Definition: JException.hh:24
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
#define ERROR(A)
Definition: JMessage.hh:66
JPoint4D error
error
Definition: JGandalf.hh:304
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double JFIT::JRegressor< JPoint4D, JGandalf >::getRmax ( ) const
inline

Get maximal road width of PDF.

Returns
road width [m]

Definition at line 245 of file JShowerBrightPointRegressor.hh.

246  {
247  using namespace JPP;
248 
249  double xmax = 0.0;
250 
251  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
252 
253  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
254  xmax = pdf[i].getXmax();
255  }
256 
257  }
258 
259  return xmax;
260  }
const double xmax
Definition: JQuadrature.cc:24
result_type JFIT::JAbstractRegressor< JPoint4D , JGandalf >::operator() ( const JPoint4D 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< JPoint4D, JGandalf > regressor_type
Definition: JRegressor.hh:81
result_type JFIT::JGandalf< JPoint4D >::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
JPoint4D value
value
Definition: JGandalf.hh:303
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
JPoint4D error
error
Definition: JGandalf.hh:304
struct JFIT::JGandalf::@11 previous
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:444

Member Data Documentation

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

Time window with respect to Cherenkov hypothesis [ns].

Default values.

Definition at line 262 of file JShowerBrightPointRegressor.hh.

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

Maximal integral of PDF [npe].

Definition at line 263 of file JShowerBrightPointRegressor.hh.

const int JFIT::JRegressor< JPoint4D, JGandalf >::NUMBER_OF_PDFS = 2
static

Definition at line 265 of file JShowerBrightPointRegressor.hh.

const JPDFType_t JFIT::JRegressor< JPoint4D, JGandalf >::pdf_t
static
Initial value:

PDF types.

Definition at line 267 of file JShowerBrightPointRegressor.hh.

PDF.

Definition at line 269 of file JShowerBrightPointRegressor.hh.

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

Energy of the shower [GeV].

Definition at line 271 of file JShowerBrightPointRegressor.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 292 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 293 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 294 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 295 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 296 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 297 of file JGandalf.hh.

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

minimal value diagonal element of Hesse matrix

minimal value diagonal element of matrix

Definition at line 298 of file JGandalf.hh.

fit parameters

Definition at line 300 of file JGandalf.hh.

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

number of iterations

Definition at line 301 of file JGandalf.hh.

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

control parameter

Definition at line 302 of file JGandalf.hh.

JPoint4D JFIT::JGandalf< JPoint4D >::value
inherited

value

Definition at line 303 of file JGandalf.hh.

JPoint4D JFIT::JGandalf< JPoint4D >::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< JPoint4D >::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.


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