10 _lifetime ("Lifetime", 1.0, 0.0),
11 _sigma ("
Sigma", 1.0, 0.0)
17 _lifetime (right._lifetime),
18 _sigma (right._sigma),
19 _punctures (right._punctures)
24 std::ostringstream mn, mx;
25 mn <<
"Min_" << _punctures.size()/2;
26 mx <<
"Max_" << _punctures.size()/2;
27 _punctures.push_back(
Parameter(mn.str(), xmin, 0.0, 10.0));
28 _punctures.push_back(
Parameter(mx.str(), xmax, 0.0, 10.0));
33 return _punctures[2*
i];
37 return _punctures[2*
i];
42 return _punctures[2*
i+1];
47 return _punctures[2*
i+1];
73 static const double sqrtTwo = sqrt(2.0);
80 std::vector<double> punctures(_punctures.size());
81 for (
size_t i=0;
i<_punctures.size();
i++) punctures[
i]=_punctures[
i].getValue();
90 for (
size_t i=0;
i<punctures.size()/2;
i++) {
91 std::sort(punctures.begin()+2*
i, punctures.begin()+2*
i+2);
92 double min1=punctures[2*
i];
93 double max1=punctures[2*
i+1];
94 for (
size_t j=
i+1;
j<punctures.size()/2;
j++) {
95 std::sort(punctures.begin()+2*
j, punctures.begin()+2*
j+2);
96 double min2=punctures[2*
j];
97 double max2=punctures[2*
j+1];
98 if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) {
99 punctures[2*
i] =std::min(min1,min2);
100 punctures[2*
i+1]=std::max(max1,max2);
101 std::vector<double>::iterator t0=punctures.begin()+2*
j, t1=t0+2;
102 punctures.erase(t0,t1);
111 double expG=0,
norm=0;
112 for (
size_t i=0;
i<punctures.size()/2;
i++) {
114 double a = punctures[2*
i];
115 double b = punctures[2*
i+1];
117 double alpha = (
a/xsigma + xsigma/tau)/sqrtTwo;
118 double beta = (
b/xsigma + xsigma/tau)/sqrtTwo;
119 double delta = 1/sqrtTwo/xsigma;
124 expG += ((erfc(alpha-
delta*
x) - erfc(beta-
delta*
x))*exp(-
x/tau) );
135 double PuncturedSmearedExp::erfc(
double x)
const {
139 z = (
x < 0) ? -
x :
x;
141 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
142 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
143 t*(-0.82215223+t*0.17087277 ))) ))) )));
144 if (
x < 0 ) ans = 2.0 - ans;
148 double PuncturedSmearedExp::pow(
double x,
int n)
const {
150 for(
int i=0;
i<
n;
i++){