Jpp  18.0.0-rc.3
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::JPolint1FunctionalMapH,
JTOOLS::JMapList
< JTOOLS::JPolint1FunctionalMap,
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 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 cosDelta, 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 238 of file JShower3EZRegressor.hh.

Member Typedef Documentation

Definition at line 243 of file JShower3EZRegressor.hh.

Definition at line 247 of file JShower3EZRegressor.hh.

Definition at line 248 of file JShower3EZRegressor.hh.

Definition at line 253 of file JShower3EZRegressor.hh.

Definition at line 254 of file JShower3EZRegressor.hh.

Definition at line 78 of file JRegressor.hh.

Definition at line 79 of file JRegressor.hh.

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

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

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

Parameters
fileDescriptorPDF file descriptor

Definition at line 265 of file JShower3EZRegressor.hh.

265  :
266  estimator(new JMEstimatorNull())
267  {
268  using namespace std;
269  using namespace JPP;
270 
271  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
272 
273  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
274 
275  try {
276 
277  JPDF_t pdf;
278 
279  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
280 
281  NOTICE("loading PDF from file " << file_name << "... " << flush);
282 
283  pdf.load(file_name.c_str());
284 
285  NOTICE("OK" << endl);
286 
287  pdf.setExceptionHandler(supervisor);
288 
289  npe[i] = JNPE_t(pdf);
290  }
291  catch(const JException& error) {
292  FATAL(error.what() << endl);
293  }
294  }
295 
296  // Add PDFs
297  for (int i = 1; i < NUMBER_OF_PDFS; i += 2) {
298 
299  npe[ i ].add(npe[i-1]);
300 
301  JNPE_t buffer;
302 
303  npe[i-1].swap(buffer);
304  }
305  }
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:292

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 315 of file JShower3EZRegressor.hh.

316  {
317  using namespace JPP;
318  using namespace std;
319 
320  JPosition3D D(pmt.getPosition());
321  JDirection3D U(pmt.getDirection());
322 
323  D.sub(shower.getPosition());
324 
325  JVersor3D shower_dir(shower.getDirection());
326 
327  const double z = D.getDot(shower_dir);
328  const double x = D.getX() - z * shower.getDX();
329  const double y = D.getY() - z * shower.getDY();
330  const double R2 = D.getLengthSquared() - z*z;
331  const double R = sqrt(R2);
332  const double d = D.getLength();
333  const double cos_th0 = z/d; // theta0 = angle between shower direction and PMT position
334 
335  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
336 
337  const double theta = U.getTheta();
338  const double phi = fabs(U.getPhi());
339 
340  JNPE_t::result_type H0 = getH0(pmt.getR()); // getH0 = Get background hypothesis value for time integrated PDF.
341  JNPE_t::result_type H1 = getH1(d, cos_th0, theta, phi, shower.getE()); // getH1 = Get signal hypothesis value for time integrated PDF.
342 
343  if (get_value(H1) >= Vmax_npe) {
344  H1 *= Vmax_npe / get_value(H1);
345  }
346 
347  H1 += H0; // now H1 is signal + background
348 
349  const bool hit = pmt.getN() != 0;
350  const double u = H1.getChi2(hit);
351 
353 
354  result.chi2 = estimator->getRho(u);
355 
356  result.gradient = JShower3EZ(JPoint4D(JVector3D(x*z / (d*d*d), // d(cos_th0)/d(x)
357  y*z / (d*d*d), // d(cos_th0)/d(y)
358  -(R*R) / (d*d*d)), // d(cos_th0)/d(z)
359  0.0), // d(cos_th0)/d(t)
360  JVersor3Z(x*z*z / (d*R), // d(cos_th0)/d(dx)
361  y*z*z / (d*R)), // d(cos_th0)/d(dy)
362  0.0); // d(cos_th0)/d(E)
363 
364  result.gradient.mul(estimator->getPsi(u));
365  result.gradient.mul(H1.getDerivativeOfChi2(hit)); // x d(chi2)/d(cos_th0)
366 
367  return result;
368  }
JNPE_t::result_type getH1(const double D, const double cosDelta, const double theta, const double phi, const double E) const
Get signal hypothesis value for time integrated PDF.
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.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:776
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
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 376 of file JShower3EZRegressor.hh.

377  {
378  return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
379  }
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  cosDelta,
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]
cosDeltaangle between shower direction and PMT position
thetaPMT zenith angle [deg]
phiPMT azimuth angle [deg]
Eshower energy [GeV]
Returns
hypothesis value

Definition at line 391 of file JShower3EZRegressor.hh.

396  {
397  JNPE_t::result_type h1 = JMATH::zero;
398 
399  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
400 
401  if (!npe[i].empty() && D <= npe[i].getXmax()) {
402 
403  try {
404 
405  JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi);
406 
407  if (get_value(y1) > 0.0) {
408  h1 += y1;
409  }
410 
411  }
412  catch(JLANG::JException& error) {
413  ERROR(error << std::endl);
414  }
415  }
416  }
417 
418  return h1;
419  }
General exception.
Definition: JException.hh:23
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:36
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:292
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 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
value
Definition: JGandalf.hh:291
JRegressor< JShower3EZ, JGandalf > regressor_type
Definition: JRegressor.hh:79
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

