Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPoint4DRegressor.hh
Go to the documentation of this file.
1 #ifndef __JPOINT4DREGRESSOR__
2 #define __JPOINT4DREGRESSOR__
3 
4 #include <string>
5 #include <iostream>
6 #include <cmath>
7 
8 #include "JPhysics/JPDFTypes.hh"
9 #include "JPhysics/JPDFTable.hh"
10 #include "JPhysics/JPDFToolkit.hh"
11 #include "JPhysics/JNPETable.hh"
12 
13 #include "JTools/JFunction1D_t.hh"
15 #include "JTools/JConstants.hh"
16 
17 #include "JTrigger/JHitL0.hh"
18 
19 #include "JGeometry3D/JVector3D.hh"
20 #include "JGeometry3D/JVersor3D.hh"
21 
22 #include "JMath/JZero.hh"
23 
24 #include "JFit/JTimeRange.hh"
25 #include "JFit/JPMTW0.hh"
26 #include "JFit/JSimplex.hh"
27 #include "JFit/JGandalf.hh"
28 #include "JFit/JMEstimator.hh"
29 #include "JFit/JRegressor.hh"
30 #include <JFit/JPoint4D.hh>
31 #include <JFit/JShower3Z.hh>
32 
33 #include "Jeep/JMessage.hh"
34 
35 #include <TH2D.h>
36 
37 /**
38  * \author adomi
39  */
40 
41 namespace JFIT{
42 
43  /**
44  * Regressor function object for JPoint4D fit using JSimplex minimiser.
45  */
46  template<>
48  public JAbstractRegressor<JPoint4D, JSimplex>
49  {
51 
52  /**
53  * Constructor.
54  */
55  JRegressor(const bool __prefit):
56  estimator(new JMEstimatorTukey(30.0)),
57  prefit(__prefit)
58  {}
59 
60 
61  /**
62  * Constructor for a 2D histogram as PDF
63  */
64  JRegressor(TH2D *hPDF2D, bool __prefit):
65  estimator(new JMEstimatorTukey(30.0)),
66  prefit(__prefit)
67  {
68  hpdf2d[0] = (TH2D*)hPDF2D->Clone();
69  hpdf2d[0]->Smooth();
70  }
71 
72 
73  /* Fit Function
74  * This method is used to determine the chi2 of given hit with respect
75  * to Shower's vertex (brightest point)
76  *
77  * \param vertex Shower's vertex
78  * \param hit L0 hit
79  * \return -log(likelihood)
80  */
81  template<class JHit_t>
82  double operator()(const JPoint4D& vertex, const JHit_t& hit) const
83  {
84  using namespace std;
85  using namespace JGEOMETRY3D;
86  using namespace JTOOLS;
87 
88  const JPosition3D hit_pos(hit.getPosition());
89  const double D = hit_pos.getDistance(vertex);
90  const double t_res = hit.getT() - vertex.getT(hit_pos);
91  const JDirection3D photonDir(hit_pos - vertex.getPosition());
92  const double cosT = photonDir.getDot(hit.getDirection());
93  const double cosT_th = d_ref / sqrt(d_ref*d_ref + D*D);
94 
95  double LogLik = 0;
96 
97  if(prefit){
98 
99  LogLik = evaluatePrefit(t_res, cosT, cosT_th);
100 
101  } else {
102 
103  LogLik = evaluateFit(t_res, D);
104 
105  }
106 
107  return LogLik;
108  }
109 
110 
111  double evaluatePrefit(const double t_res, const double cosT, const double cosT_th) const
112  {
113 
114  double p_T = 0;
115 
116  if(cosT < 0){
117  p_T = 0;
118  } else if(cosT > cosT_th){
119  p_T = p_max;
120  } else {
121  p_T = p_max * cosT/cosT_th;
122  }
123 
124  double lik = 1 / (sqrt(a*a + t_res*t_res) + p_T);
125  double u = -log(lik);
126 
127  return estimator->getRho(u);
128  }
129 
130 
131  double evaluateFit(const double t_res, const double D) const
132  {
133 
134  int x_bin = hpdf2d[0]->GetXaxis()->FindBin(t_res);
135  int y_bin = hpdf2d[0]->GetYaxis()->FindBin(D);
136 
137  int nbinsX = hpdf2d[0]->GetXaxis()->GetNbins();
138  int nbinsY = hpdf2d[0]->GetYaxis()->GetNbins();
139 
140  double lik = 0;
141 
142  if(x_bin > 1 && x_bin < nbinsX){
143  if(y_bin > 1 && y_bin < nbinsY){
144  lik = hpdf2d[0]->Interpolate(t_res, D);
145  }
146  }
147 
148  if(lik == 0){
149  lik = 0.000004; // value of K40 bg
150  }
151 
152  double u = -log(lik);
153  return estimator->getRho(u);
154 
155  }
156 
157  static const double a; //!< constant to prevent the likelihood going to infinity
158  static const double d_ref;
159  static const double p_max;
160 
161  static const int NUMBER_OF_PDFS = 1;
162  TH2D *hpdf2d[NUMBER_OF_PDFS]; //!< PDF 2D
163 
164  JLANG::JSharedPointer<JMEstimator> estimator; //!< M-Estimator function
165  const bool prefit;
166  };
167 
168  const double JRegressor<JPoint4D, JSimplex>::a = 2.0; //!< constant to prevent the likelihood going to infinity
169  const double JRegressor<JPoint4D, JSimplex>::d_ref = 10.0;
170  const double JRegressor<JPoint4D, JSimplex>::p_max = 100;
171 
172 }
173 
174 #endif
Template definition of a data regressor of given model.
Definition: JRegressor.hh:66
General purpose data regression method.
Auxiliary methods for PDF calculations.
Data structure for direction in three dimensions.
Definition: JDirection3D.hh:32
double evaluateFit(const double t_res, const double D) const
Data structure for vertex fit.
Definition: JPoint4D.hh:22
JRegressor(const bool __prefit)
Constructor.
Definition of zero value for any class.
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:266
static const double a
constant to prevent the likelihood going to infinity
double getDot(const JAngle3D &angle) const
Get dot product.
Various implementations of functional maps.
Basic data structure for L0 hit.
Numbering scheme for PDF types.
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
Constants.
The template JSharedPointer class can be used to share a pointer to an object.
JRegressor(TH2D *hPDF2D, bool __prefit)
Constructor for a 2D histogram as PDF.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
General purpose messaging.
double evaluatePrefit(const double t_res, const double cosT, const double cosT_th) const
double operator()(const JPoint4D &vertex, const JHit_t &hit) const
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:144
Abstract class for global fit method.
Definition: JRegressor.hh:73
Tukey&#39;s biweight M-estimator.
Definition: JMEstimator.hh:105
Maximum likelihood estimator (M-estimators).