herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
TwoMesonRhoKStarCurrent.h
1 // -*- C++ -*-
2 //
3 // TwoMesonRhoKStarCurrent.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_TwoMesonRhoKStarCurrent_H
10 #define HERWIG_TwoMesonRhoKStarCurrent_H
11 // This is the declaration of the TwoMesonRhoKStarCurrent class.
12 
13 #include "WeakDecayCurrent.h"
14 #include "ThePEG/PDT/EnumParticles.h"
15 #include "Herwig/Utilities/Kinematics.h"
16 #include "ThePEG/StandardModel/StandardModelBase.h"
17 
18 namespace Herwig {
19 using namespace ThePEG;
20 
54 
55 public:
56 
61 
79  virtual bool createMode(int icharge,unsigned int imode,DecayPhaseSpaceModePtr mode,
80  unsigned int iloc,unsigned int ires,
81  DecayPhaseSpaceChannelPtr phase,Energy upp);
82 
92  virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
94 
104  virtual vector<LorentzPolarizationVectorE>
105  current(const int imode, const int ichan,Energy & scale,
106  const ParticleVector & decay, DecayIntegrator::MEOption meopt) const;
107 
113  virtual bool accept(vector<int> id);
114 
120  virtual unsigned int decayMode(vector<int> id);
121 
128  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
129 
130 public:
131 
138  void persistentOutput(PersistentOStream & os) const;
139 
145  void persistentInput(PersistentIStream & is, int version);
147 
151  static void Init();
152 
153 protected:
154 
161  virtual IBPtr clone() const {return new_ptr(*this);}
162 
167  virtual IBPtr fullclone() const {return new_ptr(*this);}
169 
170 protected:
171 
174 
180  virtual void doinit();
182 
183 private:
184 
188  TwoMesonRhoKStarCurrent & operator=(const TwoMesonRhoKStarCurrent &) = delete;
189 
190 private:
191 
202  Complex BreitWigner(Energy2 q2, unsigned int imodel,unsigned int itype,
203  unsigned int ires) const {
204  // workout the index of the resonace
205  unsigned int iy(3*itype+ires);
206  // off shell mass
207  Energy q(sqrt(q2));
208  // and the running width
209  Energy pq(pcm(iy,q));
210  double ratio(pow<3,1>(pq/_mom[iy]));
211  Energy gamrun(_width[iy]*_mass[iy]*ratio/q);
212  // work out the denominator
213  complex<Energy2> denom, numer;
214  Complex ii(0.,1.);
215  if(imodel==0) {
216  denom=q2-_mass2[iy]+ii*_mass[iy]*gamrun;
217  numer=-_mass2[iy];
218  }
219  else if(imodel==1) {
220  denom = q2 - _mass2[iy] + ii*_mass[iy]*gamrun
221  - ( sqr(pq) * (GSModelhFunction(iy,q)-_hm2[iy])
222  - sqr(_mom[iy]) * (q2-_mass2[iy]) * _dhdq2[iy]);
223 
224  numer = -(_mass2[iy]+_dparam[iy]*_mass[iy]*_width[iy]);
225  }
226  else assert(false);
227  // return the answer
228  return numer/denom;
229  }
230 
236  double GSModelDParameter(const unsigned int ires)const {
237  Energy mpi(0.5*(_massa[ires]+_massb[ires]));
238  using Constants::pi;
239  return 3.*mpi*mpi/pi/sqr(_mom[ires])*log(0.5*(_mass[ires]+2.*_mom[ires])/mpi)
240  +0.5*_mass[ires]/pi/_mom[ires]-mpi*mpi*_mass[ires]/pi/pow<3,1>(_mom[ires]);
241  }
242 
249  InvEnergy2 GSModeldhdq2Parameter(const unsigned int ires)const {
250  Energy mpi = 0.5 * (_massa[ires] + _massb[ires]);
251  return _width[ires] / Constants::pi / pow<3,1>(_mom[ires]) *
252  (sqr(mpi) / _mass[ires] / _mom[ires] *
253  log(0.5 * (_mass[ires] + 2.*_mom[ires])/mpi) + 0.5);
254  }
255 
262  double GSModelhFunction(const unsigned int ires, const Energy q) const {
263  Energy pq(pcm(ires,q));
264  return _width[ires]*_mass2[ires]/pow<3,1>(_mom[ires])
265  *2.*pq/Constants::pi/q*log((q+2.*pq)/(_massa[ires]+_massb[ires]));
266  }
267 
275  Energy pcm(const unsigned int ires, const Energy q) const {
276  Energy2 q2(q*q);
277  if(q <= _massa[ires]+_massb[ires]) return ZERO;
278  return 0.5/q*sqrt((q2-sqr(_massa[ires]+_massb[ires]))*
279  (q2-sqr(_massa[ires]-_massb[ires])));
280  }
282 
283 private:
284 
292  vector<Complex> _piwgt;
293 
297  vector<double> _pimag;
298 
302  vector<double> _piphase;
304 
312  vector<Complex> _kwgt;
313 
317  vector<double> _kmag;
318 
322  vector<double> _kphase;
324 
328  int _pimodel;
329 
333  int _kmodel;
334 
339 
343  vector<Energy> _rhomasses;
344 
348  vector<Energy> _rhowidths;
349 
354 
358  vector<Energy> _kstarmasses;
359 
363  vector<Energy> _kstarwidths;
364 
372  vector<Energy> _mass;
373 
377  vector<Energy> _width;
378 
382  vector<Energy2> _mass2;
383 
387  vector<Energy2> _massw;
388 
392  vector<Energy> _massa,_massb;
393 
397  vector<Energy> _mom;
398 
403  vector<InvEnergy2> _dhdq2;
404 
409  vector<double> _hm2;
410 
415  vector<double> _dparam;
417 };
418 
419 }
420 
421 
422 #endif /* HERWIG_TwoMesonRhoKStarCurrent_H */
bool _rhoparameters
Option not to use the physical masses and widths for the .
The WeakDecayCurrent class is the base class for the hadronic currents produced in weak decays...
vector< double > _dparam
The parameter for the GS form of the Breit-Wigner.
double sqrt(int x)
vector< Energy2 > _mass2
squares of the masses
std::complex< double > Complex
Complex BreitWigner(Energy2 q2, unsigned int imodel, unsigned int itype, unsigned int ires) const
-wave breit wigner for form-factors
int _pimodel
Model to use for the propagator.
virtual IBPtr clone() const
Make a simple clone of this object.
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
vector< Energy > _rhomasses
The masses of the resonances.
constexpr double pi
vector< Complex > _piwgt
Weights for the different resonances in the current, .
vector< Energy > _rhowidths
The widths of the resonances.
constexpr auto sqr(const T &x) -> decltype(x *x)
vector< Energy > _mass
Parameters for the Breit-Wigners.
vector< double > _kmag
The magnitude for input.
double GSModelDParameter(const unsigned int ires) const
The parameter in the GS model of the propagator.
InvEnergy2 GSModeldhdq2Parameter(const unsigned int ires) const
The function in the GS model for the propagator evaluated at .
MEOption
Enum for the matrix element option.
int _kmodel
Model to use for the propagator.
vector< double > _hm2
The function at for the GS form of the Breit-Wigner.
vector< Energy > _width
The widths of the resonances.
vector< tPDPtr > tPDVector
vector< Energy > _kstarwidths
The masses of the resonances.
vector< double > _kphase
The phase for input.
vector< InvEnergy2 > _dhdq2
The function at for the GS form of the Breit-Wigner.
vector< Energy > _massa
Masses of the decay products for the momentum calculation.
vector< double > _pimag
The magnitude for input.
Energy pcm(const unsigned int ires, const Energy q) const
The momentum of the decay products for the running width.
vector< Energy2 > _massw
product of the mass and the width
-*- C++ -*-
vector< Complex > _kwgt
Weights for the different resonances in the current, .
vector< Energy > _mom
The decay for the on-shell mass momentum.
vector< PPtr > ParticleVector
bool _kstarparameters
Option not to use the physical masses and widths for the .
vector< Energy > _kstarmasses
The masses of the resonances.
double GSModelhFunction(const unsigned int ires, const Energy q) const
The function in the GS model.
Weak current for the production of two mesons via the or resonances.
constexpr ZeroUnit ZERO
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
vector< double > _piphase
The phase for input.