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)) {
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;
122 norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
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++){