Jpp
JShower3EZRegressor.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JSHOWER3EZREGRESSOR__
2 #define __JFIT__JSHOWER3EZREGRESSOR__
3 
4 #include <string>
5 
6 #include <JPhysics/JPDFTypes.hh>
7 #include <JPhysics/JPDFTable.hh>
9 #include <JPhysics/JNPETable.hh>
10 
11 #include <JTools/JFunction1D_t.hh>
13 #include <JTools/JConstants.hh>
14 
15 #include <JGeometry3D/JVector3D.hh>
16 #include <JGeometry3D/JVersor3D.hh>
17 
18 #include "JMath/JZero.hh"
19 
20 #include "JFit/JTimeRange.hh"
21 #include "JFit/JPMTW0.hh"
22 #include "JFit/JSimplex.hh"
23 #include "JFit/JMEstimator.hh"
24 #include "JFit/JRegressor.hh"
25 #include "JFit/JShower3EZ.hh"
26 
27 #include "Jeep/JMessage.hh"
28 
29 /**
30  * \file
31  * Data regression method for JFIT::JShower3Z.
32  */
33 
34 namespace JFIT {
35 
39 
40  /**
41  * Regressor function object for JShower3EZ fit using JSimplex minimiser.
42  */
43  template<>
45  public JAbstractRegressor<JShower3EZ, JSimplex>
46  {
48 
55 
61 
62  /**
63  * Default constructor
64  */
66  {}
67 
68  /**
69  * Parameterized constructor
70  *
71  * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD which
72  * will be replaced by the corresponding PDF types listed in JRegressor<JShower3Z, JGandalf>::pdf_t.
73  *
74  * \param fileDescriptor PDF file descriptor
75  */
76 
77  JRegressor(const std::string& fileDescriptor):
78  estimator(new JMEstimatorNull())
79  {
80  using namespace std;
81  using namespace JPHYSICS;
82 
83  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
84 
85  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
86 
87  try {
88 
89  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
90 
91  NOTICE("loading PDF from file " << file_name << "... " << flush);
92 
93  pdf[i].load(file_name.c_str());
94 
95  NOTICE("OK" << endl);
96 
97  pdf[i].setExceptionHandler(supervisor);
98 
99  }
100  catch(const JException& error) {
101  FATAL(error.what() << endl);
102  }
103  }
104 
105  // Add PDFs
106  for (int i = 0; i < NUMBER_OF_PDFS; i++) {
107  npe[ i ] = JNPE_t(pdf[i]);
108  }
109 
110  }
111 
112  /**
113  * Fit function.
114  * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
115  *
116  * \param shower shower
117  * \param pmt pmt
118  * \return chi2
119  */
120  double operator()(const JShower3EZ& shower, const JPMTW0& pmt) const
121  {
122  using namespace JGEOMETRY3D;
123  using namespace JPHYSICS;
124 
125  JPosition3D D(pmt.getPosition());
126  JDirection3D U(pmt.getDirection());
127 
128  D.sub(shower.getPosition());
129 
130  JVersor3D shower_dir(shower.getDirection());
131 
132  const double z = D.getDot(shower_dir);
133  const double x = D.getX() - z * shower.getDX();
134  const double y = D.getY() - z * shower.getDY();
135  const double cosDelta = z/D.getLength(); // Delta = angle between shower direction and PMT position
136 
137  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
138 
139  const double theta = U.getTheta();
140  const double phi = fabs(U.getPhi());
141 
142  JNPE_t::result_type H0 = getH0(pmt.getR()); // getH0 = Get background hypothesis value for time integrated PDF.
143  JNPE_t::result_type H1 = getH1(D.getLength(), cosDelta, theta, phi, shower.getE()); // getH1 = Get signal hypothesis value for time integrated PDF.
144 
145  if (H1.f >= Vmax_npe) {
146  H1 *= Vmax_npe / H1.f;
147  }
148 
149  H1 += H0; // now H1 is signal + background
150 
151  const bool hit = pmt.getN() != 0;
152  const double u = H1.getChi2(hit); // -log(lik)
153 
154  return estimator->getRho(u);
155  }
156 
157  /**
158  * Get background hypothesis value for time integrated PDF.
159  *
160  * \param R_Hz rate [Hz]
161  * \return hypothesis value
162  */
163  JNPE_t::result_type getH0(const double R_Hz) const
164  {
165  return JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength(), 0.0);
166  }
167 
168  /**
169  * Get signal hypothesis value for time integrated PDF.
170  *
171  * \param D PMT distance from shower [m]
172  * \param cosDelta angle between shower direction and PMT position
173  * \param theta PMT zenith angle [deg]
174  * \param phi PMT azimuth angle [deg]
175  * \param E shower energy [GeV]
176  * \return hypothesis value
177  */
178  JNPE_t::result_type getH1(const double D,
179  const double cosDelta,
180  const double theta,
181  const double phi,
182  const double E) const
183  {
185 
186  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
187 
188  if (!npe[i].empty() && D <= npe[i].getXmax()) {
189 
190  try {
191 
192  JNPE_t::result_type y1 = E * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi);
193 
194  if (y1.f > 0.0) {
195  h1 += y1;
196  }
197 
198  }
199  catch(JLANG::JException& error) {
200  ERROR(error << std::endl);
201  }
202  }
203  }
204 
205  return h1;
206  }
207 
208 
209  /**
210  * Get maximal road width of PDF.
211  *
212  * \return road width [m]
213  */
214  inline double getRmax() const
215  {
216  double xmax = 0.0;
217 
218  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
219 
220  if (!pdf[i].empty() && pdf[i].getXmax() > xmax) {
221  xmax = pdf[i].getXmax();
222  }
223 
224  }
225 
226  return xmax;
227  }
228 
229 
230  static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
231  static double Vmax_npe; //!< Maximal integral of PDF [npe]
232 
233  static const int NUMBER_OF_PDFS = 1;
234 
235  static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
236 
237  JPDF_t pdf[NUMBER_OF_PDFS]; //!< PDF
238  JNPE_t npe[NUMBER_OF_PDFS]; //!< PDF
239 
241  };
242 
243  /**
244  * PDF types.
245  */
247 
248 
249  /**
250  * Default values.
251  */
253  double JRegressor<JShower3EZ, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
254 }
255 
256 #endif
JFIT::JRegressor< JShower3EZ, JSimplex >::T_ns
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
Definition: JShower3EZRegressor.hh:230
JFIT::JRegressor< JShower3EZ, JSimplex >::getH1
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.
Definition: JShower3EZRegressor.hh:178
JFIT::JRegressor< JShower3EZ, JSimplex >::JNPE_t
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
Definition: JShower3EZRegressor.hh:60
JVector3D.hh
JFIT::JAbstractRegressor
Abstract class for global fit method.
Definition: JRegressor.hh:75
JFIT
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
JFIT::JRegressor< JShower3EZ, JSimplex >::JPDFMaplist_t
JTOOLS::JMapList< JTOOLS::JPolint1FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint1FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap > > > > JPDFMaplist_t
Definition: JShower3EZRegressor.hh:53
JMessage.hh
JZero.hh
JFIT::JShower3EZ
Data structure for fit of straight line in positive z-direction with energy.
Definition: JShower3EZ.hh:25
JFIT::JEnergy::getE
double getE() const
Get energy.
Definition: JEnergy.hh:144
JPHYSICS::JNPETable< double, double, JNPEMaplist_t >::result_type
multifunction_t::result_type result_type
Definition: JNPETable.hh:60
JFIT::JPMTW0::getR
double getR() const
Get rate.
Definition: JPMTW0.hh:56
JFIT::JRegressor< JShower3EZ, JSimplex >::JFunction1D_t
JTOOLS::JSplineFunction1S_t JFunction1D_t
Definition: JShower3EZRegressor.hh:49
JGEOMETRY3D::JDirection3D::getDirection
const JDirection3D & getDirection() const
Get direction.
Definition: JDirection3D.hh:106
JFIT::JRegressor< JShower3EZ, JSimplex >::JRegressor
JRegressor()
Default constructor.
Definition: JShower3EZRegressor.hh:65
JRegressor.hh
JFIT::ERROR
Definition: JFitStatus.hh:18
JFIT::JSimplex
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition: JSimplex.hh:42
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:706
JFIT::JRegressor< JShower3EZ, JSimplex >::JNPEMaplist_t
JTOOLS::JMapList< JTOOLS::JPolint1FunctionalMapH, JTOOLS::JMapList< JTOOLS::JPolint1FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap > > > > JNPEMaplist_t
Definition: JShower3EZRegressor.hh:59
JPHYSICS::DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition: JPDFTypes.hh:40
JSimplex.hh
JPDFToolkit.hh
JTOOLS::JRange< double >
JPHYSICS
Auxiliary classes and methods for calculation of PDF and muon energy loss.
Definition: JAbstractMedium.hh:9
JGEOMETRY3D::JDirection3D
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
JTOOLS::JPolint0FunctionalGridMap
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
Definition: JFunctionalMap_t.hh:82
JTOOLS::JSplineFunction1S_t
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.
Definition: JFunction1D_t.hh:51
JGEOMETRY3D::JVersor3D
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JFIT::JPMTW0
Auxiliary class for handling PMT geometry, rate and response.
Definition: JPMTW0.hh:22
JPHYSICS::JNPETable< double, double, JNPEMaplist_t >
JTOOLS::JMapList
Map list.
Definition: JMapList.hh:24
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
JPMTW0.hh
JFIT::JRegressor< JShower3EZ, JSimplex >::operator()
double operator()(const JShower3EZ &shower, const JPMTW0 &pmt) const
Fit function.
Definition: JShower3EZRegressor.hh:120
JFunction1D_t.hh
JFIT::JMEstimatorNull
Null M-estimator.
Definition: JMEstimator.hh:53
JFIT::JPMTW0::getN
int getN() const
Get number of hits.
Definition: JPMTW0.hh:67
JTOOLS::JPolint1FunctionalMapH
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:145
JFIT::JRegressor< JShower3EZ, JSimplex >::Vmax_npe
static double Vmax_npe
Maximal integral of PDF [npe].
Definition: JShower3EZRegressor.hh:231
JTOOLS::JPolint1FunctionalMap
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Definition: JFunctionalMap_t.hh:55
JConstants.hh
JShower3EZ.hh
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JFunctionalMap_t.hh
JGEOMETRY3D::JVersor3D::getDot
double getDot(const JVersor3D &versor) const
Get dot product.
Definition: JVersor3D.hh:153
JPHYSICS::SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition: JPDFTypes.hh:41
JFIT::JRegressor
Template definition of a data regressor of given model.
Definition: JRegressor.hh:68
JNPETable.hh
JEEP::getFilename
std::string getFilename(const std::string &file_name)
Get file name part, i.e.
Definition: JeepToolkit.hh:86
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t >
JTimeRange.hh
JGEOMETRY3D::JRotation3Z
Rotation around Z-axis.
Definition: JRotation3D.hh:85
JPHYSICS::JPDFType_t
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
JFIT::JRegressor< JShower3EZ, JSimplex >::estimator
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Definition: JShower3EZRegressor.hh:240
JFIT::JRegressor< JShower3EZ, JSimplex >::getH0
JNPE_t::result_type getH0(const double R_Hz) const
Get background hypothesis value for time integrated PDF.
Definition: JShower3EZRegressor.hh:163
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
JMEstimator.hh
std
Definition: jaanetDictionary.h:36
JGEOMETRY3D::JVersor3Z::getDY
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:156
JGEOMETRY3D
Auxiliary classes and methods for 3D geometrical objects and operations.
Definition: JAngle3D.hh:18
JLANG::JException::what
virtual const char * what() const
Get error message.
Definition: JException.hh:48
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:94
JVersor3D.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JFIT::JRegressor< JShower3EZ, JSimplex >::getRmax
double getRmax() const
Get maximal road width of PDF.
Definition: JShower3EZRegressor.hh:214
JPDFTypes.hh
JFIT::JShower3Z::getDirection
const JVersor3Z & getDirection() const
Get direction.
Definition: JVersor3Z.hh:79
JFIT::JRegressor< JShower3EZ, JSimplex >::JPDF_t
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Definition: JShower3EZRegressor.hh:54
JLANG::JException
General exception.
Definition: JException.hh:23
JPDFTable.hh
JFIT::JRegressor< JShower3EZ, JSimplex >::JRegressor
JRegressor(const std::string &fileDescriptor)
Parameterized constructor.
Definition: JShower3EZRegressor.hh:77
JGEOMETRY3D::JVersor3Z::getDX
double getDX() const
Get x direction.
Definition: JVersor3Z.hh:145
JLANG::JSharedPointer
The template JSharedPointer class can be used to share a pointer to an object.
Definition: JSharedPointer.hh:28