Jpp  18.3.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::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 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 64 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 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:304

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: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
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: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:304
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:303
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 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
JShower3EZ 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
struct JFIT::JGandalf::@12 previous
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
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
JShower3EZ error
error
Definition: JGandalf.hh:304
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:444
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 292 of file JGandalf.hh.

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

maximal distance to minimum

maximal distance to minimum.

Definition at line 293 of file JGandalf.hh.

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

minimal value control parameter

Definition at line 294 of file JGandalf.hh.

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

maximal value control parameter

Definition at line 295 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 296 of file JGandalf.hh.

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

multiplication factor control parameter

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

fit parameters

Definition at line 300 of file JGandalf.hh.

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

number of iterations

Definition at line 301 of file JGandalf.hh.

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

control parameter

Definition at line 302 of file JGandalf.hh.

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

value

Definition at line 303 of file JGandalf.hh.

JShower3EZ JFIT::JGandalf< JShower3EZ >::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< JShower3EZ >::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: