herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
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/DecayPhaseSpaceMode.h"
16 #include "Herwig/Utilities/Kinematics.h"
17 #include "ThePEG/Helicity/LorentzPolarizationVector.h"
18 
19 namespace Herwig {
20 
21 using namespace ThePEG;
22 
103 
104 public:
105 
110 
117  virtual int modeNumber(bool & cc, tcPDPtr parent,
118  const tPDVector & children) const;
119 
128  double me2(const int ichan,const Particle & part,
129  const ParticleVector & decay, MEOption meopt) const;
130 
136  virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
137 
151  virtual double threeBodyMatrixElement(const int imode , const Energy2 q2,
152  const Energy2 s3, const Energy2 s2,
153  const Energy2 s1, const Energy m1,
154  const Energy m2, const Energy m3) const;
155 
161  virtual void dataBaseOutput(ofstream & os,bool header) const;
162 
163 public:
164 
171  void persistentOutput(PersistentOStream & os) const;
172 
178  void persistentInput(PersistentIStream & is, int version);
180 
184  static void Init();
185 
186 protected:
187 
194  virtual IBPtr clone() const {return new_ptr(*this);}
195 
200  virtual IBPtr fullclone() const {return new_ptr(*this);}
202 
203 protected:
204 
212  virtual void doinit();
213 
217  virtual void doinitrun();
219 
220 private:
221 
225  a1ThreePionDecayer & operator=(const a1ThreePionDecayer &) = delete;
226 
227 private:
228 
234  Complex sigmaBreitWigner(Energy2 q2) const {
235  Energy q=sqrt(q2);
236  Energy width=_sigmawidth*Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi)/_psigma;
237  Energy2 msigma2=_sigmamass*_sigmamass;
238  Complex ii(0.,1.);
239  complex<Energy2> denom = q>2.*_mpi ? q2-msigma2+ii*msigma2*width/q :
240  q2-msigma2;
241  return msigma2/denom;
242  }
243 
249  double a1FormFactor(Energy2 q2) const {
250  return (1.+_a1mass2/_lambda2)/(1.+q2/_lambda2);
251  }
252 
259  Complex rhoBreitWigner(Energy2 q2,int ires) const {
260  Energy q=sqrt(q2);
261  Energy2 grhom = 8.*_prho[ires]*_prho[ires]*_prho[ires]/_rhomass[ires];
262  complex<Energy2> denom;
263  Complex ii(0.,1.);
264  if(q2<4.*_mpi2) {
265  denom=q2-_rhomass[ires]*_rhomass[ires]-_rhowidth[ires]*_rhomass[ires]*
266  (hFunction(q)-_hm2[ires]-(q2-_rhomass[ires]*_rhomass[ires])*_dhdq2m2[ires])
267  /grhom;
268  }
269  else {
270  Energy pcm=2.*Kinematics::pstarTwoBodyDecay(q,_mpi,_mpi);
271  Energy2 grho = pcm*pcm*pcm/q;
272  denom=q2-_rhomass[ires]*_rhomass[ires]
273  -_rhowidth[ires]*_rhomass[ires]*
274  (hFunction(q)-_hm2[ires]-(q2-_rhomass[ires]*_rhomass[ires])*_dhdq2m2[ires])/grhom
275  +ii*_rhomass[ires]*_rhowidth[ires]*grho/grhom;
276  }
277  return _rhoD[ires]/denom;
278  }
279 
285  Energy2 DParameter(int ires) const {
286  Energy2 grhom = 8.*_prho[ires]*_prho[ires]*_prho[ires]/_rhomass[ires];
287  return _rhomass[ires]*_rhomass[ires]+_rhowidth[ires]*_rhomass[ires]*
288  (hFunction(ZERO)-_hm2[ires]+sqr(_rhomass[ires])*_dhdq2m2[ires])/grhom;
289  }
290 
296  double dhdq2Parameter(int ires) const {
297  Energy2 mrho2(sqr(_rhomass[ires]));
298  double root = sqrt(1.-4.*_mpi2/mrho2);
299  using Constants::pi;
300  return root/pi*(root+(1.+2*_mpi2/mrho2)*log((1+root)/(1-root)));
301  }
302 
308  Energy2 hFunction(const Energy q) const {
309  static const Energy2 eps(0.01*MeV2);
310  Energy2 q2=sqr(q), output;
311  double root = sqrt(1.-4.*_mpi2/q2);
312  if(q2>4*_mpi2) {
313  output=root*log((1.+root)/(1.-root))*(q2-4*_mpi2)/Constants::pi;
314  }
315  else if(q2>eps) output=ZERO;
316  else output=-8.*_mpi2/Constants::pi;
317  return output;
318  }
319 
323  Energy4 lambda(Energy2 a, Energy2 b, Energy2 c) const {
324  return sqr(a)+sqr(b)+sqr(c)-2.*a*b-2.*a*c-2.*b*c;
325  }
326 
327 private:
328 
332  vector<Energy> _rhomass;
333 
337  vector<Energy> _rhowidth;
338 
342  vector<Energy> _prho;
343 
347  vector<Energy2> _hm2;
348 
352  vector<Energy2> _rhoD;
353 
358  vector<double> _dhdq2m2;
359 
363  Energy _sigmamass;
364 
368  Energy _sigmawidth;
369 
373  Energy _psigma;
374 
378  Energy _mpi;
379 
383  Energy2 _mpi2;
384 
388  Energy2 _lambda2;
389 
393  Energy2 _a1mass2;
394 
399 
403  double _zmag;
404 
408  double _zphase;
409 
413  vector<Complex> _rhocoupling;
414 
418  vector<double> _rhomag;
419 
423  vector<double> _rhophase;
424 
428  double _coupling;
429 
434 
438  mutable vector<double> _zerowgts;
439 
443  mutable vector<double> _onewgts;
444 
448  mutable vector<double> _twowgts;
449 
453  mutable vector<double> _threewgts;
454 
458  mutable double _zeromax;
459 
463  mutable double _onemax;
464 
468  mutable double _twomax;
469 
473  mutable double _threemax;
474 
478  mutable RhoDMatrix _rho;
479 
483  mutable vector<Helicity::LorentzPolarizationVector> _vectors;
484 
485 };
486 
487 }
488 
489 
490 #endif /* HERWIG_a1ThreePionDecayer_H */
double _coupling
The overall coupling for the decay.
Energy2 _lambda2
The parameter for the form factor.
double sqrt(int x)
std::complex< double > Complex
Energy _sigmawidth
The width of the meson.
Energy2 _a1mass2
The mass squared of the meson, .
Energy _sigmamass
The mass of the meson.
double _zeromax
Maximum weight for the zero charged pion channel.
Energy4 lambda(Energy2 a, Energy2 b, Energy2 c) const
Momentum Function.
double dhdq2Parameter(int ires) const
The function in the rho propagator evaluated at .
double _threemax
Maximum weight for the three charged pion channel.
Energy2 DParameter(int ires) const
Normalisation factor for the propagator to ensure .
virtual IBPtr clone() const
Make a simple clone of this object.
Complex sigmaBreitWigner(Energy2 q2) const
Breit-wigner for the , this is .
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
vector< double > _dhdq2m2
The function in the rho propagator evaluated at for the different multiplets.
vector< Energy > _rhomass
Mass of the rho resonances.
constexpr double pi
vector< Energy > _prho
Momentum of the pions produced in the decay.
constexpr auto sqr(const T &x) -> decltype(x *x)
Energy2 hFunction(const Energy q) const
The function in the propagator.
Complex _zsigma
The coupling for the resonance.
vector< double > _rhophase
Phase of the rho coupling.
double _zphase
The phase of the coupling.
vector< double > _rhomag
Magnitude of the rho coupling.
MEOption
Enum for the matrix element option.
The a1ThreePionDecayer class is designed to implement the decay of the a_1 to three pions...
Energy2 _mpi2
The mass of the pion, .
vector< tPDPtr > tPDVector
double _zmag
The magnitude of the coupling.
vector< double > _onewgts
Weights for the channels for the one charged pion channel.
Complex rhoBreitWigner(Energy2 q2, int ires) const
Breit-Wigner for the , this is .
Energy _mpi
The mass of the pion, .
vector< double > _threewgts
Weights for the channels for the three charged pion channel.
vector< double > _twowgts
Weights for the channels for the two charged pion channel.
Main class for Decayers implementing multi-channel phase space integration.
ThePEG::Ptr< ParticleData >::transient_const_pointer tcPDPtr
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
vector< Complex > _rhocoupling
is the coupling of the th multiplet.
vector< Energy2 > _hm2
The function evaluated at .
vector< Energy > _rhowidth
Width of the rho resonaces.
vector< Helicity::LorentzPolarizationVector > _vectors
Polarization vectors.
vector< double > _zerowgts
Weights for the channels for the zero charged pion channel.
double a1FormFactor(Energy2 q2) const
The form factor, .
RhoDMatrix _rho
Spin density matrix.
-*- C++ -*-
double _twomax
Maximum weight for the two charged pion channel.
vector< PPtr > ParticleVector
Energy _psigma
The momenta of the pions produced in the meson decay.
bool _localparameters
use local values of the mass parameters
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
constexpr ZeroUnit ZERO
vector< Energy2 > _rhoD
The normalization factor for the propagator factor.
double _onemax
Maximum weight for the one charged pion channel.