Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEnergyRegressor.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JENERGYREGRESSOR__
2 #define __JFIT__JENERGYREGRESSOR__
3 
4 #include <string>
5 
7 #include "JPhysics/JPDFTypes.hh"
8 #include "JPhysics/JPDFTable.hh"
9 #include "JPhysics/JNPETable.hh"
10 #include "JPhysics/JPDFToolkit.hh"
11 #include "JPhysics/JConstants.hh"
12 #include "JTools/JFunction1D_t.hh"
14 #include "JGeometry3D/JAxis3D.hh"
15 #include "JFit/JRegressor.hh"
16 #include "JFit/JEnergy.hh"
17 #include "JFit/JNPE.hh"
18 #include "JFit/JNPEHit.hh"
19 #include "JFit/JMEstimator.hh"
20 #include "JFit/JFitToolkit.hh"
21 #include "JFit/JTimeRange.hh"
22 #include "Jeep/JMessage.hh"
23 
24 
25 /**
26  * \file
27  * Data regression method for JFIT::JEnergy.
28  * \author mdejong
29  */
30 namespace JFIT {}
31 namespace JPP { using namespace JFIT; }
32 
33 namespace JFIT {
34 
41 
42 
43  /**
44  * Regressor function object for JEnergy fit.
45  */
46  template<>
47  struct JRegressor<JEnergy> :
48  public JAbstractRegressor<JEnergy>
49  {
52 
57 
58 
59  /**
60  * Constructor.
61  *
62  * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD
63  * which will be replaced by the PDF types listed in JRegressor<JEnergy, JSimplex>::pdf_t.
64  *
65  * \param fileDescriptor PDF file descriptor
66  */
67  JRegressor(const std::string& fileDescriptor) :
68  estimator(new JMEstimatorNormal())
69  {
70  using namespace std;
71  using namespace JPP;
72 
73  typedef JSplineFunction1D_t JFunction1D_t;
74  typedef JPDFTable<JFunction1D_t, JNPEMaplist_t> JPDF_t;
75 
76  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
77 
78  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
79 
80  try {
81 
82  JPDF_t pdf;
83 
84  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
85 
86  NOTICE("loading PDF from file " << file_name << "... " << flush);
87 
88  pdf.load(file_name.c_str());
89 
90  NOTICE("OK" << endl);
91 
92  pdf.setExceptionHandler(supervisor);
93 
94  if (!is_bremsstrahlung(pdf_t[i]))
95  Y1.push_back(JNPE_t(pdf));
96  else
97  YB.push_back(JNPE_t(pdf));
98  }
99  catch(const JLANG::JException& error) {
100  FATAL(error.what() << endl);
101  }
102  }
103 
104  // Add PDFs
105 
106  if (Y1.size() == 2) { Y1[1].add(Y1[0]); Y1.erase(Y1.begin()); }
107  if (YB.size() == 2) { YB[1].add(YB[0]); YB.erase(YB.begin()); }
108  }
109 
110 
111  /**
112  * Fit function.
113  * This method is used to determine chi2 of given number of photo-electrons for given energy of muon.
114  *
115  * \param x energy
116  * \param npe npe
117  * \return chi2
118  */
119  inline double operator()(const JEnergy& x, const JNPEHit& npe) const
120  {
121  const double E = x.getE();
122  const double u = npe.getChi2(E);
123 
124  return estimator->getRho(u);
125  }
126 
127 
128  /**
129  * Create data structure for handling light yields for PMT.
130  * Note that the PMT geometry should be relative to the muon trajectory,
131  * conform method JGEOMETRY3D::JAxis3D::transform.
132  *
133  * \param axis PMT axis
134  * \param R_Hz singles rate [Hz]
135  * \return light yields
136  */
137  inline JNPE getNPE(const JAxis3D& axis,
138  const double R_Hz) const
139  {
140  using namespace std;
141  using namespace JPP;
142 
143  const double x = axis.getX();
144  const double y = axis.getY();
145  const double R = sqrt(x*x + y*y);
146  const double z = axis.getZ() - R / getTanThetaC();
147 
148  const double theta = axis.getTheta();
149  const double phi = fabs(axis.getPhi());
150 
151  const double yA = getNPE(Y1, R, theta, phi);
152  const double yB = getNPE(YB, R, theta, phi);
153 
154  return JNPE(JK40(T_ns, R_Hz), yA, yB, z);
155  }
156 
157 
158  /**
159  * Create data structure for handling light yields for set of PMTs.
160  * Note that the PMT geometry should be relative to the muon trajectory,
161  * conform method JGEOMETRY3D::JAxis3D::transform.
162  *
163  * \param __begin begin of PMTs
164  * \param __end end of PMTs
165  * \param rates_Hz K40 singles and multiples rates [Hz]
166  * \return light yields
167  */
168  template<class T>
169  inline JNPE getNPE(T __begin,
170  T __end,
171  const JK40Rates& rates_Hz) const
172  {
173  using namespace std;
174  using namespace JPP;
175 
176  double yA = 0.0;
177  double yB = 0.0;
178  double z = 0.0;
179 
180  for (T i = __begin; i != __end; ++i) {
181 
182  const JNPE npe = this->getNPE(*i, rates_Hz.getSinglesRate());
183  const double x = i->getX();
184  const double y = i->getY();
185  const double R = sqrt(x*x + y*y);
186 
187  yA += npe.getYA();
188  yB += npe.getYB();
189  z += i->getZ() - R / getTanThetaC();
190  }
191 
192  const int N = distance(__begin, __end);
193 
194  return JNPE(JK40(T_ns, N, rates_Hz), yA, yB, z / N);
195  }
196 
197 
198  /**
199  * Get maximal road width of NPE.
200  *
201  * \return road width [m]
202  */
203  inline double getRmax() const
204  {
205  return std::max(getRmax(Y1),
206  getRmax(YB));
207  }
208 
209 
210  std::vector<JNPE_t> Y1; //!< light from muon
211  std::vector<JNPE_t> YB; //!< light from EM showers
212 
213  static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
214 
215 
216  static const int NUMBER_OF_PDFS = 4;
217 
218  static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
219 
221 
222 
223  protected:
224  /**
225  * Get maximal road width of PDF.
226  *
227  * \param NPE NPE tables
228  * \return road width [m]
229  */
230  static inline double getRmax(const std::vector<JNPE_t>& NPE)
231  {
232  double xmax = 0.0;
233 
234  for (std::vector<JNPE_t>::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
235  if (!i->empty() && i->getXmax() > xmax) {
236  xmax = i->getXmax();
237  }
238  }
239 
240  return xmax;
241  }
242 
243 
244  /**
245  * Get number of photo-electrons.
246  *
247  * \param NPE NPE tables
248  * \param R distance between muon and PMT [m]
249  * \param theta zenith angle orientation PMT [rad]
250  * \param phi azimuth angle orientation PMT [rad]
251  * \return number of photo-electrons
252  */
253  static inline double getNPE(const std::vector<JNPE_t>& NPE,
254  const double R,
255  const double theta,
256  const double phi)
257  {
258  using JTOOLS::get_value;
259 
260  double npe = 0.0;
261 
262  for (std::vector<JNPE_t>::const_iterator i = NPE.begin(); i != NPE.end(); ++i) {
263 
264  if (R <= i->getXmax()) {
265 
266  try {
267 
268  const double y = get_value((*i)(std::max(R, i->getXmin()), theta, phi));
269 
270  if (y > 0.0) {
271  npe += y;
272  }
273  }
274  catch(const JLANG::JException& error) {
275  ERROR(error << std::endl);
276  }
277  }
278  }
279 
280  return npe;
281  }
282  };
283 
284 
285  /**
286  * PDF types.
287  */
292 
293 
294  /**
295  * Time range.
296  */
298 }
299 
300 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
Template definition of a data regressor of given model.
Definition: JRegressor.hh:68
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
double getE() const
Get energy.
Definition: JEnergy.hh:170
General exception.
Definition: JException.hh:23
General purpose data regression method.
Auxiliary methods for PDF calculations.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
static double getNPE(const std::vector< JNPE_t > &NPE, const double R, const double theta, const double phi)
Get number of photo-electrons.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
double getYB() const
Get light yield due to bremsstrahlung.
Definition: JNPE.hh:80
direct light from EM showers
Definition: JPDFTypes.hh:32
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JNPE getNPE(T __begin, T __end, const JK40Rates &rates_Hz) const
Create data structure for handling light yields for set of PMTs.
Auxiliary class for converting various rates to expectation values of the number of hits within a giv...
Definition: JK40.hh:29
double getChi2(const double E_GeV) const
Get chi2.
Definition: JNPEHit.hh:87
double getPhi() const
Get phi angle.
Definition: JVersor3D.hh:144
direct light from muon
Definition: JPDFTypes.hh:29
Various implementations of functional maps.
double getYA() const
Get light yield due to muon itself.
Definition: JNPE.hh:69
Numbering scheme for PDF types.
Axis object.
Definition: JAxis3D.hh:38
The template JSinglePointer class can be used to hold a pointer to an object.
bool is_bremsstrahlung(const int pdf)
Test if given PDF type corresponds to Cherenkov light from Bremsstrahlung.
Definition: JPDFTypes.hh:163
JNPE getNPE(const JAxis3D &axis, const double R_Hz) const
Create data structure for handling light yields for PMT.
scattered light from muon
Definition: JPDFTypes.hh:30
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
JLANG::JSinglePointer< JMEstimator > estimator
M-Estimator function.
JRegressor(const std::string &fileDescriptor)
Constructor.
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
std::vector< JNPE_t > Y1
light from muon
do set_variable OUTPUT_DIRECTORY $WORKDIR T
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
double operator()(const JEnergy &x, const JNPEHit &npe) const
Fit function.
double getTheta() const
Get theta angle.
Definition: JVersor3D.hh:128
Auxiliary class for handling various light yields.
Definition: JNPE.hh:27
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
Physics constants.
Auxiliary class for simultaneously handling light yields and response of PMT.
Definition: JNPEHit.hh:19
scattered light from EM showers
Definition: JPDFTypes.hh:33
double getY() const
Get y position.
Definition: JVector3D.hh:104
Normal M-estimator.
Definition: JMEstimator.hh:66
int debug
debug level
Definition: JSirene.cc:63
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:40
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
double getRmax() const
Get maximal road width of NPE.
double getSinglesRate() const
Get singles rate.
Definition: JK40Rates.hh:97
virtual const char * what() const override
Get error message.
Definition: JException.hh:48
static double getRmax(const std::vector< JNPE_t > &NPE)
Get maximal road width of PDF.
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint1FunctionalGridMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JNPEMaplist_t
double getX() const
Get x position.
Definition: JVector3D.hh:94
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:88
Data structure for fit of energy.
Definition: JEnergy.hh:28
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
double getNPE(const Hit &hit)
Get true charge of hit.
double u[N+1]
Definition: JPolint.hh:739
Abstract class for global fit method.
Definition: JRegressor.hh:75
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:998
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Maximum likelihood estimator (M-estimators).
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
then usage $script[input file[working directory[option]]] nWhere option can be E
Definition: JMuonPostfit.sh:37
Auxiliary class for handling debug parameter within a class.
Definition: JMessage.hh:44
std::vector< JNPE_t > YB
light from EM showers