herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
18namespace Herwig {
19using namespace ThePEG;
20
28
29public:
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
61public:
62
70
76 void persistentInput(PersistentIStream & is, int version);
78
85 static void Init();
86
87protected:
88
95 virtual IBPtr clone() const {return new_ptr(*this);}
96
101 virtual IBPtr fullclone() const {return new_ptr(*this);}
103
104protected:
105
112 virtual void doinitrun();
114
115private:
116
122
123private:
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
263private:
264
269
273 vector<double> _a0;
274
279 vector<double> _a1;
280};
281
282}
283
284#endif /* HERWIG_MamboDecayer_H */
The HwDecayerBase class is the base class for Decayers in Herwig.
Definition: HwDecayerBase.h:37
The MamboDecayer class inherits from the Decayer class in ThePEG and implements the algorithm of R....
Definition: MamboDecayer.h:27
long double BesselK0(const long double x) const
Modified Bessel function of second kind .
Definition: MamboDecayer.h:224
static void Init()
The standard Init function used to initialize the interfaces.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
Definition: MamboDecayer.h:101
virtual void dataBaseOutput(ofstream &os, bool header) const
Output the setup information for the particle database.
double _maxweight
Maximum weight.
Definition: MamboDecayer.h:268
vector< double > _a0
Store coefficents for aysymptotic expansion of .
Definition: MamboDecayer.h:273
vector< double > _a1
Store data for aysymptotic expansion of the first derivative .
Definition: MamboDecayer.h:279
double calculateMomentum(vector< Lorentz5Momentum > &mom, Energy comEn) const
Set array of mometum to particles.
MamboDecayer & operator=(const MamboDecayer &)=delete
The assignment operator is private and must never be called.
virtual ParticleVector decay(const Particle &parent, const tPDVector &children) const
Perform the decay of the particle to the specified decay products.
MamboDecayer()
The default constructor.
Definition: MamboDecayer.h:34
long double BesselI0(const long double x) const
Modified Bessel function of first kind .
Definition: MamboDecayer.h:184
virtual void doinitrun()
Initialize this object.
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
virtual bool accept(tcPDPtr parent, const tPDVector &children) const
Check if this decayer can perfom the decay for a particular mode.
virtual IBPtr clone() const
Make a simple clone of this object.
Definition: MamboDecayer.h:95
void BesselFns(const long double x, long double &f, long double &fp) const
Compute the values and it's derivative using asymptotic expansion for large x values.
Definition: MamboDecayer.h:151
long double BesselK1(const long double x) const
Modified Bessel function of second kind .
Definition: MamboDecayer.h:245
void colourConnections(const Particle &parent, ParticleVector &out) const
Set up the colour connections for the decay.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
void BesselIExpand(const long double x, long double &f, long double &fp) const
Compute the values and it's derivative using asymptotic expansion.
Definition: MamboDecayer.h:169
long double BesselI1(const long double x) const
Modified Bessel function of first kind .
Definition: MamboDecayer.h:204
-*- C++ -*-
double sqrt(int x)
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
vector< tPDPtr > tPDVector
vector< PPtr > ParticleVector
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
constexpr auto sqr(const T &x) -> decltype(x *x)