Jpp
 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 "JTools/JConstants.hh"
8 #include "JLang/JException.hh"
9 
10 
11 namespace JMATH {}
12 namespace JPP { using namespace JMATH; }
13 
14 namespace JMATH {
15 
17 
18 
19  /**
20  * Gauss function (normalised to 1 at x = 0).
21  *
22  * \param x x
23  * \param sigma sigma
24  * \return function value
25  */
26  inline double gauss(const double x, const double sigma)
27  {
28  const double u = x / sigma;
29 
30  if (fabs(u) < 10.0)
31  return exp(-0.5*u*u);
32  else
33  return 0.0;
34  }
35 
36 
37  /**
38  * Gauss function (normalised to 1 at x = x0).
39  *
40  * \param x x
41  * \param x0 central value
42  * \param sigma sigma
43  * \return function value
44  */
45  inline double gauss(const double x, const double x0, const double sigma)
46  {
47  return gauss(x - x0, sigma);
48  }
49 
50 
51  /**
52  * Normalised Gauss function.
53  *
54  * \param x x
55  * \param sigma sigma
56  * \return function value
57  */
58  inline double Gauss(const double x, const double sigma)
59  {
60  return gauss(x, sigma) / sqrt(2.0*JTOOLS::PI) / sigma;
61  }
62 
63 
64  /**
65  * Normalised Gauss function.
66  *
67  * \param x x
68  * \param x0 central value
69  * \param sigma sigma
70  * \return function value
71  */
72  inline double Gauss(const double x, const double x0, const double sigma)
73  {
74  return Gauss(x - x0, sigma);
75  }
76 
77 
78  /**
79  * Incomplete gamma function.
80  *
81  * \param a a
82  * \param x x
83  * \return function value
84  */
85  inline double Gamma(const double a, const double x)
86  {
87  using namespace std;
88 
89  const int max = 100;
90 
91  if (x < 0.0) { THROW(JValueOutOfRange, "x < 0 " << x); }
92  if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
93 
94  const double gln = lgamma(a);
95 
96  if (x < a + 1.0) {
97 
98  if (x <= 0.0) {
99  THROW(JValueOutOfRange, "x <= 0 " << x);
100  }
101 
102  double ap = a;
103  double sum = 1.0 /a;
104  double del = sum;
105 
106  for (int i = 1; i != max; ++i) {
107 
108  ap += 1.0;
109  del *= x/ap;
110  sum += del;
111 
112  if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
113  return sum*exp(-x + a*log(x) - gln);
114  }
115  }
116 
117  } else {
118 
119  double b = x + 1.0 - a;
120  double c = numeric_limits<double>::epsilon() / numeric_limits<double>::min();
121  double d = 1.0 / b;
122  double h = d;
123 
124  for (int i = 1; i != max; ++i) {
125 
126  const double an = -i * (i-a);
127 
128  b += 2.0;
129  d = an*d + b;
130 
131  if (fabs(d) < numeric_limits<double>::min()) {
132  d = numeric_limits<double>::min();
133  }
134 
135  c = b + an/c;
136 
137  if (fabs(c) < numeric_limits<double>::min()) {
138  c = numeric_limits<double>::min();
139  }
140 
141  d = 1.0/d;
142 
143  const double del = d*c;
144 
145  h *= del;
146 
147  if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
148  return 1.0 - exp(-x + a*log(x) - gln) * h;
149  }
150  }
151 
152  THROW(JValueOutOfRange, "a " << a);
153  }
154 
155  return 0.0;
156  }
157 }
158 
159 #endif
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
static const double PI
Constants.
Definition: JConstants.hh:20
double Gamma(const double a, const double x)
Incomplete gamma function.
fi JEventTimesliceWriter a
Constants.
then print_variable DETECTOR INPUT_FILE INTERMEDIATE_FILE check_input_file $DETECTOR $INPUT_FILE check_output_file $INTERMEDIATE_FILE $OUTPUT_FILE JMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JPath.sh:52
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
double Gauss(const double x, const double sigma)
Normalised Gauss function.
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:706
then set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable OUTPUT_FILE histogram.root JHistogram1D -o $WORKDIR/$OUTPUT_FILE -F "$FORMULA" -