herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
MamboDecayer.h
1 // -*- C++ -*-
2 //
3 // MamboDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2019 The Herwig Collaboration
5 //
6 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef HERWIG_MamboDecayer_H
10 #define HERWIG_MamboDecayer_H
11 //
12 // This is the declaration of the MamboDecayer class.
13 //
14 
15 #include "HwDecayerBase.h"
16 #include "ThePEG/PDT/DecayMode.h"
17 
18 namespace Herwig {
19 using namespace ThePEG;
20 
27 class MamboDecayer: public HwDecayerBase {
28 
29 public:
30 
34  MamboDecayer() : _maxweight(10.), _a0(10,0.), _a1(10,0.) {}
35 
42  virtual bool accept(tcPDPtr parent, const tPDVector & children) const;
43 
50  virtual ParticleVector decay(const Particle & parent,
51  const tPDVector & children) const;
52 
58  virtual void dataBaseOutput(ofstream & os,bool header) const;
59 
60 
61 public:
62 
69  void persistentOutput(PersistentOStream & os) const;
70 
76  void persistentInput(PersistentIStream & is, int version);
78 
85  static void Init();
86 
87 protected:
88 
95  virtual IBPtr clone() const {return new_ptr(*this);}
96 
101  virtual IBPtr fullclone() const {return new_ptr(*this);}
103 
104 protected:
105 
112  virtual void doinitrun();
114 
115 private:
116 
121  MamboDecayer & operator=(const MamboDecayer &) = delete;
122 
123 private:
124 
131  double calculateMomentum(vector<Lorentz5Momentum> & mom,
132  Energy comEn) const;
133 
139  void colourConnections(const Particle & parent,
140  ParticleVector & out) const;
141 
151  void BesselFns(const long double x,
152  long double & f, long double & fp) const {
153  assert(x>=0.);
154  if( x < 10. ) {
155  f = BesselK0(x)/BesselK1(x);
156  fp = ( sqr(f)*x + f - x )/x;
157  }
158  else
159  BesselIExpand(-x, f, fp);
160  }
161 
169  void BesselIExpand(const long double x,
170  long double & f, long double & fp) const {
171  long double y = 1./x;
172  f = 1.+ y*(_a0[0] + y*(_a0[1] + y*(_a0[2] + y*(_a0[3]
173  + y*(_a0[4] + y*(_a0[5] + y*(_a0[6] + y*(_a0[7]
174  + y*(_a0[8] + y*_a0[9] )))))))));
175  fp = -y*y*(_a1[0] + y*(_a1[1] + y*(_a1[2] + y*(_a1[3]
176  + y*(_a1[4] + y*(_a1[5] + y*(_a1[6] + y*(_a1[7]
177  + y*(_a1[8] + y*_a1[9] )))))))));
178  }
179 
184  long double BesselI0(const long double x) const {
185  long double y,ans;
186  if(x < 3.75) {
187  y = sqr(x/3.75);
188  ans = 1. + y*(3.5156229 + y*(3.0899424 + y*(1.2067492
189  + y*(0.2659732 + y*(0.0360768+y*0.0045813)))));
190  }
191  else {
192  y = (3.75/x);
193  ans = (exp(x)/sqrt(x))*(0.39894228 + y*(0.01328592
194  + y*(0.00225319 + y*(-0.00157565 + y*(0.00916281
195  + y*(-0.02057706+y*(0.02635537+y*(-0.01647633+y*0.00392377))))))));
196  }
197  return ans;
198  }
199 
204  long double BesselI1(const long double x) const {
205  long double y,ans;
206  if(x < 3.75) {
207  y = sqr(x/3.75);
208  ans = x*(0.5 + y*(0.87890594 + y*(0.51498869 + y*(0.15084934
209  + y*(0.02658733 + y*(0.00301532 + y*0.00032411))))));
210  }
211  else {
212  y = 3.75/x;
213  ans = (0.39894228 + y*(-0.03988024 + y*(-0.00362018
214  + y*(0.00163801 + y*(-0.01031555 + y*(0.02282967
215  + y*(-0.02895312 + y*(0.01787654-y*0.00420059))))))))*(exp(x)/sqrt(x));
216  }
217  return ans;
218  }
219 
224  long double BesselK0(const long double x) const {
225  long double y,ans;
226  if(x <= 2.0) {
227  y = x*x/4.0;
228  ans = -log(x/2.0)*BesselI0(x) - 0.57721566
229  + y*(0.42278420 + y*(0.23069756
230  + y*(0.03488590 + y*(0.00262698 + y*(0.00010750+y*0.00000740)))));
231  }
232  else {
233  y = 2.0/x;
234  ans = (1.25331414 + y*(-0.07832358 + y*(+0.02189568
235  + y*(-0.01062446 + y*(0.00587872
236  + y*(-0.00251540 + y*0.00053208))))))*(exp(-x)/sqrt(x));
237  }
238  return ans;
239  }
240 
245  long double BesselK1(const long double x) const {
246  long double y,ans;
247  if(x <= 2.0) {
248  y = x*x/4.;
249  ans = log(x/2.)*BesselI1(x) + (1./x)*(1. + y*(0.15443144
250  + y*(-0.67278579 + y*(-0.18156897
251  + y*(-0.01919402+y*(-0.00110404-(y*0.00004686)))))));
252  }
253  else {
254  y = 2./x;
255  ans = (exp(-x)/sqrt(x))*(1.25331414 + y*(0.23498619
256  + y*(-0.03655620 + y*(0.01504268 + y*(-0.00780353
257  + y*(0.00325614+y*(-0.00068245)))))));
258  }
259  return ans;
260  }
262 
263 private:
264 
268  double _maxweight;
269 
273  vector<double> _a0;
274 
279  vector<double> _a1;
280 };
281 
282 }
283 
284 #endif /* HERWIG_MamboDecayer_H */
double _maxweight
Maximum weight.
Definition: MamboDecayer.h:268
double sqrt(int x)
long double BesselI0(const long double x) const
Modified Bessel function of first kind .
Definition: MamboDecayer.h:184
long double BesselI1(const long double x) const
Modified Bessel function of first kind .
Definition: MamboDecayer.h:204
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
constexpr auto sqr(const T &x) -> decltype(x *x)
The MamboDecayer class inherits from the Decayer class in ThePEG and implements the algorithm of R...
Definition: MamboDecayer.h:27
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
Definition: MamboDecayer.h:101
void BesselIExpand(const long double x, long double &f, long double &fp) const
Compute the values and it&#39;s derivative using asymptotic expansion.
Definition: MamboDecayer.h:169
long double BesselK1(const long double x) const
Modified Bessel function of second kind .
Definition: MamboDecayer.h:245
long double BesselK0(const long double x) const
Modified Bessel function of second kind .
Definition: MamboDecayer.h:224
vector< tPDPtr > tPDVector
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
MamboDecayer()
The default constructor.
Definition: MamboDecayer.h:34
vector< double > _a0
Store coefficents for aysymptotic expansion of .
Definition: MamboDecayer.h:273
virtual IBPtr clone() const
Make a simple clone of this object.
Definition: MamboDecayer.h:95
The HwDecayerBase class is the base class for Decayers in Herwig.
Definition: HwDecayerBase.h:37
-*- C++ -*-
vector< double > _a1
Store data for aysymptotic expansion of the first derivative .
Definition: MamboDecayer.h:279
void BesselFns(const long double x, long double &f, long double &fp) const
Compute the values and it&#39;s derivative using asymptotic expansion for large x values.
Definition: MamboDecayer.h:151
vector< PPtr > ParticleVector