Jpp  18.6.0-rc.1
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< JShower3EZ, JGandalf > Struct Template Reference

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

#include <JShower3EZRegressor.hh>

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

Public Types

typedef JTOOLS::JSplineFunction1S_t JFunction1D_t
 
typedef JTOOLS::JMapList
< JTOOLS::JPolint0FunctionalMap,
JTOOLS::JMapList
< JTOOLS::JPolint0FunctionalMap,
JTOOLS::JMapList
< JTOOLS::JPolint0FunctionalGridMap,
JTOOLS::JMapList
< JTOOLS::JPolint0FunctionalGridMap > > > > 
JPDFMaplist_t
 
typedef JPHYSICS::JPDFTable
< JFunction1D_t, JPDFMaplist_t
JPDF_t
 
typedef JTOOLS::JMapList
< JTOOLS::JPolint1FunctionalMap,
JTOOLS::JMapList
< JTOOLS::JPolint1FunctionalMapH,
JTOOLS::JMapList
< JTOOLS::JPolint1FunctionalGridMap,
JTOOLS::JMapList
< JTOOLS::JPolint1FunctionalGridMap > > > > 
JNPEMaplist_t
 
typedef JPHYSICS::JNPETable
< double, double,
JNPEMaplist_t
JNPE_t
 
typedef JGandalf< JShower3EZminimiser_type
 
typedef JRegressor< JShower3EZ,
JGandalf
regressor_type
 
typedef minimiser_type::result_type result_type
 
typedef JFIT_LOCAL::JTypedef_t
< JShower3EZ >::parameter_type 
parameter_type
 Data type of fit parameter. More...
 

Public Member Functions

 JRegressor (const std::string &fileDescriptor)
 Parameterized constructor. More...
 
result_type operator() (const JShower3EZ &shower, const JPMTW0 &pmt) const
 Fit function. More...
 
JNPE_t::result_type getH0 (const double R_Hz) const
 Get background hypothesis value for time integrated PDF. More...
 
JNPE_t::result_type getH1 (const double D, const double cd, const double theta, const double phi, const double E) const
 Get signal hypothesis value for time integrated PDF. More...
 
result_type operator() (const JShower3EZ &value, T __begin, T __end)
 Global fit. More...
 
result_type operator() (const JFunction_t &fit, T __begin, T __end, Args...args)
 Multi-dimensional fit of multiple data sets. More...
 

Public Attributes

JNPE_t npe [NUMBER_OF_PDFS]
 PDF. More...
 
JLANG::JSharedPointer
< JMEstimator
estimator
 M-Estimator function. More...
 
std::vector< parameter_typeparameters
 fit parameters More...
 
int numberOfIterations
 number of iterations More...
 
double lambda
 control parameter More...
 
JShower3EZ value
 value More...
 
JShower3EZ 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 JFIT::JPDFType_t pdf_t [NUMBER_OF_PDFS]
 
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< JShower3EZ, JGandalf >

Regressor function object for JShower3EZ fit using JGandalf minimiser.

Definition at line 276 of file JShower3EZRegressor.hh.

Member Typedef Documentation

Definition at line 281 of file JShower3EZRegressor.hh.

Definition at line 285 of file JShower3EZRegressor.hh.

Definition at line 286 of file JShower3EZRegressor.hh.

Definition at line 291 of file JShower3EZRegressor.hh.

Definition at line 292 of file JShower3EZRegressor.hh.

Definition at line 80 of file JRegressor.hh.

Definition at line 81 of file JRegressor.hh.

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

Definition at line 82 of file JRegressor.hh.

Data type of fit parameter.

Definition at line 95 of file JGandalf.hh.

Constructor & Destructor Documentation

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

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

Parameters
fileDescriptorPDF file descriptor

Definition at line 303 of file JShower3EZRegressor.hh.

303  :
304  estimator(new JMEstimatorNull())
305  {
306  using namespace std;
307  using namespace JPP;
308 
309  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
310 
311  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
312 
313  try {
314 
315  JPDF_t pdf;
316 
317  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
318 
319  NOTICE("loading PDF from file " << file_name << "... " << flush);
320 
321  pdf.load(file_name.c_str());
322 
323  NOTICE("OK" << endl);
324 
325  pdf.setExceptionHandler(supervisor);
326 
327  npe[i] = JNPE_t(pdf);
328  }
329  catch(const JException& error) {
330  FATAL(error.what() << endl);
331  }
332  }
333 
334  // Add PDFs
335  for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
336 
337  npe[ i ].add(npe[i-1]);
338 
339  JNPE_t buffer;
340 
341  npe[i-1].swap(buffer);
342  }
343  }
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
static const JFIT::JPDFType_t pdf_t[NUMBER_OF_PDFS]
void add(const JNPETable &input)
Add NPE table.
Definition: JNPETable.hh:124
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
#define NOTICE(A)
Definition: JMessage.hh:64
#define FATAL(A)
Definition: JMessage.hh:67
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
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
JShower3EZ error
error
Definition: JGandalf.hh:334

