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

AnalyticConvolution.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: AnalyticConvolution.cc,v 1.8 2010/07/22 21:55:10 garren Exp $
6 #include <cmath> // for isfinite
7 #if (defined _WIN32)
8 #include <float.h> // Visual C++ _finite
9 #endif
10 namespace Genfun {
11 FUNCTION_OBJECT_IMP(AnalyticConvolution)
12 
14  _lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default
15  _frequency("Frequency", 0.0, 0.0), // Bounded from below by zero, by default
16  _sigma ("Sigma", 1.0, 0.0), // Bounded from below by zero, by default
17  _offset ("Offset", 0.0), // Offset is unbounded
18  _type(type)
19 {
20 }
21 
23  AbsFunction(right),
24  _lifetime (right._lifetime),
25  _frequency (right._frequency),
26  _sigma (right._sigma),
27  _offset (right._offset),
28  _type(right._type)
29 {
30 }
31 
33 }
34 
36  return _frequency;
37 }
38 
40  return _frequency;
41 }
42 
44  return _lifetime;
45 }
46 
48  return _lifetime;
49 }
50 
52  return _sigma;
53 }
54 
55 const Parameter & AnalyticConvolution::sigma() const {
56  return _sigma;
57 }
58 
60  return _offset;
61 }
62 
63 const Parameter & AnalyticConvolution::offset() const {
64  return _offset;
65 }
66 double AnalyticConvolution::operator() (double argument) const {
67  // Fetch the paramters. This operator does not convolve numerically.
68  static const double sqrtTwo = sqrt(2.0);
69  double xsigma = _sigma.getValue();
70  double tau = _lifetime.getValue();
71  double xoffset = _offset.getValue();
72  double x = argument-xoffset;
73  double freq = _frequency.getValue();
74 
75  // smeared exponential an its asymmetry.
76  double expG=0.0, asymm=0.0;
77 
78  if (_type==SMEARED_NEG_EXP) {
79  expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/x))/(2.0*tau*tau)) *
80  erfc((xsigma*xsigma+tau*(/*xoffset*/x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
81 #if (defined _WIN32)
82  if (!_finite(expG)) {
83  expG=0.0;
84  }
85 #else
86  if (!std::isfinite(expG)) {
87  expG=0.0;
88  }
89 #endif
90  return expG;
91  }
92  else {
93  expG = exp((xsigma*xsigma +2*tau*(/*xoffset*/-x))/(2.0*tau*tau)) *
94  erfc((xsigma*xsigma+tau*(/*xoffset*/-x))/(sqrtTwo*xsigma*tau))/(2.0*tau);
95  }
96 
97  // Both sign distribution=> return smeared exponential:
98  if (_type==SMEARED_EXP) {
99 #if (defined _WIN32)
100  if (!_finite(expG)) {
101  expG=0.0;
102  }
103 #else
104  if (!std::isfinite(expG)) {
105  expG=0.0;
106  }
107 #endif
108  return expG;
109  }
110 
111 
112  // Calcualtion of asymmetry:
113 
114  // First, if the mixing frequency is too high compared with the lifetime, we
115  // cannot see this oscillation. We abandon the complicated approach and just do
116  // this instead:
117  if (xsigma>6.0*tau) {
118  asymm = expG*(1/(1+tau*tau*freq*freq));
119  }
120  else if (xsigma==0.0) {
121  if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) {
122  if (x>=0) asymm= (expG*cos(freq*x));
123  }
124  else if (_type==SMEARED_SIN_EXP) {
125  if (x>=0) asymm= (expG*sin(freq*x));
126  }
127  }
128  else {
129  std::complex<double> z(freq*xsigma/sqrtTwo, (xsigma/tau-x/xsigma)/sqrtTwo);
130  if (x<0) {
131  if (_type==SMEARED_COS_EXP|| _type==MIXED || _type==UNMIXED ) {
132  asymm= 2.0*nwwerf(z).real()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
133  }
134  else if (_type==SMEARED_SIN_EXP) {
135  asymm= 2.0*nwwerf(z).imag()/tau/4.0*exp(-x*x/2.0/xsigma/xsigma);
136  }
137  }
138  else {
139  if (_type==SMEARED_COS_EXP||_type==MIXED || _type==UNMIXED) {
140  asymm= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
141  exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*cos(freq*x - freq/tau*xsigma*xsigma);
142  }
143  else if (_type==SMEARED_SIN_EXP) {
144  asymm= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-x*x/2.0/xsigma/xsigma) +
145  exp(xsigma*xsigma/2 *(1/tau/tau - freq*freq) - x/tau)*(1./tau)*sin(freq*x - freq/tau*xsigma*xsigma);
146  }
147  }
148  }
149 
150  // Return either MIXED, UNMIXED, or ASYMMETRY function.
151  if (_type==UNMIXED) {
152  double retVal = (expG+asymm)/2.0;
153  if (retVal<0)
154  std::cerr
155  << "Warning in AnalyticConvolution: negative probablity"
156  << std::endl;
157  if (retVal<0)
158  std::cerr
159  << xsigma << ' ' << tau << ' ' << xoffset << ' '
160  << freq << ' '<< argument
161  << std::endl;
162  if (retVal<0)
163  std::cerr << retVal << std::endl;
164  return retVal;
165  }
166  else if (_type==MIXED) {
167  double retVal = (expG-asymm)/2.0;
168  if (retVal<0)
169  std::cerr
170  << "Warning in AnalyticConvolution: negative probablity"
171  << std::endl;
172  if (retVal<0)
173  std::cerr
174  << xsigma << ' ' << tau << ' ' << xoffset << ' '
175  << freq << ' ' << argument
176  << std::endl;
177  if (retVal<0)
178  std::cerr << retVal << std::endl;
179  return retVal;
180  }
181  else if (_type==SMEARED_COS_EXP || _type==SMEARED_SIN_EXP) {
182  return asymm;
183  }
184  else {
185  std::cerr
186  << "Unknown sign parity. State is not allowed"
187  << std::endl;
188  exit(0);
189  return 0.0;
190  }
191 
192 }
193 
194 double AnalyticConvolution::erfc(double x) const {
195 // This is accurate to 7 places.
196 // See Numerical Recipes P 221
197  double t, z, ans;
198  z = (x < 0) ? -x : x;
199  t = 1.0/(1.0+.5*z);
200  ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
201  t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
202  t*(-0.82215223+t*0.17087277 ))) ))) )));
203  if ( x < 0 ) ans = 2.0 - ans;
204  return ans;
205 }
206 
207 double AnalyticConvolution::pow(double x,int n) const {
208  double val=1;
209  for(int i=0;i<n;i++){
210  val=val*x;
211  }
212  return val;
213 }
214 
215 std::complex<double> AnalyticConvolution::nwwerf(std::complex<double> z) const {
216  std::complex<double> zh,r[38],s,t,v;
217 
218  const double z1 = 1;
219  const double hf = z1/2;
220  const double z10 = 10;
221  const double c1 = 74/z10;
222  const double c2 = 83/z10;
223  const double c3 = z10/32;
224  const double c4 = 16/z10;
225  const double c = 1.12837916709551257;
226  const double p = pow(2.0*c4,33);
227 
228  double x=z.real();
229  double y=z.imag();
230  double xa=(x >= 0) ? x : -x;
231  double ya=(y >= 0) ? y : -y;
232  if(ya < c1 && xa < c2){
233  zh = std::complex<double>(ya+c4,xa);
234  r[37]= std::complex<double>(0,0);
235  // do 1 n = 36,1,-1
236  for(int n = 36; n>0; n--){
237  t=zh+double(n)*std::conj(r[n+1]);
238  r[n]=hf*t/std::norm(t);
239  }
240  double xl=p;
241  s=std::complex<double>(0,0);
242  // do 2 n = 33,1,-1
243  for(int k=33; k>0; k--){
244  xl=c3*xl;
245  s=r[k]*(s+xl);
246  }
247  v=c*s;
248  }
249  else{
250  zh=std::complex<double>(ya,xa);
251  r[1]=std::complex<double>(0,0);
252  // do 3 n = 9,1,-1
253  for(int n=9;n>0;n--){
254  t=zh+double(n)*std::conj(r[1]);
255  r[1]=hf*t/std::norm(t);
256  }
257  v=c*r[1];
258  }
259  if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
260  if(y < 0) {
261  v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
262  if(x > 0) v=std::conj(v);
263  }
264  else{
265  if(x < 0) v=std::conj(v);
266  }
267  return v;
268 }
269 } // namespace Genfun
Genfun::Parameter::getValue
virtual double getValue() const
Definition: Parameter.cc:27
double
#define double(obj)
Definition: excDblThrow.cc:32
Genfun::AnalyticConvolution::UNMIXED
@ UNMIXED
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:35
Genfun::AbsFunction
Definition: CLHEP/GenericFunctions/AbsFunction.hh:48
Genfun::AnalyticConvolution::offset
Parameter & offset()
Definition: AnalyticConvolution.cc:59
AnalyticConvolution.hh
Genfun::AnalyticConvolution
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:27
Genfun::AnalyticConvolution::SMEARED_NEG_EXP
@ SMEARED_NEG_EXP
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:39
Genfun::AnalyticConvolution::SMEARED_COS_EXP
@ SMEARED_COS_EXP
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:37
Genfun::AnalyticConvolution::SMEARED_SIN_EXP
@ SMEARED_SIN_EXP
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:38
Genfun::AnalyticConvolution::SMEARED_EXP
@ SMEARED_EXP
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:36
CLHEP::detail::n
n
Definition: Ranlux64Engine.cc:85
Genfun::AnalyticConvolution::operator()
virtual double operator()(double argument) const
Definition: AnalyticConvolution.cc:66
type
Signatures of Hep3Vector::rotate For equivalent ZOOM axis There is no harm in leaving this axis CLHEP has implemented a first forming an identity then rotating that by axis and I leave the CLHEP code alone people are of course free to use the ZOOM originated method with signature which I believe will be faster Return types for rotateZ CLHEP and PhysicsVectors each have these three and they are identical except that the ZOOM version returns a reference to while in CLHEP they return void Having methods that alter an object return a reference to that object is convenient for certain chained and costs nothing I don t wish to potentially break ZOOM user code for no good so I have made these CLHEP method conform to this convention There are a couple of other CLHEP rotate and which use the void return type
Definition: minorMergeIssues.doc:113
Genfun::AnalyticConvolution::lifetime
Parameter & lifetime()
Definition: AnalyticConvolution.cc:43
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
Genfun::AnalyticConvolution::frequency
Parameter & frequency()
Definition: AnalyticConvolution.cc:35
Genfun::AnalyticConvolution::MIXED
@ MIXED
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:34
Genfun::AnalyticConvolution::AnalyticConvolution
AnalyticConvolution(Type=SMEARED_EXP)
Definition: AnalyticConvolution.cc:13
Genfun::AnalyticConvolution::~AnalyticConvolution
virtual ~AnalyticConvolution()
Definition: AnalyticConvolution.cc:32
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
Genfun::Sigma
Definition: CLHEP/GenericFunctions/Sigma.hh:19
i
long i
Definition: JamesRandomSeeding.txt:27
Gaussian.hh
Genfun::AnalyticConvolution::sigma
Parameter & sigma()
Definition: AnalyticConvolution.cc:51
Exponential.hh
Genfun::Parameter
Definition: CLHEP/GenericFunctions/Parameter.hh:35
exit
it has advantages For I leave the ZMthrows but substitute I replaced ZMthrow with ZMthrowA in this package when will issue its message and exit(-1). Also
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
k
long k
Definition: JamesRandomSeeding.txt:29
Genfun::AnalyticConvolution::Type
Type
Definition: CLHEP/GenericFunctions/AnalyticConvolution.hh:34
FUNCTION_OBJECT_IMP
#define FUNCTION_OBJECT_IMP(classname)
Definition: CLHEP/GenericFunctions/AbsFunction.hh:156
CLHEP::norm
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57
Genfun
Definition: CLHEP/GenericFunctions/Abs.hh:14