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>
6 namespace Genfun {
7 FUNCTION_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  _lifetime (right._lifetime),
17  _sigma (right._sigma),
18  _punctures (right._punctures)
19 {
20 }
21 
22 void PuncturedSmearedExp::puncture(double min, double max) {
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));
28 }
29 
30 
32  return _punctures[2*i];
33 }
34 
35 const Parameter & PuncturedSmearedExp::min(unsigned int i) const {
36  return _punctures[2*i];
37 }
38 
39 
41  return _punctures[2*i+1];
42 
43 }
44 
45 const Parameter & PuncturedSmearedExp::max(unsigned int i) const {
46  return _punctures[2*i+1];
47 }
48 
49 
51 }
52 
54  return _lifetime;
55 }
56 
58  return _lifetime;
59 }
60 
62  return _sigma;
63 }
64 
66  return _sigma;
67 }
68 
69 
70 double PuncturedSmearedExp::operator() (double argument) const {
71  // Fetch the paramters. This operator does not convolve numerically.
72  static const double sqrtTwo = sqrt(2.0);
73 
74  double sigma = _sigma.getValue();
75  double tau = _lifetime.getValue();
76  double x = argument;
77 
78  // Copy:
79  std::vector<double> punctures(_punctures.size());
80  for (size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
81 
82  // Overlap elimination:
83  bool overlap=true;
84 
85  while (overlap) {
86 
87  overlap=false;
88 
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)) {
98  punctures[2*i] =std::min(min1,min2);
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);
102  overlap=true;
103  break;
104  }
105  }
106  if (overlap) break;
107  }
108  }
109 
110  double expG=0,norm=0;
111  for (size_t i=0;i<punctures.size()/2;i++) {
112 
113  double a = punctures[2*i];
114  double b = punctures[2*i+1];
115 
116  double alpha = (a/sigma + sigma/tau)/sqrtTwo;
117  double beta = (b/sigma + sigma/tau)/sqrtTwo;
118  double delta = 1/sqrtTwo/sigma;
119 
120 
121  norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
122 
123  expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
124 
125 
126  }
127 
128  // protection:
129  if (norm==0) return norm;
130 
131  return expG/norm;
132 }
133 
134 double PuncturedSmearedExp::erfc(double x) const {
135  // This is accurate to 7 places.
136  // See Numerical Recipes P 221
137  double t, z, ans;
138  z = (x < 0) ? -x : x;
139  t = 1.0/(1.0+.5*z);
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;
144  return ans;
145 }
146 
147 double PuncturedSmearedExp::pow(double x,int n) const {
148  double val=1;
149  for(int i=0;i<n;i++){
150  val=val*x;
151  }
152  return val;
153 }
154 
155 } // namespace Genfun
156 
157