Member Function Documentation

result_type JFIT::JRegressor< JShower3EZ, JGandalf >::operator() ( const JShower3EZ shower,
const JPMTW0 pmt 
) const
inline

Fit function.

This method is used to determine the chi2 of given PMT with respect to shower hypothesis.

Parameters
showershower
pmtpmt
Returns
chi2

Definition at line 353 of file JShower3EZRegressor.hh.

354  {
355  using namespace JPP;
356  using namespace std;
357 
358  JPosition3D D(pmt.getPosition());
359  JDirection3D U(pmt.getDirection());
360 
361  D.sub(shower.getPosition());
362 
363  const double x = D.getX();
364  const double y = D.getY();
365  const double d = D.getLength();
366  const double cd = D.getDot(shower.getDirection())/d; // cosine angle between shower direction and PMT position
367 
368  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
369 
370  const double theta = getPMTAngle(U.getTheta());
371  const double phi = getPMTAngle(U.getPhi());
372 
373  JNPE_t::result_type H0 = getH0(pmt.getR()); // background hypothesis value for time integrated PDF.
374  JNPE_t::result_type H1 = getH1(d, cd, theta, phi, shower.getE()); // signal hypothesis value for time integrated PDF.
375 
376  if (get_value(H1) >= Vmax_npe) {
377  H1 *= Vmax_npe / get_value(H1);
378  }
379 
380  double signal_npe = get_value(H1);
381 
382  H1 += H0; // now H1 is signal + background
383 
384  double expectation = get_value(H1);
385 
386  const bool hit = pmt.getN() != 0;
387  const double u = H1.getChi2(hit);
388 
390 
391  result.chi2 = estimator->getRho(u);
392 
393  double energy_gradient = signal_npe/shower.getE(); // dP/dE
394  if(hit) energy_gradient *= -exp(-expectation)/(1-exp(-expectation)); //dchi2/d(H1), if !hit is 1
395 
396  result.gradient = JShower3EZ(JPoint4D(JVector3D(0, // d(cos_th0)/d(x)
397  0, // d(cos_th0)/d(y)
398  0), // d(cos_th0)/d(z)
399  0.0), // d(cos_th0)/d(t)
400  JVersor3Z(x/d, // d(cos_th0)/d(dx)
401  y/d), // d(cos_th0)/d(dy)
402  energy_gradient); // d(chi2)/d(E)
403 
404  result.gradient.mul(estimator->getPsi(u));
405  static_cast<JShower3Z&>(result.gradient).mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(cos_th0)
406 
407  return result;
408  }
double getPMTAngle(const double angle)
Constrain PMT angle to [0,pi].
JNPE_t::result_type getH1(const double D, const double cd, const double theta, const double phi, const double E) const
Get signal hypothesis value for time integrated PDF.
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonStart.sh:47
static double Vmax_npe
Maximal integral of PDF [npe].
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
double u[N+1]
Definition: JPolint.hh:865
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
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
JNPE_t::result_type getH0(const double R_Hz) const
Get background hypothesis value for time integrated PDF.
JNPE_t::result_type JFIT::JRegressor< JShower3EZ, 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 416 of file JShower3EZRegressor.hh.

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

Get signal hypothesis value for time integrated PDF.

Parameters
DPMT distance from shower [m]
cdcosine angle between shower direction and PMT position
thetaPMT zenith angle [deg]
phiPMT azimuth angle [deg]
Eshower energy [GeV]
Returns
hypothesis value

Definition at line 431 of file JShower3EZRegressor.hh.

