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

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:347
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:346
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  // normal end
204 
205  const result_type result = current.result;
206 
207  lambda = LAMBDA_MIN;
208 
209  reset();
210 
211  update(fit, __begin, __end, args...);
212 
213  try {
214  V.invert();
215  }
216  catch (const exception& error) {
217  V.reset();
218  }
219 
220  for (size_t i = 0; i != N; ++i) {
221  get(error, parameters[i]) = sqrt(V(i,i));
222  }
223 
224  return result;
225  }
226 
227  if (lambda > LAMBDA_MIN) {
228  lambda /= LAMBDA_DOWN;
229  }
230  }
231 
232  // store current values
233 
234  previous.value = value;
235  previous.result = current.result;
236 
237  } else {
238 
239  value = previous.value; // restore value
240 
241  lambda *= LAMBDA_UP;
242 
243  if (lambda > LAMBDA_MAX) {
244  break;
245  }
246 
247  reset();
248 
249  update(fit, __begin, __end, args...);
250  }
251 
252  DEBUG("Hesse matrix:" << endl << V << endl);
253 
254  // force definite positiveness
255 
256  for (size_t i = 0; i != N; ++i) {
257 
258  if (V(i,i) < PIVOT) {
259  V(i,i) = PIVOT;
260  }
261 
262  h[i] = 1.0 / sqrt(V(i,i));
263  }
264 
265  // normalisation
266 
267  for (size_t row = 0; row != N; ++row) {
268  for (size_t col = 0; col != row; ++col) {
269  V(row,col) *= h[row] * h[col];
270  V(col,row) = V(row,col);
271  }
272  }
273 
274  for (size_t i = 0; i != N; ++i) {
275  V(i,i) = 1.0 + lambda;
276  }
277 
278  // solve A x = b
279 
280  for (size_t col = 0; col != N; ++col) {
281  x[col] = h[col] * get(current.result.gradient, parameters[col]);
282  }
283 
284  try {
285  V.solve(x);
286  }
287  catch (const exception& error) {
288 
289  ERROR("JGandalf: " << error.what() << endl << V << endl);
290 
291  break;
292  }
293 
294  // update value
295 
296  for (size_t row = 0; row != N; ++row) {
297 
298  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
299 
300  get(value, parameters[row]) -= h[row] * x[row];
301 
302  DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
303  }
304 
305  model(value);
306  }
307 
308  // abnormal end
309 
310  const result_type result = previous.result;
311 
312  value = previous.value; // restore value
313 
314  lambda = LAMBDA_MIN;
315 
316  reset();
317 
318  update(fit, __begin, __end, args...);
319 
320  try {
321  V.invert();
322  }
323  catch (const exception& error) {
324  V.reset();
325  }
326 
327  for (size_t i = 0; i != N; ++i) {
328  get(error, parameters[i]) = sqrt(V(i,i));
329  }
330 
331  return result;
332  }
double lambda
control parameter
Definition: JGandalf.hh:345
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:336
JShower3EZ value
value
Definition: JGandalf.hh:346
static double PIVOT
minimal value diagonal element of Hesse matrix
Definition: JGandalf.hh:341
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:348
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:343
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:337
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:340
void reset()
Reset current parameters.
Definition: JGandalf.hh:354
#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:374
struct JFIT::JGandalf::@10 current
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:339
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:335
int numberOfIterations
number of iterations
Definition: JGandalf.hh:344
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:338
JShower3EZ error
error
Definition: JGandalf.hh:347
struct JFIT::JGandalf::@11 previous
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:487

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 335 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 336 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 337 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 338 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 339 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 340 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 341 of file JGandalf.hh.

fit parameters

Definition at line 343 of file JGandalf.hh.

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

number of iterations

Definition at line 344 of file JGandalf.hh.

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

control parameter

Definition at line 345 of file JGandalf.hh.

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

value

Definition at line 346 of file JGandalf.hh.

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

error

Definition at line 347 of file JGandalf.hh.

Hesse matrix.

Definition at line 348 of file JGandalf.hh.

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

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