Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JShowerBrightPointRegressor.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__
2 #define __JFIT__JSHOWERBRIGHTPOINTREGRESSOR__
3 
4 #include "JPhysics/JPDFTypes.hh"
5 #include "JPhysics/JPDFTable.hh"
7 #include "JPhysics/JNPETable.hh"
8 #include "JPhysics/JConstants.hh"
9 
10 #include "JTools/JResult.hh"
11 
12 #include "JMath/JZero.hh"
13 
14 #include "JFit/JGandalf.hh"
15 #include "JFit/JPoint4E.hh"
16 #include "JFit/JMEstimator.hh"
17 #include "JFit/JRegressor.hh"
18 #include "JFit/JFitToolkit.hh"
19 #include "JFit/JTimeRange.hh"
20 
21 #include "Jeep/JMessage.hh"
22 
23 /**
24  * \file
25  * Data regression method for JFIT::JPoint4E from a bright point isoptropic emission PDF.
26  *
27  * \author adomi, vcarretero
28  */
29 
30 namespace JFIT {}
31 namespace JPP { using namespace JFIT; }
32 
33 namespace JFIT {
34 
38  using JTOOLS::get_value;
39 
40  /**
41  * Regressor function object for JPoint4E fit using JGandalf minimiser.
42  */
43  template<>
45  public JAbstractRegressor<JPoint4E, JGandalf>
46  {
48 
53 
54 
55  /**
56  * Parameterized constructor
57  *
58  * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
59  * will be replaced by the corresponding PDF types.
60  *
61  * \param fileDescriptor PDF file descriptor
62  * \param TTS TTS [ns]
63  * \param numberOfPoints number of points for Gauss-Hermite integration of TTS
64  * \param epsilon precision for Gauss-Hermite integration of TTS
65  */
66  JRegressor(const std::string& fileDescriptor,
67  const double TTS,
68  const int numberOfPoints = 25,
69  const double epsilon = 1.0e-10)
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  }
111 
112  /**
113  * Fit function.
114  * This method is used to determine the chi2 and gradient of given hit with respect a bright point emitting isotropically
115  *
116  * JHit_t refers to a data structure which should have the following member methods:
117  * - double getX(); // [m]
118  * - double getY(); // [m]
119  * - double getZ(); // [m]
120  * - double getDX(); // [u]
121  * - double getDY(); // [u]
122  * - double getDZ(); // [u]
123  * - double getT(); // [ns]
124  *
125  * \param vx shower vertex
126  * \param hit hit
127  * \return chi2 and gradient
128  */
129  template<class JHit_t>
130  result_type operator()(const JPoint4E& vx, const JHit_t& hit) const
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, vx.getE()); // 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 = JPoint4E(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  0); // d(ct)/d(E) //Set to 0 for now, not a parameter to fit
169 
170  result.gradient.mul(getInverseSpeedOfLight() * (H1.getDerivativeOfChi2() -
171  H0.getDerivativeOfChi2())); // x d(chi2)/d(ct1)
172 
173  return result;
174 
175  }
176 
177  /**
178  * Get background hypothesis value for time differentiated PDF.
179  *
180  * \param R_Hz rate [Hz]
181  * \param t1 time [ns]
182  * \return hypothesis value
183  */
184  JPDF_t::result_type getH0(const double R_Hz,
185  const double t1) const
186  {
187  using namespace JPP;
188 
189  return JPDF_t::result_type(R_Hz * 1e-9, t1, T_ns);
190  }
191 
192  /**
193  * Get signal hypothesis value for bright point emission PDF.
194  *
195  * \param D hit distance from shower vertex [m]
196  * \param ct cosine of the HIT angle
197  * \param t arrival time of the light
198  * \param E shower energy [GeV]
199  * \return hypothesis value
200  */
202  const double ct,
203  const double t,
204  const double E) const
205  {
206  using namespace JPP;
207 
209 
210  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
211 
212  if (!pdf[i].empty() && D <= pdf[i].getXmax()) {
213 
214  try {
215 
216  JPDF_t::result_type y1 = E * pdf[i](std::max(D, pdf[i].getXmin()), ct, t);
217 
218  // safety measures
219 
220  if (y1.f <= 0.0) {
221  y1.f = 0.0;
222  y1.fp = 0.0;
223  }
224 
225  if (y1.v <= 0.0) {
226  y1.v = 0.0;
227  }
228 
229  h1 += y1;
230 
231  }
232  catch(JLANG::JException& error) {
233  ERROR(error << std::endl);
234  }
235  }
236  }
237 
238  return h1;
239  }
240 
241  /**
242  * Get maximal road width of PDF.
243  *
244  * \return road width [m]
245  */
246  inline double getRmax() const
247  {
248  using namespace JPP;
249 
250  double xmax = 0.0;
251 
252  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
253 
254  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
255  xmax = pdf[i].getXmax();
256  }
257 
258  }
259 
260  return xmax;
261  }
262 
263  static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
264  static double Vmax_npe; //!< Maximal integral of PDF [npe]
265 
266  static const int NUMBER_OF_PDFS = 2;
267 
268  static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
269 
270  JPDF_t pdf[NUMBER_OF_PDFS]; //!< PDF
271  };
272 
273 
274  /**
275  * PDF types.
276  */
279 
280  /**
281  * Default values.
282  */
284  double JRegressor<JPoint4E, JGandalf>::Vmax_npe = std::numeric_limits<double>::max();
285 
286 }
287 
288 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMapList_t > JPDF_t
Template definition of a data regressor of given model.
Definition: JRegressor.hh:70
const double xmax
Definition: JQuadrature.cc:24
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type...
General exception.
Definition: JException.hh:24
Type definition of a 2nd degree polynomial interpolation based on a JMap implementation.
General purpose data regression method.
Auxiliary methods for PDF calculations.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
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.
This include file containes various data structures that can be used as specific return types for the...
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Data structure for vertex fit.
Definition: JPoint4D.hh:22
double getE() const
Get energy.
Definition: JPoint4E.hh:52
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
direct light from bright point
Definition: JPDFTypes.hh:42
static double Vmax_npe
Maximal integral of PDF [npe].
Definition of zero value for any class.
Numbering scheme for PDF types.
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
JPDF_t::result_type getH0(const double R_Hz, const double t1) const
Get background hypothesis value for time differentiated PDF.
minimiser_type::result_type result_type
Definition: JRegressor.hh:82
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
Data structure for vertex fit.
Definition: JPoint4E.hh:22
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Physics constants.
double getRmax() const
Get maximal road width of PDF.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:24
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
JRegressor(const std::string &fileDescriptor, const double TTS, const int numberOfPoints=25, const double epsilon=1.0e-10)
Parameterized constructor.
General purpose messaging.
scattered light from bright point
Definition: JPDFTypes.hh:43
JTOOLS::JMAPLIST< JTOOLS::JPolint2FunctionalMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JPDFMapList_t
#define FATAL(A)
Definition: JMessage.hh:67
Fit method based on the Levenberg-Marquardt method.
Definition: JGandalf.hh:84
virtual const char * what() const override
Get error message.
Definition: JException.hh:64
const double getSpeedOfLight()
Get speed of light.
const double getInverseSpeedOfLight()
Get inverse speed of light.
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
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
transformablemultifunction_type::result_type result_type
Definition: JPDFTable.hh:50
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:147
Abstract class for global fit method.
Definition: JRegressor.hh:77
result_type operator()(const JPoint4E &vx, const JHit_t &hit) const
Fit function.
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
Maximum likelihood estimator (M-estimators).
const double epsilon
Definition: JQuadrature.cc:21
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).