herwig
is hosted by
Hepforge
,
IPPP Durham
Herwig
7.3.0
Decay
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 */
Herwig::HwDecayerBase
The HwDecayerBase class is the base class for Decayers in Herwig.
Definition:
HwDecayerBase.h:37
Herwig::MamboDecayer
The MamboDecayer class inherits from the Decayer class in ThePEG and implements the algorithm of R....
Definition:
MamboDecayer.h:27
Herwig::MamboDecayer::BesselK0
long double BesselK0(const long double x) const
Modified Bessel function of second kind .
Definition:
MamboDecayer.h:224
Herwig::MamboDecayer::Init
static void Init()
The standard Init function used to initialize the interfaces.
Herwig::MamboDecayer::fullclone
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
Definition:
MamboDecayer.h:101
Herwig::MamboDecayer::dataBaseOutput
virtual void dataBaseOutput(ofstream &os, bool header) const
Output the setup information for the particle database.
Herwig::MamboDecayer::_maxweight
double _maxweight
Maximum weight.
Definition:
MamboDecayer.h:268
Herwig::MamboDecayer::_a0
vector< double > _a0
Store coefficents for aysymptotic expansion of .
Definition:
MamboDecayer.h:273
Herwig::MamboDecayer::_a1
vector< double > _a1
Store data for aysymptotic expansion of the first derivative .
Definition:
MamboDecayer.h:279
Herwig::MamboDecayer::calculateMomentum
double calculateMomentum(vector< Lorentz5Momentum > &mom, Energy comEn) const
Set array of mometum to particles.
Herwig::MamboDecayer::operator=
MamboDecayer & operator=(const MamboDecayer &)=delete
The assignment operator is private and must never be called.
Herwig::MamboDecayer::decay
virtual ParticleVector decay(const Particle &parent, const tPDVector &children) const
Perform the decay of the particle to the specified decay products.
Herwig::MamboDecayer::MamboDecayer
MamboDecayer()
The default constructor.
Definition:
MamboDecayer.h:34
Herwig::MamboDecayer::BesselI0
long double BesselI0(const long double x) const
Modified Bessel function of first kind .
Definition:
MamboDecayer.h:184
Herwig::MamboDecayer::doinitrun
virtual void doinitrun()
Initialize this object.
Herwig::MamboDecayer::persistentInput
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
Herwig::MamboDecayer::accept
virtual bool accept(tcPDPtr parent, const tPDVector &children) const
Check if this decayer can perfom the decay for a particular mode.
Herwig::MamboDecayer::clone
virtual IBPtr clone() const
Make a simple clone of this object.
Definition:
MamboDecayer.h:95
Herwig::MamboDecayer::BesselFns
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
Herwig::MamboDecayer::BesselK1
long double BesselK1(const long double x) const
Modified Bessel function of second kind .
Definition:
MamboDecayer.h:245
Herwig::MamboDecayer::colourConnections
void colourConnections(const Particle &parent, ParticleVector &out) const
Set up the colour connections for the decay.
Herwig::MamboDecayer::persistentOutput
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
Herwig::MamboDecayer::BesselIExpand
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
Herwig::MamboDecayer::BesselI1
long double BesselI1(const long double x) const
Modified Bessel function of first kind .
Definition:
MamboDecayer.h:204
ThePEG::Particle
ThePEG::PersistentIStream
ThePEG::PersistentOStream
Herwig
-*- C++ -*-
Definition:
BasicConsistency.h:17
ThePEG
ThePEG::sqrt
double sqrt(int x)
ThePEG::IBPtr
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
ThePEG::tPDVector
vector< tPDPtr > tPDVector
ThePEG::ParticleVector
vector< PPtr > ParticleVector
ThePEG::tcPDPtr
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
ThePEG::sqr
constexpr auto sqr(const T &x) -> decltype(x *x)
Generated on Thu Jun 20 2024 17:50:52 for Herwig by
1.9.6