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  AbsFunction(right),
17  _lifetime (right._lifetime),
18  _sigma (right._sigma),
19  _punctures (right._punctures)
20 {
21 }
22 
23 void 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 
36 const Parameter & PuncturedSmearedExp::min(unsigned int i) const {
37  return _punctures[2*i];
38 }
39 
40 
42  return _punctures[2*i+1];
43 
44 }
45 
46 const 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 
71 double 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 
135 double 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 
148 double 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