10 _lifetime ("Lifetime", 1.0, 0.0),
11 _sigma ("
Sigma", 1.0, 0.0)
16 _lifetime (right._lifetime),
17 _sigma (right._sigma),
18 _punctures (right._punctures)
23 std::ostringstream mn, mx;
24 mn <<
"Min_" << _punctures.size()/2;
25 mx <<
"Max_" << _punctures.size()/2;
26 _punctures.push_back(
Parameter(mn.str(),
min, 0.0, 10.0));
27 _punctures.push_back(
Parameter(mx.str(),
max, 0.0, 10.0));
32 return _punctures[2*i];
36 return _punctures[2*i];
41 return _punctures[2*i+1];
46 return _punctures[2*i+1];
72 static const double sqrtTwo = sqrt(2.0);
79 std::vector<double> punctures(_punctures.size());
80 for (
size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
89 for (
size_t i=0;i<punctures.size()/2;i++) {
90 std::sort(punctures.begin()+2*i, punctures.begin()+2*i+2);
91 double min1=punctures[2*i];
92 double max1=punctures[2*i+1];
93 for (
size_t j=i+1;j<punctures.size()/2;j++) {
94 std::sort(punctures.begin()+2*j, punctures.begin()+2*j+2);
95 double min2=punctures[2*j];
96 double max2=punctures[2*j+1];
97 if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) {
99 punctures[2*i+1]=
std::max(max1,max2);
100 std::vector<double>::iterator t0=punctures.begin()+2*j, t1=t0+2;
101 punctures.erase(t0,t1);
110 double expG=0,
norm=0;
111 for (
size_t i=0;i<punctures.size()/2;i++) {
113 double a = punctures[2*i];
114 double b = punctures[2*i+1];
116 double alpha = (a/sigma + sigma/tau)/sqrtTwo;
117 double beta = (b/sigma + sigma/tau)/sqrtTwo;
118 double delta = 1/sqrtTwo/
sigma;
121 norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
123 expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
134 double PuncturedSmearedExp::erfc(
double x)
const {
138 z = (x < 0) ? -x : x;
140 ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
141 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
142 t*(-0.82215223+t*0.17087277 ))) ))) )));
143 if ( x < 0 ) ans = 2.0 - ans;
147 double PuncturedSmearedExp::pow(
double x,
int n)
const {
149 for(
int i=0;i<
n;i++){