CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

PuncturedSmearedExp.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id: PuncturedSmearedExp.cc,v 1.4 2010/06/16 18:22:01 garren Exp $
4#include <sstream>
5#include <cmath>
6namespace Genfun {
7FUNCTION_OBJECT_IMP(PuncturedSmearedExp)
8
10 _lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default
11 _sigma ("Sigma", 1.0, 0.0) // Bounded from below by zero, by default
12{
13}
14
16 AbsFunction(right),
17 _lifetime (right._lifetime),
18 _sigma (right._sigma),
19 _punctures (right._punctures)
20{
21}
22
23void PuncturedSmearedExp::puncture(double xmin, double xmax) {
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));
29}
30
31
33 return _punctures[2*i];
34}
35
36const Parameter & PuncturedSmearedExp::min(unsigned int i) const {
37 return _punctures[2*i];
38}
39
40
42 return _punctures[2*i+1];
43
44}
45
46const Parameter & PuncturedSmearedExp::max(unsigned int i) const {
47 return _punctures[2*i+1];
48}
49
50
52}
53
55 return _lifetime;
56}
57
59 return _lifetime;
60}
61
63 return _sigma;
64}
65
67 return _sigma;
68}
69
70
71double PuncturedSmearedExp::operator() (double argument) const {
72 // Fetch the paramters. This operator does not convolve numerically.
73 static const double sqrtTwo = sqrt(2.0);
74
75 double xsigma = _sigma.getValue();
76 double tau = _lifetime.getValue();
77 double x = argument;
78
79 // Copy:
80 std::vector<double> punctures(_punctures.size());
81 for (size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
82
83 // Overlap elimination:
84 bool overlap=true;
85
86 while (overlap) {
87
88 overlap=false;
89
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);
103 overlap=true;
104 break;
105 }
106 }
107 if (overlap) break;
108 }
109 }
110
111 double expG=0,norm=0;
112 for (size_t i=0;i<punctures.size()/2;i++) {
113
114 double a = punctures[2*i];
115 double b = punctures[2*i+1];
116
117 double alpha = (a/xsigma + xsigma/tau)/sqrtTwo;
118 double beta = (b/xsigma + xsigma/tau)/sqrtTwo;
119 double delta = 1/sqrtTwo/xsigma;
120
121
122 norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
123
124 expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
125
126
127 }
128
129 // protection:
130 if (norm==0) return norm;
131
132 return expG/norm;
133}
134
135double PuncturedSmearedExp::erfc(double x) const {
136 // This is accurate to 7 places.
137 // See Numerical Recipes P 221
138 double t, z, ans;
139 z = (x < 0) ? -x : x;
140 t = 1.0/(1.0+.5*z);
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;
145 return ans;
146}
147
148double PuncturedSmearedExp::pow(double x,int n) const {
149 double val=1;
150 for(int i=0;i<n;i++){
151 val=val*x;
152 }
153 return val;
154}
155
156} // namespace Genfun
157
158
#define FUNCTION_OBJECT_IMP(classname)
virtual double getValue() const
Definition: Parameter.cc:27
void puncture(double min, double max)
Parameter & min(unsigned int i)
virtual double operator()(double argument) const
Parameter & max(unsigned int i)
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57
@ b
@ a