436  {
437  JNPE_t::result_type h1 = JMATH::zero;
438 
439  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
440 
441  if (!npe[i].empty() && D <= npe[i].getXmax()) {
442 
443  try {
444 
445  JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cd, theta, phi);
446 
447  if (get_value(y1) > 0.0) {
448  h1 += y1;
449  }
450 
451  }
452  catch(JLANG::JException& error) {
453  ERROR(error << std::endl);
454  }
455  }
456  }
457 
458  return h1;
459  }
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
JShower3EZ error
error
Definition: JGandalf.hh:334
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
result_type JFIT::JAbstractRegressor< JShower3EZ , JGandalf >::operator() ( const JShower3EZ 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:333
JRegressor< JShower3EZ, JGandalf > regressor_type
Definition: JRegressor.hh:81
result_type JFIT::JGandalf< JShower3EZ >::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 160 of file JGandalf.hh.

161  {
162  using namespace std;
163  using namespace JPP;
164 
165  // note that all model values should be assigned to the start value of the model before use
166  // because the actual list of model parameters can vary from fit to fit
167  // (e.g. if model consists of a container).
168 
169  const size_t N = parameters.size();
170 
171  V.resize(N);
172  h.resize(N);
173  x.resize(N);
174 
175  previous.result.chi2 = numeric_limits<double>::max();
176 
177  current.result.chi2 = numeric_limits<double>::max();
178  current.result.gradient = value;
179  current.result.gradient = zero;
180 
181  error = value;
182  error = zero;
183 
184  lambda = LAMBDA_MIN;
185 
187 
188  DEBUG("step: " << numberOfIterations << endl);
189 
190  reset();
191 
192  update(fit, __begin, __end, args...);
193 
194  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
195  DEBUG("chi2: " << FIXED(12,5) << current.result.chi2 << endl);
196 
197  if (current.result.chi2 < previous.result.chi2) {
198 
199  if (numberOfIterations != 0) {
200 
201  if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
202 
203  const result_type result = current.result;
204 
205  reset();
206 
207  update(fit, __begin, __end, args...);
208 
209  try {
210  V.invert();
211  }
212  catch (const exception& error) {
213  V.reset();
214  }
215 
216  return result; // normal end
217  }
218 
219  if (lambda > LAMBDA_MIN) {
220  lambda /= LAMBDA_DOWN;
221  }
222  }
223  // store current values
224 
225  previous.value = value;
226  previous.error = error;
227  previous.result = current.result;
228 
229  } else {
230 
231  value = previous.value; // restore value
232 
233  lambda *= LAMBDA_UP;
234 
235  if (lambda > LAMBDA_MAX) {
236  break;
237  }
238 
239  reset();
240 
241  update(fit, __begin, __end, args...);
242  }
243 
244  DEBUG("Hesse matrix:" << endl << V << endl);
245 
246  // force definite positiveness
247 
248  for (size_t i = 0; i != N; ++i) {
249 
250  if (V(i,i) < PIVOT) {
251  V(i,i) = PIVOT;
252  }
253 
254  h[i] = 1.0 / sqrt(V(i,i));
255  }
256 
257  // normalisation
258 
259  for (size_t row = 0; row != N; ++row) {
260  for (size_t col = 0; col != row; ++col) {
261  V(row,col) *= h[row] * h[col];
262  V(col,row) = V(row,col);
263  }
264  }
265 
266  for (size_t i = 0; i != N; ++i) {
267  V(i,i) = 1.0 + lambda;
268  }
269 
270  // solve A x = b
271 
272  for (size_t col = 0; col != N; ++col) {
273  x[col] = h[col] * get(current.result.gradient, parameters[col]);
274  }
275 
276  try {
277  V.solve(x);
278  }
279  catch (const exception& error) {
280 
281  ERROR("JGandalf: " << error.what() << endl << V << endl);
282 
283  break;
284  }
285 
286  // update value and error
287 
288  for (size_t row = 0; row != N; ++row) {
289 
290  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
291 
292  get(value, parameters[row]) -= h[row] * x[row];
293  get(error, parameters[row]) = h[row];
294 
295  DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
296  }
297  model(value);
298  }
299 
300  // abnormal end
301 
302  const result_type result = previous.result;
303 
304  value = previous.value; // restore value
305  error = previous.error; // restore error
306 
307  reset();
308 
309  update(fit, __begin, __end, args...);
310 
311  try {
312  V.invert();
313  }
314  catch (const exception& error) {
315  V.reset();
316  }
317 
318  return result;
319  }
double lambda
control parameter
Definition: JGandalf.hh:332
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:323
JShower3EZ value
value
Definition: JGandalf.hh:333
static double PIVOT
minimal value diagonal element of Hesse matrix
Definition: JGandalf.hh:328
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:335
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:459
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:330
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:324
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:327
void reset()
Reset current parameters.
Definition: JGandalf.hh:341
#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:361
struct JFIT::JGandalf::@10 current
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:326
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:322
int numberOfIterations
number of iterations
Definition: JGandalf.hh:331
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:325
JShower3EZ error
error
Definition: JGandalf.hh:334
struct JFIT::JGandalf::@11 previous
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:474

Member Data Documentation

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

Time window with respect to Cherenkov hypothesis [ns].

Definition at line 462 of file JShower3EZRegressor.hh.

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

Maximal integral of PDF [npe].

Definition at line 463 of file JShower3EZRegressor.hh.

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

Definition at line 465 of file JShower3EZRegressor.hh.

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

Definition at line 467 of file JShower3EZRegressor.hh.

PDF.

Definition at line 469 of file JShower3EZRegressor.hh.

M-Estimator function.

Definition at line 471 of file JShower3EZRegressor.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 322 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 323 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 324 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 325 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 326 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 327 of file JGandalf.hh.

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

minimal value diagonal element of Hesse matrix

minimal value diagonal element of matrix

Definition at line 328 of file JGandalf.hh.

fit parameters

Definition at line 330 of file JGandalf.hh.

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

number of iterations

Definition at line 331 of file JGandalf.hh.

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

control parameter

Definition at line 332 of file JGandalf.hh.

JShower3EZ JFIT::JGandalf< JShower3EZ >::value
inherited

value

Definition at line 333 of file JGandalf.hh.

JShower3EZ JFIT::JGandalf< JShower3EZ >::error
inherited

error

Definition at line 334 of file JGandalf.hh.

Hesse matrix.

Definition at line 335 of file JGandalf.hh.

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

Definition at line 478 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: