Jpp  16.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JShowerBjorkenYRegressor.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JSHOWERBJORKENYREGRESSOR__
2 #define __JFIT__JSHOWERBJORKENYREGRESSOR__
3 
4 #include "JPhysics/JPDFTypes.hh"
5 #include "JPhysics/JPDFTable.hh"
7 #include "JPhysics/JNPETable.hh"
8 
11 #include "JPhysics/JConstants.hh"
12 #include "JTools/JRange.hh"
13 #include "JTools/JResult.hh"
14 #include "JTools/JFunction1D_t.hh"
17 
18 #include "JGeometry3D/JVector3D.hh"
19 #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/JMEstimator.hh"
28 #include "JFit/JRegressor.hh"
29 #include "JFit/JShowerEH.hh"
30 #include "JFit/JFitToolkit.hh"
31 
32 #include "Jeep/JMessage.hh"
33 
34 /**
35  * \file
36  * Data regression method for JFIT::JShowerEH.
37  * \author adomi
38  */
39 
40 namespace JFIT {
41 
47  using JTOOLS::get_value;
48 
49  /**
50  * Regressor function object for JShowerEH fit using JSimplex minimiser.
51  */
52  template<>
54  public JAbstractRegressor<JShowerEH, JSimplex>
55  {
57 
64 
70 
71 
72  // Part for the Isotropic PDF
76 
80 
81 
82  /**
83  * Parameterized constructor
84  *
85  * The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD which
86  * will be replaced by the corresponding PDF types listed in JRegressor<JShower3Z, JGandalf>::pdf_t.
87  *
88  * \param fileDescriptor PDF file descriptor
89  */
90 
91  JRegressor(const std::string& fileDescriptor):
92  estimator(new JMEstimatorNull())
93  {
94  using namespace std;
95  using namespace JPP;
96 
97  const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(JMATH::zero));
98  const JPDF_t2::JSupervisor supervisor2(new JPDF_t2::JDefaultResult(JMATH::zero));
99 
100  for (int i = 0; i != NUMBER_OF_PDFS; ++i) {
101 
102  try {
103 
104  JPDF_t pdf;
105  JPDF_t2 pdf2;
106 
107  const string file_name = getFilename(fileDescriptor, pdf_t[i]);
108 
109  NOTICE("loading PDF from file " << file_name << "... " << flush);
110 
111  if(i < 2){
112 
113  pdf.load(file_name.c_str());
114 
115  pdf.setExceptionHandler(supervisor);
116 
117  npe[ i ] = JNPE_t(pdf);
118 
119  } else {
120 
121  pdf2.load(file_name.c_str());
122 
123  NOTICE("OK" << endl);
124 
125  pdf2.setExceptionHandler(supervisor2);
126 
127  npe2[ i-2 ] = JNPE_t2(pdf2);
128 
129  }
130  }
131  catch(const JException& error) {
132  FATAL(error.what() << endl);
133  }
134  }
135 
136  // Add PDFs
137  for (int i = 1; i < (NUMBER_OF_PDFS-2); i += 2) {
138 
139  npe[ i ].add(npe[i-1]);
140 
141  JNPE_t buffer;
142 
143  npe[i-1].swap(buffer);
144 
145  npe2[ i ].add(npe2[i-1]);
146 
147  JNPE_t2 buffer2;
148 
149  npe2[i-1].swap(buffer2);
150 
151  }
152  }
153 
154  /**
155  * Fit function.
156  * This method is used to determine the chi2 of given PMT with respect to shower hypothesis.
157  *
158  * \param shower shower
159  * \param pmt pmt
160  * \return chi2
161  */
162  double operator()(const JShowerEH& shower, const JPMTW0& pmt) const
163  {
164  using namespace JPP;
165  using namespace std;
166 
167  JPosition3D D(pmt.getPosition());
168  JDirection3D U(pmt.getDirection());
169 
170  D.sub(shower.getPosition());
171 
172  double ct = U.getDot(D) / D.getLength();
173 
174  JVersor3D shower_dir(shower.getDirection());
175 
176  const double z = D.getDot(shower_dir);
177  const double x = D.getX() - z * shower.getDX();
178  const double y = D.getY() - z * shower.getDY();
179  const double cosDelta = z/D.getLength(); // Delta = angle between shower direction and PMT position
180 
181  U.rotate(JRotation3Z(-atan2(y,x))); // rotate PMT axis to x-z plane
182 
183  const double theta = U.getTheta();
184  const double phi = fabs(U.getPhi());
185 
186  double H0 = getH0(pmt.getR()); // background
187  double H1 = getH1(D.getLength(), ct, cosDelta, theta, phi,
188  shower.getEem(), shower.getEh(), shower.getBy()); // signal
189 
190  if (H1 >= Vmax_npe) {
191  H1 *= Vmax_npe / H1;
192  }
193 
194  H1 += H0; // now H1 is signal + background
195 
196  const bool hit = pmt.getN() != 0;
197  const double u = getChi2(H1, hit); // -log(lik)
198 
199  return estimator->getRho(u);
200  }
201 
202  /**
203  * Get background hypothesis value for time integrated PDF.
204  *
205  * \param R_Hz rate [Hz]
206  * \return hypothesis value
207  */
208  double getH0(const double R_Hz) const
209  {
210  return get_value(JNPE_t::result_type(R_Hz * 1e-9 * T_ns.getLength()));
211  }
212 
213  /**
214  * Get signal hypothesis value for time integrated PDF.
215  *
216  * \param D PMT distance from shower [m]
217  * \param ct angle between shower direction and PMT position
218  * \param cosDelta angle between shower direction and PMT position
219  * \param theta PMT zenith angle [deg]
220  * \param phi PMT azimuth angle [deg]
221  * \param Eem EM shower energy [GeV]
222  * \param Eh H shower energy [GeV]
223  * \param Y Bjorken Y
224  * \return hypothesis value
225  */
226  double getH1(const double D,
227  const double ct,
228  const double cosDelta,
229  const double theta,
230  const double phi,
231  const double Eem,
232  const double Eh,
233  const double Y) const
234  {
235 
236  double h1 = 0;
237 
238  for (int i = 0; i != (NUMBER_OF_PDFS-1); ++i) {
239 
240  if (!npe[i].empty() && D <= npe[i].getXmax() && !npe2[i].empty() && D <= npe2[i].getXmax()) {
241 
242  try {
243 
244  JNPE_t::result_type P_em;
246 
247  P_em = abs(Eem) * npe[i](std::max(D, npe[i].getXmin()), cosDelta, theta, phi);
248 
249  P_h = abs(Eh) * npe2[i](std::max(D, npe2[i].getXmin()), ct);
250 
251  double y1 = get_value(P_em) + get_value(P_h);
252 
253  if(y1 > 0.0){
254  h1 += y1;
255  }
256 
257  }
258  catch(JLANG::JException& error) {
259  ERROR(error << std::endl);
260  }
261  }
262  }
263 
264  return h1;
265  }
266 
267  static JTimeRange T_ns; //!< Time window with respect to Cherenkov hypothesis [ns]
268  static double Vmax_npe; //!< Maximal integral of PDF [npe]
269 
270  static const int NUMBER_OF_PDFS = 4;
271 
272  static const JPDFType_t pdf_t[NUMBER_OF_PDFS];
273 
274  JNPE_t npe[NUMBER_OF_PDFS-2]; //!< PDF
275  JNPE_t2 npe2[NUMBER_OF_PDFS-2]; //!< PDF
276 
278  };
279 
280  /**
281  * PDF types.
282  */
287 
288  /**
289  * Default values.
290  */
292  double JRegressor<JShowerEH, JSimplex>::Vmax_npe = std::numeric_limits<double>::max();
293 
294 }
295 
296 #endif
Auxiliary methods to evaluate Poisson probabilities and chi2.
Template definition of a data regressor of given model.
Definition: JRegressor.hh:68
Type definition of a spline interpolation method based on a JCollection with JResultPDF result type...
General exception.
Definition: JException.hh:23
JRegressor(const std::string &fileDescriptor)
Parameterized constructor.
void load(const char *file_name)
Load from input file.
General purpose data regression method.
double getR() const
Get rate.
Definition: JPMTW0.hh:56
Null M-estimator.
Definition: JMEstimator.hh:53
Auxiliary methods for PDF calculations.
This include file containes various data structures that can be used as specific return types for the...
Auxiliary class for handling PMT geometry, rate and response.
Definition: JPMTW0.hh:22
scattered light from EM shower
Definition: JPDFTypes.hh:41
const JDirection3D & getDirection() const
Get direction.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
JLANG::JSharedPointer< JMEstimator > estimator
M-Estimator function.
double getEem() const
Get EM energy.
Definition: JShowerEH.hh:209
JPHYSICS::JPDFTable< JFunction1D_t, JMapList_t2 > JPDF_t2
double getH1(const double D, const double ct, const double cosDelta, const double theta, const double phi, const double Eem, const double Eh, const double Y) const
Get signal hypothesis value for time integrated PDF.
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
direct light from bright point
Definition: JPDFTypes.hh:45
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
Definition of zero value for any class.
double getH0(const double R_Hz) const
Get background hypothesis value for time integrated PDF.
double operator()(const JShowerEH &shower, const JPMTW0 &pmt) const
Fit function.
double getDY() const
Get y direction.
Definition: JVersor3Z.hh:158
Various implementations of functional maps.
Numbering scheme for PDF types.
void add(const JNPETable &input)
Add NPE table.
Definition: JNPETable.hh:124
Type definition of a 1st degree polynomial interpolation based on a JGridMap implementation.
The template JSharedPointer class can be used to share a pointer to an object.
JPHYSICS::JNPETable< double, double, JNPEMapList_t2 > JNPE_t2
Map list.
Definition: JMapList.hh:24
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
static JTimeRange T_ns
Time window with respect to Cherenkov hypothesis [ns].
#define NOTICE(A)
Definition: JMessage.hh:64
static double Vmax_npe
Maximal integral of PDF [npe].
Physics constants.
direct light from EM shower
Definition: JPDFTypes.hh:40
JTOOLS::JMapList< JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap > > > > JPDFMaplist_t
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JNPEMapList_t2
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
JPDFType_t
PDF types.
Definition: JPDFTypes.hh:27
Auxiliary class for recursive map list generation.
Definition: JMapList.hh:108
General purpose messaging.
scattered light from bright point
Definition: JPDFTypes.hh:46
const JVersor3Z & getDirection() const
Get direction.
Definition: JVersor3Z.hh:81
#define FATAL(A)
Definition: JMessage.hh:67
Data structure for fit of straight line in positive z-direction with energy.
Definition: JShowerEH.hh:28
void setExceptionHandler(const typename function_type::supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
virtual const char * what() const override
Get error message.
Definition: JException.hh:48
JPHYSICS::JNPETable< double, double, JNPEMaplist_t > JNPE_t
double getEh() const
Get Hadronic energy.
Definition: JShowerEH.hh:249
double getBy() const
Get bjorken y.
Definition: JShowerEH.hh:107
int getN() const
Get number of hits.
Definition: JPMTW0.hh:67
Auxiliary class to define a range between two values.
Simple fit method based on Powell&#39;s algorithm, see reference: Numerical Recipes in C++...
Definition: JSimplex.hh:42
double getDX() const
Get x direction.
Definition: JVersor3Z.hh:147
Type definition of a zero degree polynomial interpolation based on a JGridMap implementation.
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
JPHYSICS::JPDFTable< JFunction1D_t, JPDFMaplist_t > JPDF_t
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
double u[N+1]
Definition: JPolint.hh:755
JTOOLS::JMAPLIST< JTOOLS::JPolint1FunctionalMap, JTOOLS::JPolint1FunctionalGridMap >::maplist JMapList_t2
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
Abstract class for global fit method.
Definition: JRegressor.hh:75
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).
Type definition of a zero degree polynomial interpolation based on a JMap implementation.
JTOOLS::JMapList< JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap, JTOOLS::JMapList< JTOOLS::JPolint0FunctionalGridMap > > > > JNPEMaplist_t