Jpp  17.1.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMath/JMathSupportkit.hh
Go to the documentation of this file.
1 #ifndef __JMATHSUPPORTKIT__
2 #define __JMATHSUPPORTKIT__
3 
4 #include <limits>
5 #include <cmath>
6 
7 #include "JMath/JConstants.hh"
8 #include "JLang/JException.hh"
9 
10 
11 /**
12  * \file
13  *
14  * Auxiliary methods for mathematics.
15  * \author mdejong
16  */
17 
18 namespace JMATH {}
19 namespace JPP { using namespace JMATH; }
20 
21 namespace JMATH {
22 
24 
25 
26  /**
27  * Gauss function (normalised to 1 at x = 0).
28  *
29  * \param x x
30  * \param sigma sigma
31  * \return function value
32  */
33  inline double gauss(const double x, const double sigma)
34  {
35  const double u = x / sigma;
36 
37  if (fabs(u) < 20.0)
38  return exp(-0.5*u*u);
39  else
40  return 0.0;
41  }
42 
43 
44  /**
45  * Gauss function (normalised to 1 at x = x0).
46  *
47  * \param x x
48  * \param x0 central value
49  * \param sigma sigma
50  * \return function value
51  */
52  inline double gauss(const double x, const double x0, const double sigma)
53  {
54  return gauss(x - x0, sigma);
55  }
56 
57 
58  /**
59  * Normalised Gauss function.
60  *
61  * \param x x
62  * \param sigma sigma
63  * \return function value
64  */
65  inline double Gauss(const double x, const double sigma)
66  {
67  return gauss(x, sigma) / sqrt(2.0*PI) / sigma;
68  }
69 
70 
71  /**
72  * Normalised Gauss function.
73  *
74  * \param x x
75  * \param x0 central value
76  * \param sigma sigma
77  * \return function value
78  */
79  inline double Gauss(const double x, const double x0, const double sigma)
80  {
81  return Gauss(x - x0, sigma);
82  }
83 
84 
85  /**
86  * Incomplete gamma function.
87  *
88  * Source code is taken from reference:
89  * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
90  * Cambridge University Press.
91  *
92  * \param a a
93  * \param x x
94  * \return function value
95  */
96  inline double Gamma(const double a, const double x)
97  {
98  using namespace std;
99 
100  const int max = 1000000;
101 
102  if (x < 0.0) { THROW(JValueOutOfRange, "x < 0 " << x); }
103  if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
104 
105  const double gln = lgamma(a);
106 
107  if (x < a + 1.0) {
108 
109  if (x < 0.0) {
110  THROW(JValueOutOfRange, "x <= 0 " << x);
111  }
112 
113  if (x == 0.0) {
114  return 0.0;
115  }
116 
117  double ap = a;
118  double sum = 1.0 / a;
119  double del = sum;
120 
121  for (int i = 0; i != max; ++i) {
122 
123  ap += 1.0;
124  del *= x/ap;
125  sum += del;
126 
127  if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
128  return sum*exp(-x + a*log(x) - gln);
129  }
130  }
131 
132  THROW(JValueOutOfRange, "i == " << max);
133 
134  } else {
135 
136  const double FPMIN = numeric_limits<double>::min() / numeric_limits<double>::epsilon();
137 
138  double b = x + 1.0 - a;
139  double c = 1.0 / FPMIN;
140  double d = 1.0 / b;
141  double h = d;
142 
143  for (int i = 1; i != max; ++i) {
144 
145  const double an = -i * (i-a);
146 
147  b += 2.0;
148  d = an*d + b;
149 
150  if (fabs(d) < FPMIN) {
151  d = FPMIN;
152  }
153 
154  c = b + an/c;
155 
156  if (fabs(c) < FPMIN) {
157  c = FPMIN;
158  }
159 
160  d = 1.0/d;
161 
162  const double del = d*c;
163 
164  h *= del;
165 
166  if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
167  return 1.0 - exp(-x + a*log(x) - gln) * h;
168  }
169  }
170 
171  THROW(JValueOutOfRange, "i == " << max);
172  }
173  }
174 
175 
176  /**
177  * Legendre polynome.
178  *
179  * \param n degree
180  * \param x x
181  * \return function value
182  */
183  inline double legendre(const size_t n, const double x)
184  {
185  switch (n) {
186 
187  case 0:
188  return 1.0;
189 
190  case 1:
191  return x;
192 
193  default:
194  {
195  double p0;
196  double p1 = 1.0;
197  double p2 = x;
198 
199  for (size_t i = 2; i <= n; ++i) {
200  p0 = p1;
201  p1 = p2;
202  p2 = ((2*i-1) * x*p1 - (i-1) * p0) / i;
203  }
204 
205  return p2;
206  }
207  }
208  }
209 
210 
211  /**
212  * Binomial function.
213  *
214  * \param n n
215  * \param k k
216  * \return function value
217  */
218  inline double binomial(const size_t n, const size_t k)
219  {
220  if (n == 0 || n < k) {
221  return 0.0;
222  }
223 
224  if (k == 0 || n == k) {
225  return 1.0;
226  }
227 
228  const int k1 = std::min(k, n - k);
229  const int k2 = n - k1;
230 
231  double value = k2 + 1;
232 
233  for (int i = k1; i != 1; --i) {
234  value *= (double) (k2 + i) / (double) i;
235  }
236 
237  return value;
238  }
239 
240 
241  /**
242  * Poisson probability density distribition.
243  *
244  * \param n number of occurences
245  * \param mu expectation value
246  * \return probability
247  */
248  inline double poisson(const size_t n, const double mu)
249  {
250  using namespace std;
251 
252  if (mu > 0.0) {
253 
254  if (n > 0)
255  return exp(n*log(mu) - lgamma(n+1) - mu);
256  else
257  return exp(-mu);
258  } else if (mu == 0.0) {
259 
260  return (n == 0 ? 1.0 : 0.0);
261  }
262 
263  THROW(JValueOutOfRange, "mu <= 0 " << mu);
264  }
265 
266 
267  /**
268  * Poisson cumulative density distribition.
269  *
270  * \param n number of occurences
271  * \param mu expectation value
272  * \return probability
273  */
274  inline double Poisson(const size_t n, const double mu)
275  {
276  return 1.0 - Gamma(n + 1, mu);
277  }
278 }
279 
280 #endif
then cat $TRIPOD_INITIAL<< EOF1 256877.5 4743716.7-2438.42 256815.5 4743395.0-2435.53 257096.2 4743636.0-2439.5EOFfiJEditDetector-a $DETECTOR_INITIAL-s"-1 addz -6.9"-o $DETECTOReval`JPrintDetector-a $DETECTOR-O SUMMARY`for STRING in ${STRINGS[*]};do set_variable MODULE`getModule-a $DETECTOR-L"$STRING 0"`JEditDetector-a $DETECTOR-M"$MODULE setz -2.9"-o $DETECTORdonecp-p $TRIPOD_INITIAL $TRIPODJAcoustics.sh $DETECTOR_IDcat > acoustics_trigger_parameters txt<< EOFQ=0.0;TMax_s=0.020;numberOfHits=90;EOFJAcousticsEventBuilder.sh $DETECTOR $RUNS[*]INPUT_FILES=(`ls KM3NeT_ ${(l:8::0::0:) DETECTOR_ID}_0 *${^RUNS}_event.root`) cd $WORKDIRif[!$HOMEDIR-ef $WORKDIR];then cp-p $HOMEDIR/$DETECTOR $WORKDIR cp-p $HOMEDIR/$TRIPOD $WORKDIR cp-p $HOMEDIR/${^INPUT_FILES}$WORKDIR cp-p $HOMEDIR/{acoustics_fit_parameters, acoustics_trigger_parameters, disable, hydrophone, mechanics, sound_velocity, tripod, waveform}.txt $WORKDIRfisource $JPP_DIR/examples/JAcoustics/acoustics-fit-toolkit.shtimer_startinitialise stage_1B > &stage log
Exceptions.
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
double Gamma(const double a, const double x)
Incomplete gamma function.
const double sigma[]
Definition: JQuadrature.cc:74
const int n
Definition: JPolint.hh:676
Mathematical constants.
static const double PI
Mathematical constants.
p2
Definition: module-Z:fit.sh:74
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
double Poisson(const size_t n, const double mu)
Poisson cumulative density distribition.
then JCalibrateToT a
Definition: JTuneHV.sh:116
then set_variable NUMBER_OF_TESTS else set_variable NUMBER_OF_TESTS fi function gauss()
double legendre(const size_t n, const double x)
Legendre polynome.
double poisson(const size_t n, const double mu)
Poisson probability density distribition.
double Gauss(const double x, const double sigma)
Normalised Gauss function.
$WORKDIR ev_configure_domsimulator txt echo process $DOM_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DOM_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
double u[N+1]
Definition: JPolint.hh:755
double binomial(const size_t n, const size_t k)
Binomial function.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` source JAcousticsToolkit.sh typeset -A TRIPODS get_tripods $WORKDIR/tripod.txt TRIPODS XMEAN
const double epsilon
Definition: JQuadrature.cc:21