herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
ThreeMesonDefaultCurrent.h
1 // -*- C++ -*-
2 //
3 // ThreeMesonDefaultCurrent.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_ThreeMesonDefaultCurrent_H
10 #define HERWIG_ThreeMesonDefaultCurrent_H
11 //
12 // This is the declaration of the ThreeMesonDefaultCurrent class.
13 //
14 #include "ThreeMesonCurrentBase.h"
15 #include "Herwig/Utilities/Interpolator.h"
16 #include "Herwig/Utilities/Kinematics.h"
17 #include "ThePEG/StandardModel/StandardModelBase.h"
18 
19 namespace Herwig {
20 using namespace ThePEG;
21 
49 
54  friend class Defaulta1MatrixElement;
55 
56 public:
57 
62 
69  void persistentOutput(PersistentOStream & os) const;
70 
76  void persistentInput(PersistentIStream & is, int version);
78 
82  static void Init();
83 
84 public:
85 
88 
103  virtual bool createMode(int icharge,unsigned int imode,DecayPhaseSpaceModePtr mode,
104  unsigned int iloc,unsigned int ires,
105  DecayPhaseSpaceChannelPtr phase,Energy upp);
107 
114  virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
115 
128  double threeBodyMatrixElement(const int imode, const Energy2 q2,
129  const Energy2 s3, const Energy2 s2,
130  const Energy2 s1, const Energy m1,
131  const Energy m2, const Energy m3) const;
132 
133 protected:
134 
139  virtual bool acceptMode(int) const;
140 
151  virtual FormFactors calculateFormFactors(const int ichan, const int imode,
152  Energy2 q2,
153  Energy2 s1, Energy2 s2, Energy2 s3) const;
154 
155 protected:
156 
163  virtual IBPtr clone() const {return new_ptr(*this);}
164 
169  virtual IBPtr fullclone() const {return new_ptr(*this);}
171 
172 protected:
173 
181  virtual void doinit();
182 
186  virtual void doinitrun();
187 
191  virtual void doupdate();
193 
194 private:
195 
199  ThreeMesonDefaultCurrent & operator=(const ThreeMesonDefaultCurrent &) = delete;
200 
201 private:
202 
209  Complex BrhoF123(Energy2 q2,int ires) const {
210  Complex output(0.),norm(0.);
211  for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
212  norm+=_rhoF123wgts[ix];
213  }
214  if(ires<0) {
215  for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
216  output+=_rhoF123wgts[ix]*rhoKBreitWigner(q2,0,ix);
217  }
218  }
219  else {
220  unsigned int temp(ires);
221  if(temp<_rhoF123wgts.size()&&temp<3)
222  output=_rhoF123wgts[temp]*rhoKBreitWigner(q2,0,temp);
223  else
224  output=0.;
225  }
226  return output/norm;
227  }
228 
235  Complex BrhoF5(Energy2 q2,int ires) const {
236  Complex output(0.),norm(0.);
237  for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
238  norm+=_rhoF5wgts[ix];
239  }
240  if(ires<0) {
241  for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
242  output+=_rhoF5wgts[ix]*rhoKBreitWigner(q2,1,ix);
243  }
244  }
245  else {
246  unsigned int temp(ires);
247  if(temp<_rhoF5wgts.size()&&temp<3) {
248  output=_rhoF5wgts[temp]*rhoKBreitWigner(q2,1,temp);
249  }
250  }
251  return output/norm;
252  }
253 
260  Complex BKstarF123(Energy2 q2,int ires) const {
261  Complex output(0.),norm(0.);
262  for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
263  norm+=_kstarF123wgts[ix];
264  }
265  if(ires<0) {
266  for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
267  output+=_kstarF123wgts[ix]*rhoKBreitWigner(q2,2,ix);
268  }
269  }
270  else {
271  unsigned int temp(ires);
272  if(temp<_kstarF123wgts.size()&&temp<3) {
273  output=_kstarF123wgts[temp]*rhoKBreitWigner(q2,2,temp);
274  }
275  }
276  return output/norm;
277  }
278 
285  Complex BKstarF5(Energy2 q2,int ires) const {
286  Complex output(0.),norm(0.);
287  for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
288  norm+=_kstarF5wgts[ix];
289  }
290  if(ires<0) {
291  for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
292  output+=_kstarF5wgts[ix]*rhoKBreitWigner(q2,3,ix);
293  }
294  }
295  else {
296  unsigned int temp(ires);
297  if(temp<_kstarF5wgts.size()&&temp<3) {
298  output=_kstarF5wgts[ires]*rhoKBreitWigner(q2,3,temp);
299  }
300  }
301  return output/norm;
302  }
303 
311  Complex FKrho(Energy2 si,Energy2 sj,int ires) const {
312  Complex output;
313  if(ires<0){output = _rhoKstarwgt*BKstarF123(si,-1)+BrhoF123(sj,-1);}
314  else if(ires%2==0){output= _rhoKstarwgt*BKstarF123(si,ires/2);}
315  else if(ires%2==1){output=BrhoF123(sj,ires/2);}
316  output /=(1.+_rhoKstarwgt);
317  return output;
318  }
319 
325  Complex a1BreitWigner(Energy2 q2) const {
326  Complex ii(0.,1.);
327  Energy2 m2(_a1mass*_a1mass);
328  Energy q(sqrt(q2));
329  return m2/(m2-q2-ii*q*a1Width(q2));
330  }
331 
337  Complex K1BreitWigner(Energy2 q2) const {
338  Energy2 m2 = sqr(_k1mass);
339  Complex ii(0.,1.);
340  complex<Energy2> fact(m2 - ii*_k1mass*_k1width);
341  return fact/(fact-q2);
342  }
343 
349  Energy a1Width(Energy2 q2) const {
350  Energy output;
351  if(!_a1opt) output = _a1mass*_a1width*g(q2)/g(sqr(_a1mass))/sqrt(q2);
352  else output = (*_a1runinter)(q2);
353  return output;
354  }
355 
359  double g(Energy2 q2) const {
360  double output;
361  if(q2 < 9.*sqr(_mpi)) {
362  output=0.;
363  }
364  else if(q2 < sqr(_rhoF123masses[0]+_mpi)) {
365  double diff = (q2-9.*sqr(_mpi))/GeV2;
366 
367  output = 4.1*sqr(diff)*diff*(1.-3.3*diff+5.8*sqr(diff));
368  }
369  else {
370  double ratio = q2/GeV2;
371  output = ratio*(1.623+10.38/ratio-9.32/sqr(ratio)+0.65/(ratio*sqr(ratio)));
372  }
373  return output;
374  }
375 
380  void inita1Width(int iopt);
381 
388  Complex rhoKBreitWigner(Energy2 q2,unsigned int itype,unsigned int ires) const;
389 
390 private:
391 
396  vector<double> _rhoF123wgts;
397 
402  vector<double> _kstarF123wgts;
403 
408  vector<double> _rhoF5wgts;
409 
414  vector<double> _kstarF5wgts;
415 
419  double _rhoKstarwgt;
420 
424  vector<Energy> _a1runwidth;
425 
429  vector<Energy2> _a1runq2;
430 
431 
436 
441 
445  Energy _a1mass;
446 
450  Energy _a1width;
451 
455  Energy _k1mass;
456 
460  Energy _k1width;
461 
465  Energy _fpi;
466 
470  Energy _mpi;
471 
475  Energy _mK;
476 
481 
485  vector<Energy> _rhoF123masses;
486 
490  vector<Energy> _rhoF5masses;
491 
495  vector<Energy> _rhoF123widths;
496 
500  vector<Energy> _rhoF5widths;
501 
506 
510  vector<Energy> _kstarF123masses;
511 
515  vector<Energy> _kstarF5masses;
516 
520  vector<Energy> _kstarF123widths;
521 
525  vector<Energy> _kstarF5widths;
526 
531 
536 
540  bool _a1opt;
541 
545  Energy _maxmass;
546 
550  Energy _maxcalc;
551 
552 };
553 
554 }
555 
556 #endif /* THEPEG_ThreeMesonDefaultCurrent_H */
Energy _k1mass
The mass of the resonances.
Energy a1Width(Energy2 q2) const
The running width.
double sqrt(int x)
std::complex< double > Complex
Energy _maxcalc
The maximum mass when the running width was calculated.
vector< Energy > _kstarF123masses
The masses for the form factors.
vector< double > _rhoF5wgts
Parameters for the Breit-Wigner in the form factors.
Interpolator< Energy, Energy2 >::Ptr _a1runinter
The interpolator for the running width calculation.
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
vector< Energy > _kstarF123widths
The widths for the form factors.
vector< Energy > _rhoF123masses
The masses for the form factors.
constexpr auto sqr(const T &x) -> decltype(x *x)
bool _a1parameters
Use local values of the parameters.
Complex BKstarF123(Energy2 q2, int ires) const
The Breit-Wigner for the form factors.
virtual IBPtr fullclone() const
Make a clone of this object, possibly modifying the cloned object to make it sane.
double _rhoKstarwgt
The relative weight of the and where needed.
Energy _a1mass
The mass of the resonances.
vector< Energy > _rhoF5masses
The masses for the form factors.
vector< double > _rhoF123wgts
Parameters for the Breit-Wigner in the form factors.
vector< Energy > _a1runwidth
The width for the running width calculation.
double g(Energy2 q2) const
The function of Kuhn and Santamaria.
vector< Energy > _kstarF5masses
The masses for the form factors.
Complex K1BreitWigner(Energy2 q2) const
The Breit-Wigner.
Energy _fpi
The pion decay constant, .
virtual IBPtr clone() const
Make a simple clone of this object.
Complex BrhoF5(Energy2 q2, int ires) const
The Breit-Wigner for the form factors.
bool _rhoparameters
use local values of the masses and widths
Complex BKstarF5(Energy2 q2, int ires) const
The Breit-Wigner for the form factors.
bool _kstarparameters
use local values of the resonances masses and widths
vector< Energy > _kstarF5widths
The widths for the form factors.
Energy _a1width
The width of the resonances.
-*- C++ -*-
vector< double > _kstarF123wgts
Parameters for the Breit-Wigner in the form factors.
vector< double > _kstarF5wgts
Parameters for the Breit-Wigner in the form factors.
Complex a1BreitWigner(Energy2 q2) const
Breit-Wigner
vector< Energy > _rhoF123widths
The widths for the form factors.
The ThreeMesonDefaultCurrent class implements the currents from Z.Phys.C58:445 (1992), this paper uses the form from Z.Phys.C48:445 (1990) for the width and is the default model in TAUOLA.
vector< Energy > _rhoF5widths
The widths for the form factors.
Complex FKrho(Energy2 si, Energy2 sj, int ires) const
Mixed Breit Wigner for the form factor.
Energy _maxmass
The maximum mass of the hadronic system.
vector< Energy2 > _a1runq2
The for the running width calculation.
Energy _k1width
The width of the resonances.
This is the base class for the three meson decays of the weak current.
bool _initializea1
Initialize the running width.
Complex BrhoF123(Energy2 q2, int ires) const
The Breit-Wigner for the form factors.
bool _k1parameters
Use local values of the parameters.