Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
JShowerBjorkenYRegressor.hh
Go to the documentation of this file.
1#ifndef __JFIT__JSHOWERBJORKENYREGRESSOR__
2#define __JFIT__JSHOWERBJORKENYREGRESSOR__
3
8
12#include "JTools/JRange.hh"
13#include "JTools/JResult.hh"
17
21
23#include "JMath/JZero.hh"
24
25#include "JFit/JTimeRange.hh"
26#include "JFit/JPMTW0.hh"
27#include "JFit/JSimplex.hh"
28#include "JFit/JMEstimator.hh"
29#include "JFit/JRegressor.hh"
30#include "JFit/JShowerEH.hh"
31#include "JFit/JFitToolkit.hh"
32
33#include "Jeep/JMessage.hh"
34
35/**
36 * \file
37 * Data regression method for JFIT::JShowerEH.
38 * \author adomi
39 */
40
41namespace JFIT {
42
49
50 /**
51 * Regressor function object for JShowerEH fit using JSimplex minimiser.
52 */
53 template<>
54 struct JRegressor<JShowerEH, JSimplex> :
55 public JAbstractRegressor<JShowerEH, JSimplex>
56 {
57 using JAbstractRegressor<JShowerEH, JSimplex>::operator();
58
59 typedef JTOOLS::JSplineFunction1S_t JFunction1D_t;
65
71
72
73 // Part for the Isotropic PDF
75 JTOOLS::JPolint1FunctionalGridMap>::maplist JMapList_t2;
77
79 JTOOLS::JPolint1FunctionalGridMap>::maplist JNPEMapList_t2;
81
82
83 /**
84 * Parameterized constructor
85 *
86 * The PDF file descriptor should contain the wild card character JPHYSICS::WILDCARD which
87 * will be replaced by the corresponding PDF types listed in JRegressor<JShower3Z, JGandalf>::pdf_t.
88 *
89 * \param fileDescriptor PDF file descriptor
90 */
91
92 JRegressor(const std::string& fileDescriptor):
93 estimator(new JMEstimatorNull())
94 {
95 using namespace std;
96 using namespace JPP;
97
98 const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
99 const JPDF_t2::JSupervisor supervisor2(new JPDF_t2::JDefaultResult(JMATH::zero));
100
101 for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
102
103 try {
104
105 JPDF_t pdf;
106 JPDF_t2 pdf2;
107
108 const string file_name = getFilename(fileDescriptor, pdf_t[i]);
109
110 NOTICE("loading PDF from file " << file_name << "... " << flush);
111
112 if(i < 2){
113
114 pdf.load(file_name.c_str());
115
116 pdf.setExceptionHandler(supervisor);
117
118 npe[ i ] = JNPE_t(pdf);
119
120 } else {
121
122 pdf2.load(file_name.c_str());
123
124 NOTICE("OK" << endl);
125
126 pdf2.setExceptionHandler(supervisor2);
127
128 npe2[ i-2 ] = JNPE_t2(pdf2);
129
130 }
131 }
132 catch(const JException& error) {
133 FATAL(error.what() << endl);
134 }
135 }
136
137 // Add PDFs
138 for (int i = 1; i < (NUMBER_OF_PDFS-2); i += 2) {
139
140 npe[ i ].add(npe[i-1]);
141
142 JNPE_t buffer;
143
144 npe[i-1].swap(buffer);
145
146 npe2[ i ].add(npe2[i-1]);
147
148 JNPE_t2 buffer2;
149
150 npe2[i-1].swap(buffer2);
151
152 }
153 }
154
155 /**
156 * Fit function.
157 * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
158 *
159 * \param shower shower
160 * \param pmt pmt
161 * \return chi2
162 */
163 double operator()(const JShowerEH& shower, const JPMTW0& pmt) const
164 {
165 using namespace JPP;
166 using namespace std;
167
168 JPosition3D D(pmt.getPosition());
169 JDirection3D U(pmt.getDirection());
170
171 D.sub(shower.getPosition());
172
173 double ct = U.getDot(D) / D.getLength();
174
175 JVersor3D shower_dir(shower.getDirection());
176
177 const double z = D.getDot(shower_dir);
178 const double x = D.getX() - z * shower.getDX();
179 const double y = D.getY() - z * shower.getDY();
180 const double cosDelta = z/D.getLength(); // Delta = angle between shower direction and PMT position
181
182 U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
183
184 const double theta = U.getTheta();
185 const double phi = fabs(U.getPhi());
186
187 double H0 = getH0(pmt.getR()); // background
188 double H1 = getH1(D.getLength(), ct, cosDelta, theta, phi,
189 shower.getEem(), shower.getEh(), shower.getBy()); // signal
190
191 if (H1 >= Vmax_npe) {
192 H1 *= Vmax_npe / H1;
193 }
194
195 H1 += H0; // now H1 is signal + background
196
197 const bool hit = pmt.getN() != 0;
198 const double u = getChi2(H1, hit); // -log(lik)
199
200 return estimator->getRho(u);
201 }
202
203 /**
204 * Get background hypothesis value for time integrated PDF.
205 *
206 * \param R_Hz rate [Hz]
207 * \return hypothesis value
208 */
209 double getH0(const double R_Hz) const
210 {
211 return get_value(JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength()));
212 }
213
214 /**
215 * Get signal hypothesis value for time integrated PDF.
216 *
217 * \param D PMT distance from shower [m]
218 * \param ct angle between shower direction and PMT position
219 * \param cosDelta angle between shower direction and PMT position
220 * \param theta PMT zenith angle [deg]
221 * \param phi PMT azimuth angle [deg]
222 * \param Eem EM shower energy [GeV]
223 * \param Eh H shower energy [GeV]
224 * \param Y Bjorken Y
225 * \return hypothesis value
226 */
227 double getH1(const double D,
228 const double ct,
229 const double cosDelta,
230 const double theta,
231 const double phi,
232 const double Eem,
233 const double Eh,
234 const double Y) const
235 {
236
237 double h1 = 0;
238
239 for (int i = 0; i != (NUMBER_OF_PDFS-1); ++i) {
240
241 if (!npe[i].empty() && D <= npe[i].getXmax() && !npe2[i].empty() && D <= npe2[i].getXmax()) {
242
243 try {
244
245 JNPE_t::result_type P_em;
246 JNPE_t2::result_type P_h;
247
248 P_em = fabs(Eem) * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi);
249
250 P_h = fabs(Eh) * npe2[i](std::max(D, npe2[i].getXmin()), ct);
251
252 double y1 = get_value(P_em) + get_value(P_h);
253
254 if(y1 > 0.0){
255 h1 += y1;
256 }
257
258 }
259 catch(JLANG::JException& error) {
260 ERROR(error << std::endl);
261 }
262 }
263 }
264
265 return h1;
266 }
267
268 static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
269 static double Vmax_npe; //!< Maximal integral of PDF [npe]
270
271 static const int NUMBER_OF_PDFS = 4;
272
273 static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
274
275 JNPE_t npe[NUMBER_OF_PDFS-2]; //!< PDF
276 JNPE_t2 npe2[NUMBER_OF_PDFS-2]; //!< PDF
277
278 JLANG::JSharedPointer<JMEstimator> estimator; //!< M-Estimator function
279 };
280
281 /**
282 * PDF types.
283 */
284 const JPDFType_t JRegressor<JShowerEH, JSimplex>::pdf_t[] = { DIRECT_LIGHT_FROM_EMSHOWER,
285 SCATTERED_LIGHT_FROM_EMSHOWER,
286 DIRECT_LIGHT_FROM_BRIGHT_POINT,
287 SCATTERED_LIGHT_FROM_BRIGHT_POINT };
288
289 /**
290 * Default values.
291 */
293 double JRegressor<JShowerEH, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
294
295}
296
297#endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
Various implementations of functional maps.
Maximum likelihood estimator (M-estimators).
General purpose messaging.
#define ERROR(A)
Definition JMessage.hh:66
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
Auxiliary methods for PDF calculations.
Numbering scheme for PDF types.
Physics constants.
Auxiliary class to define a range between two values.
General purpose data regression method.
This include file containes various data structures that can be used as specific return types for the...
Definition of zero value for any class.
const JVersor3Z & getDirection() const
Get direction.
Definition JVersor3Z.hh:81
Data structure for fit of straight line in positive z-direction with energy.
Definition JShowerEH.hh:32
double getEem() const
Get EM energy.
Definition JShowerEH.hh:209
double getEh() const
Get Hadronic energy.
Definition JShowerEH.hh:249
double getBy() const
Get bjorken y.
Definition JShowerEH.hh:107
Data structure for direction in three dimensions.
const JDirection3D & getDirection() const
Get direction.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
Rotation around Z-axis.
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
double getDY() const
Get y direction.
Definition JVersor3Z.hh:158
double getDX() const
Get x direction.
Definition JVersor3Z.hh:147
General exception.
Definition JException.hh:24
virtual const char * what() const override
Get error message.
Definition JException.hh:64
The template JSharedPointer class can be used to share a pointer to an object.
Custom class for integrated values of the PDF of the arrival time of Cherenkov light.
Definition JNPETable.hh:46
Multi-dimensional PDF table for arrival time of Cherenkov light.
Definition JPDFTable.hh:44
Auxiliary classes and methods for linear and iterative data regression.
Definition JEnergy.hh:15
double getChi2(const double P)
Get chi2 corresponding to given probability.
static const JZero zero
Function object to assign zero value.
Definition JZero.hh:105
JPDFType_t
PDF types.
Definition JPDFTypes.hh:24
@ SCATTERED_LIGHT_FROM_EMSHOWER
scattered light from EM shower
Definition JPDFTypes.hh:38
@ SCATTERED_LIGHT_FROM_BRIGHT_POINT
scattered light from bright point
Definition JPDFTypes.hh:43
@ DIRECT_LIGHT_FROM_BRIGHT_POINT
direct light from bright point
Definition JPDFTypes.hh:42
@ DIRECT_LIGHT_FROM_EMSHOWER
direct light from EM shower
Definition JPDFTypes.hh:37
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition JResult.hh:998
Auxiliary class for handling PMT geometry, rate and response.
Definition JPMTW0.hh:24
int getN() const
Get number of hits.
Definition JPMTW0.hh:67
double getR() const
Get rate.
Definition JPMTW0.hh:56
Template definition of a data regressor of given model.
Definition JRegressor.hh:70
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Map list.
Definition JMapList.hh:25
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
Type definition of a zero degree polynomial interpolation based on a JMap implementation.
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type.