herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
a1ThreePionDecayer.h
1// -*- C++ -*-
2//
3// a1ThreePionDecayer.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_a1ThreePionDecayer_H
10#define HERWIG_a1ThreePionDecayer_H
11//
12// This is the declaration of the a1ThreePionDecayer class.
13//
14#include "Herwig/Decay/DecayIntegrator.h"
15#include "Herwig/Decay/PhaseSpaceMode.h"
16#include "Herwig/Utilities/Kinematics.h"
17#include "ThePEG/Helicity/LorentzPolarizationVector.h"
18
19namespace Herwig {
20
21using namespace ThePEG;
22
103
104public:
105
110
117 virtual int modeNumber(bool & cc, tcPDPtr parent,
118 const tPDVector & children) const;
119
129 double me2(const int ichan,const Particle & part,
130 const tPDVector & outgoing,
131 const vector<Lorentz5Momentum> & momenta,
132 MEOption meopt) const;
133
137 virtual void constructSpinInfo(const Particle & part,
139
145 virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
146
160 virtual double threeBodyMatrixElement(const int imode , const Energy2 q2,
161 const Energy2 s3, const Energy2 s2,
162 const Energy2 s1, const Energy m1,
163 const Energy m2, const Energy m3) const;
164
170 virtual void dataBaseOutput(ofstream & os,bool header) const;
171
172public:
173
181
187 void persistentInput(PersistentIStream & is, int version);
189
193 static void Init();
194
195protected:
196
203 virtual IBPtr clone() const {return new_ptr(*this);}
204
209 virtual IBPtr fullclone() const {return new_ptr(*this);}
211
212protected:
213
221 virtual void doinit();
222
226 virtual void doinitrun();
228
229private:
230
235
236private:
237
243 Complex sigmaBreitWigner(Energy2 q2) const {
244 Energy q=sqrt(q2);
246 Energy2 msigma2=_sigmamass*_sigmamass;
247 Complex ii(0.,1.);
248 complex<Energy2> denom = q>2.*_mpi ? q2-msigma2+ii*msigma2*width/q :
249 q2-msigma2;
250 return msigma2/denom;
251 }
252
258 double a1FormFactor(Energy2 q2) const {
259 return (1.+_a1mass2/_lambda2)/(1.+q2/_lambda2);
260 }
261
268 Complex rhoBreitWigner(Energy2 q2,int ires) const {
269 Energy q=sqrt(q2);
270 Energy2 grhom = 8.*_prho[ires]*_prho[ires]*_prho[ires]/_rhomass[ires];
271 complex<Energy2> denom;
272 Complex ii(0.,1.);
273 if(q2<4.*_mpi2) {
274 denom=q2-_rhomass[ires]*_rhomass[ires]-_rhowidth[ires]*_rhomass[ires]*
275 (hFunction(q)-_hm2[ires]-(q2-_rhomass[ires]*_rhomass[ires])*_dhdq2m2[ires])
276 /grhom;
277 }
278 else {
279 Energy pcm=2.*Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
280 Energy2 grho = pcm*pcm*pcm/q;
281 denom=q2-_rhomass[ires]*_rhomass[ires]
282 -_rhowidth[ires]*_rhomass[ires]*
283 (hFunction(q)-_hm2[ires]-(q2-_rhomass[ires]*_rhomass[ires])*_dhdq2m2[ires])/grhom
284 +ii*_rhomass[ires]*_rhowidth[ires]*grho/grhom;
285 }
286 return _rhoD[ires]/denom;
287 }
288
294 Energy2 DParameter(int ires) const {
295 Energy2 grhom = 8.*_prho[ires]*_prho[ires]*_prho[ires]/_rhomass[ires];
296 return _rhomass[ires]*_rhomass[ires]+_rhowidth[ires]*_rhomass[ires]*
297 (hFunction(ZERO)-_hm2[ires]+sqr(_rhomass[ires])*_dhdq2m2[ires])/grhom;
298 }
299
305 double dhdq2Parameter(int ires) const {
306 Energy2 mrho2(sqr(_rhomass[ires]));
307 double root = sqrt(1.-4.*_mpi2/mrho2);
308 using Constants::pi;
309 return root/pi*(root+(1.+2*_mpi2/mrho2)*log((1+root)/(1-root)));
310 }
311
317 Energy2 hFunction(const Energy q) const {
318 static const Energy2 eps(0.01*MeV2);
319 Energy2 q2=sqr(q), output;
320 if(q2>4*_mpi2) {
321 double root = sqrt(1.-4.*_mpi2/q2);
322 output=root*log((1.+root)/(1.-root))*(q2-4*_mpi2)/Constants::pi;
323 }
324 else if(q2>eps) output=ZERO;
325 else output=-8.*_mpi2/Constants::pi;
326 return output;
327 }
328
332 Energy4 lambda(Energy2 a, Energy2 b, Energy2 c) const {
333 return sqr(a)+sqr(b)+sqr(c)-2.*a*b-2.*a*c-2.*b*c;
334 }
335
336private:
337
341 vector<Energy> _rhomass;
342
346 vector<Energy> _rhowidth;
347
351 vector<Energy> _prho;
352
356 vector<Energy2> _hm2;
357
361 vector<Energy2> _rhoD;
362
367 vector<double> _dhdq2m2;
368
373
378
382 Energy _psigma;
383
387 Energy _mpi;
388
392 Energy2 _mpi2;
393
397 Energy2 _lambda2;
398
402 Energy2 _a1mass2;
403
408
412 double _zmag;
413
417 double _zphase;
418
422 vector<Complex> _rhocoupling;
423
427 vector<double> _rhomag;
428
432 vector<double> _rhophase;
433
437 double _coupling;
438
443
447 mutable vector<double> _zerowgts;
448
452 mutable vector<double> _onewgts;
453
457 mutable vector<double> _twowgts;
458
462 mutable vector<double> _threewgts;
463
467 mutable double _zeromax;
468
472 mutable double _onemax;
473
477 mutable double _twomax;
478
482 mutable double _threemax;
483
488
492 mutable vector<Helicity::LorentzPolarizationVector> _vectors;
493
494};
495
496}
497
498
499#endif /* HERWIG_a1ThreePionDecayer_H */
Main class for Decayers implementing multi-channel phase space integration.
int imode() const
The mode being used for this decay.
MEOption
Enum for the matrix element option.
The a1ThreePionDecayer class is designed to implement the decay of the a_1 to three pions.
Energy2 _a1mass2
The mass squared of the meson, .
static void Init()
Standard Init function used to initialize the interfaces.
Energy _sigmawidth
The width of the meson.
bool _localparameters
use local values of the mass parameters
vector< double > _threewgts
Weights for the channels for the three charged pion channel.
Energy2 hFunction(const Energy q) const
The function in the propagator.
double dhdq2Parameter(int ires) const
The function in the rho propagator evaluated at .
Energy2 _mpi2
The mass of the pion, .
double me2(const int ichan, const Particle &part, const tPDVector &outgoing, const vector< Lorentz5Momentum > &momenta, MEOption meopt) const
Return the matrix element squared for a given mode and phase-space channel.
vector< Energy > _rhowidth
Width of the rho resonaces.
virtual IBPtr clone() const
Make a simple clone of this object.
vector< double > _dhdq2m2
The function in the rho propagator evaluated at for the different multiplets.
Energy2 DParameter(int ires) const
Normalisation factor for the propagator to ensure .
virtual void constructSpinInfo(const Particle &part, ParticleVector outgoing) const
Construct the SpinInfos for the particles produced in the decay.
vector< double > _rhophase
Phase of the rho coupling.
vector< Helicity::LorentzPolarizationVector > _vectors
Polarization vectors.
Energy2 _lambda2
The parameter for the form factor.
Energy _psigma
The momenta of the pions produced in the meson decay.
virtual void doinitrun()
Initialize this object to the begining of the run phase.
virtual double threeBodyMatrixElement(const int imode, const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy m1, const Energy m2, const Energy m3) const
The matrix element to be integrated for the three-body decays as a function of the invariant masses o...
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
virtual void doinit()
Initialize this object after the setup phase before saving and EventGenerator to disk.
double _twomax
Maximum weight for the two charged pion channel.
a1ThreePionDecayer()
Default constructor.
vector< double > _twowgts
Weights for the channels for the two charged pion channel.
double _zmag
The magnitude of the coupling.
vector< Energy > _rhomass
Mass of the rho resonances.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
virtual void dataBaseOutput(ofstream &os, bool header) const
Output the setup information for the particle database.
vector< double > _onewgts
Weights for the channels for the one charged pion channel.
virtual int modeNumber(bool &cc, tcPDPtr parent, const tPDVector &children) const
Which of the possible decays is required.
double _zphase
The phase of the coupling.
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode &dm) const
Method to return an object to calculate the 3 body partial width.
Energy _sigmamass
The mass of the meson.
double _threemax
Maximum weight for the three charged pion channel.
vector< Energy2 > _rhoD
The normalization factor for the propagator factor.
RhoDMatrix _rho
Spin density matrix.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
Complex rhoBreitWigner(Energy2 q2, int ires) const
Breit-Wigner for the , this is .
vector< double > _rhomag
Magnitude of the rho coupling.
double a1FormFactor(Energy2 q2) const
The form factor, .
Energy _mpi
The mass of the pion, .
double _onemax
Maximum weight for the one charged pion channel.
Complex _zsigma
The coupling for the resonance.
vector< Energy2 > _hm2
The function evaluated at .
Complex sigmaBreitWigner(Energy2 q2) const
Breit-wigner for the , this is .
double _zeromax
Maximum weight for the zero charged pion channel.
a1ThreePionDecayer & operator=(const a1ThreePionDecayer &)=delete
Private and non-existent assignment operator.
vector< Energy > _prho
Momentum of the pions produced in the decay.
vector< Complex > _rhocoupling
is the coupling of the th multiplet.
double _coupling
The overall coupling for the decay.
Energy4 lambda(Energy2 a, Energy2 b, Energy2 c) const
Momentum Function.
vector< double > _zerowgts
Weights for the channels for the zero charged pion channel.
Energy pstarTwoBodyDecay(const Energy M, const Energy m1, const Energy m2)
For the two body decay M -> m1 + m2 it gives the module of the 3-momentum of the decay product in the...
Definition: Kinematics.h:93
-*- C++ -*-
constexpr double pi
double sqrt(int x)
std::complex< double > Complex
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
vector< tPDPtr > tPDVector
vector< PPtr > ParticleVector
constexpr ZeroUnit ZERO
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
constexpr auto sqr(const T &x) -> decltype(x *x)