Definition at line 128 of file JGandalf.hh.

129  {
130  using namespace std;
131  using namespace JPP;
132 
133  // note that all model values should be assigned to the start value of the model before use
134  // because the actual list of model parameters can vary from fit to fit
135  // (e.g. if model consists of a container).
136 
137  const size_t N = parameters.size();
138 
139  V.resize(N);
140  h.resize(N);
141 
142  previous.result.chi2 = numeric_limits<double>::max();
143 
144  previous.V.resize(N);
145 
146  current.result.chi2 = numeric_limits<double>::max();
147  current.result.gradient = value;
148  current.result.gradient = zero;
149 
150  current.V.resize(N);
151 
152  error = value;
153  error = zero;
154 
155  lambda = LAMBDA_MIN;
156 
157 
159 
160  DEBUG("step: " << numberOfIterations << endl);
161 
162  reset();
163 
164  update(fit, __begin, __end, args...);
165 
166  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
167  DEBUG("chi2: " << FIXED(12,5) << current.result.chi2 << endl);
168 
169  if (current.result.chi2 < previous.result.chi2) {
170 
171  if (numberOfIterations != 0) {
172 
173  if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
174  return current.result; // normal end
175  }
176 
177  if (lambda > LAMBDA_MIN) {
178  lambda /= LAMBDA_DOWN;
179  }
180 
181  // store current values
182 
183  previous.value = value;
184  previous.error = error;
185  previous.result = current.result;
186 
187  previous.V.swap(V);
188  }
189 
190  } else {
191 
192  value = previous.value; // restore value
193 
194  lambda *= LAMBDA_UP;
195 
196  if (lambda > LAMBDA_MAX) {
197  break;
198  }
199 
200  reset();
201 
202  update(fit, __begin, __end, args...);
203  }
204 
205 
206  V.swap(current.V);
207 
208  DEBUG("Hesse matrix:" << endl << V << endl);
209 
210 
211  // force definite positiveness
212 
213  for (size_t i = 0; i != N; ++i) {
214 
215  if (V(i,i) < PIVOT) {
216  V(i,i) = PIVOT;
217  }
218 
219  h[i] = 1.0 / sqrt(V(i,i));
220  }
221 
222 
223  // normalisation
224 
225  for (size_t row = 0; row != N; ++row) {
226  for (size_t col = 0; col != row; ++col) {
227  V(row,col) *= h[row] * h[col];
228  V(col,row) = V(row,col);
229  }
230  }
231 
232  for (size_t i = 0; i != N; ++i) {
233  V(i,i) = 1.0 + lambda;
234  }
235 
236 
237  // invert Hesse matrix
238 
239  try {
240  V.invert();
241  }
242  catch (const JException& error) {
243 
244  ERROR("JGandalf: " << error.what() << endl << V << endl);
245 
246  break;
247  }
248 
249 
250  // update value and error
251 
252  for (size_t row = 0; row != N; ++row) {
253 
254  DEBUG("u[" << noshowpos << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
255 
256  for (size_t col = 0; col != N; ++col) {
257 
258  V(row,col) *= h[row] * h[col];
259 
260  get(value, parameters[row]) -= V(row,col) * get(current.result.gradient, parameters[col]);
261  }
262 
263  DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
264 
265  get(error, parameters[row]) = h[row];
266  }
267  }
268 
269  // abnormal end
270 
271  value = previous.value; // restore value
272  error = previous.error; // restore error
273 
274  V.swap(previous.V); // restore Hesse matrix
275 
276  return previous.result;
277  }
double lambda
control parameter
Definition: JGandalf.hh:290
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:281
JShower3EZ value
value
Definition: JGandalf.hh:291
static double PIVOT
minimal value diagonal element of Hesse matrix
Definition: JGandalf.hh:286
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:293
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
fit parameters
Definition: JGandalf.hh:288
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:421
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:282
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:285
struct JFIT::JGandalf::@12 previous
void reset()
Reset current parameters.
Definition: JGandalf.hh:299
#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:319
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:284
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:80
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:143
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:280
int numberOfIterations
number of iterations
Definition: JGandalf.hh:289
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:283
JShower3EZ error
error
Definition: JGandalf.hh:292
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:432
struct JFIT::JGandalf::@11 current

Member Data Documentation

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

Time window with respect to Cherenkov hypothesis [ns].

Definition at line 422 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 423 of file JShower3EZRegressor.hh.

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

Definition at line 425 of file JShower3EZRegressor.hh.

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

Definition at line 427 of file JShower3EZRegressor.hh.

PDF.

Definition at line 429 of file JShower3EZRegressor.hh.

M-Estimator function.

Definition at line 431 of file JShower3EZRegressor.hh.

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

maximal number of iterations

maximal number of iterations.

Definition at line 280 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 281 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 282 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 283 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 284 of file JGandalf.hh.

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

multiplication factor control parameter

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

fit parameters

Definition at line 288 of file JGandalf.hh.

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

number of iterations

Definition at line 289 of file JGandalf.hh.

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

control parameter

Definition at line 290 of file JGandalf.hh.

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

value

Definition at line 291 of file JGandalf.hh.

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

error

Definition at line 292 of file JGandalf.hh.

Hesse matrix.

Definition at line 293 of file JGandalf.hh.

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

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