Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-06 09:24:05

0001 // -*- C++ -*-
0002 //
0003 // TwoKaonOnePionCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_TwoKaonOnePionCurrent_H
0010 #define HERWIG_TwoKaonOnePionCurrent_H
0011 //
0012 // This is the declaration of the TwoKaonOnePionCurrent class.
0013 //
0014 
0015 #include "WeakCurrent.h"
0016 
0017 namespace Herwig {
0018 
0019 using namespace ThePEG;
0020 
0021 /**
0022  * The TwoKaonOnePionCurrent class implements the model of M. Finkemeier 
0023  * and E.~Mirkes, Z. Phys. C 69 (1996) 243 [arXiv:hep-ph/9503474],
0024  * for the weak current for three mesons where at least one of the mesons is
0025  * a kaon.
0026  *
0027  * \ingroup Decay
0028  *
0029  *  This is the base class for the three meson decays of the weak current.
0030  *  It is designed so that the currents for the following modes can be implemented
0031  *  in classes inheriting from this
0032  * - \f$    K^-   \pi^-    K^+ \f$, (imode=0)
0033  * - \f$    K^0   \pi^-    \bar{K}^0\f$, (imode=1)
0034  * - \f$    K^-   \pi^0    K^0 \f$, (imode=2)
0035  * - \f$    \pi^0  \pi^0    K^- \f$, (imode=3)
0036  * - \f$    K^-   \pi^-    \pi^+ \f$, (imode=4)
0037  * - \f$    \pi^-  \bar{K}^0  \pi^0 \f$, (imode=5)
0038  * - \f$    \pi^-  \pi^0    \eta \f$, (imode=6)
0039  *
0040  * obviously there are other modes with three pseudoscalar mesons for the decay
0041  * of the weak current but this model original came from \f$\tau\f$ decay where
0042  * these are the only modes. However one case which is important is the inclusion
0043  * of the mixing in the neutral kaon sector for which we include the additional
0044  * currents
0045  * - \f$    K^0_S \pi^- K^0_S\f$, (imode=9)
0046  * - \f$    K^0_L \pi^- K^0_L\f$, (imode=10)
0047  * - \f$    K^0_S \pi^- K^0_L\f$, (imode=11)
0048  *
0049  *  In this case the current is given by
0050  *  \f[ J^\mu = \left(g^{\mu\nu}-\frac{q^\mu q^\nu}{q^2}\right)
0051  *   \left[F_1(p_2-p_3)^\mu +F_2(p_3-p_1)^\mu+F_3(p_1-p_2)^\mu\right]
0052  *  +q^\mu F_4
0053  *  +F_5\epsilon^{\mu\alpha\beta\gamma}p_1^\alpha p_2^\beta p_3^\gamma
0054  *  \f]
0055  * where
0056  * - \f$p_{1,2,3}\f$ are the momenta of the mesons in the order given above.
0057  * - \f$F_1,F_2,F_3,F_4,F_5\f$ are the form factors which must be 
0058  *  calculated in the calculateFormFactors member which should be implemented
0059  * in classes inheriting from this.
0060  *
0061  * @see WeakCurrent.
0062  *  
0063  * \author Peter Richardson
0064  * @see \ref TwoKaonOnePionCurrentInterfaces "The interfaces"
0065  * defined for TwoKaonOnePionCurrent.
0066  */
0067 class TwoKaonOnePionCurrent: public WeakCurrent {
0068 
0069 public:
0070 
0071   /**
0072    * The default constructor.
0073    */
0074   TwoKaonOnePionCurrent();
0075 
0076   /** @name Methods for the construction of the phase space integrator. */
0077   //@{
0078   /**
0079    * Complete the construction of the decay mode for integration.classes inheriting
0080    * from this one.
0081    * This method is purely virtual and must be implemented in the classes inheriting
0082    * from WeakCurrent.
0083    * @param icharge   The total charge of the outgoing particles in the current.
0084    * @param resonance If specified only include terms with this particle
0085    * @param flavour Information on the required flavours of the quarks
0086    * @param imode     The mode in the current being asked for.
0087    * @param mode      The phase space mode for the integration
0088    * @param iloc      The location of the of the first particle from the current in
0089    *                  the list of outgoing particles.
0090    * @param ires      The location of the first intermediate for the current.
0091    * @param phase     The prototype phase space channel for the integration.
0092    * @param upp       The maximum possible mass the particles in the current are
0093    *                  allowed to have.
0094    * @return Whether the current was sucessfully constructed.
0095    */
0096   virtual bool createMode(int icharge, tcPDPtr resonance,
0097               FlavourInfo flavour,
0098               unsigned int imode,PhaseSpaceModePtr mode,
0099               unsigned int iloc,int ires,
0100               PhaseSpaceChannel phase, Energy upp );
0101   //@}
0102 
0103 
0104   /**
0105    * Hadronic current. This method is purely virtual and must be implemented in
0106    * all classes inheriting from this one.
0107    * @param resonance If specified only include terms with this particle
0108    * @param flavour Information on the required flavours of the quarks
0109    * @param imode The mode
0110    * @param ichan The phase-space channel the current is needed for.
0111    * @param scale The invariant mass of the particles in the current.
0112    * @param outgoing The particles produced in the decay
0113    * @param momenta  The momenta of the particles produced in the decay
0114    * @param meopt Option for the calculation of the matrix element
0115    * @return The current. 
0116    */
0117   virtual vector<LorentzPolarizationVectorE> 
0118   current(tcPDPtr resonance,
0119       FlavourInfo flavour,
0120       const int imode, const int ichan,Energy & scale,
0121       const tPDVector & outgoing,
0122       const vector<Lorentz5Momentum> & momenta,
0123       DecayIntegrator::MEOption meopt) const;
0124 
0125   /**
0126    * Accept the decay. Checks the mesons against the list.
0127    * @param id The id's of the particles in the current.
0128    * @return Can this current have the external particles specified.
0129    */
0130   virtual bool accept(vector<int> id);
0131 
0132   /**
0133    * Return the decay mode number for a given set of particles in the current. 
0134    * Checks the mesons against the list.
0135    * @param id The id's of the particles in the current.
0136    * @return The number of the mode
0137    */
0138   virtual unsigned int decayMode(vector<int> id);
0139 
0140   /**
0141    * The particles produced by the current. This returns the mesons for the mode.
0142    * @param icharge The total charge of the particles in the current.
0143    * @param imode The mode for which the particles are being requested
0144    * @param iq The PDG code for the quark
0145    * @param ia The PDG code for the antiquark
0146    * @return The external particles for the current.
0147    */
0148   virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
0149   
0150   /**
0151    * Output the setup information for the particle database
0152    * @param os The stream to output the information to
0153    * @param header Whether or not to output the information for MySQL
0154    * @param create Whether or not to add a statement creating the object
0155    */
0156   virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
0157   
0158   /**
0159    * the matrix element for the \f$a_1\f$ decay to calculate the running width
0160    * @param imode The mode for which the matrix element is needed.
0161    * @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
0162    * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
0163    * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
0164    * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
0165    * @param m1 The mass of the first  outgoing particle.
0166    * @param m2 The mass of the second outgoing particle.
0167    * @param m3 The mass of the third  outgoing particle.
0168    * @return The matrix element squared summed over spins.
0169    */
0170   double threeBodyMatrixElement(const int imode,  const Energy2 q2,
0171                 const Energy2 s3, const Energy2 s2, 
0172                 const Energy2 s1, const Energy  m1, 
0173                 const Energy  m2, const Energy  m3) const;
0174   
0175 public:
0176 
0177   /** @name Functions used by the persistent I/O system. */
0178   //@{
0179   /**
0180    * Function used to write out object persistently.
0181    * @param os the persistent output stream written to.
0182    */
0183   void persistentOutput(PersistentOStream & os) const;
0184 
0185   /**
0186    * Function used to read in object persistently.
0187    * @param is the persistent input stream read from.
0188    * @param version the version number of the object when written.
0189    */
0190   void persistentInput(PersistentIStream & is, int version);
0191   //@}
0192 
0193   /**
0194    * The standard Init function used to initialize the interfaces.
0195    * Called exactly once for each class by the class description system
0196    * before the main function starts or
0197    * when this class is dynamically loaded.
0198    */
0199   static void Init();
0200 
0201 protected:
0202 
0203   /** @name Clone Methods. */
0204   //@{
0205   /**
0206    * Make a simple clone of this object.
0207    * @return a pointer to the new object.
0208    */
0209   virtual IBPtr clone() const;
0210 
0211   /** Make a clone of this object, possibly modifying the cloned object
0212    * to make it sane.
0213    * @return a pointer to the new object.
0214    */
0215   virtual IBPtr fullclone() const;
0216   //@}
0217 
0218 protected:
0219 
0220   /** @name Standard Interfaced functions. */
0221   //@{
0222   /**
0223    * Initialize this object after the setup phase before saving and
0224    * EventGenerator to disk.
0225    * @throws InitException if object could not be initialized properly.
0226    */
0227   virtual void doinit();
0228 
0229   /**
0230    * Initialize this object to the begining of the run phase.
0231    */
0232   virtual void doinitrun();
0233 
0234   /**
0235    * Check sanity of the object during the setup phase.
0236    */
0237   virtual void doupdate();
0238   //@}
0239 
0240 private:
0241 
0242   /**
0243    * The assignment operator is private and must never be called.
0244    * In fact, it should not even be implemented.
0245    */
0246   TwoKaonOnePionCurrent & operator=(const TwoKaonOnePionCurrent &) = delete;
0247 
0248 private:
0249 
0250   /**
0251    *  The \f$\rho\f$ lineshape for the axial-vector terms
0252    * @param q2 The scale \f$q^2\f$ for the lineshape
0253    * @param ires Which \f$\rho\f$ multiplet
0254    */
0255   Complex Trho1(Energy2 q2,int ires) const {
0256     if(ires>=int(_rho1wgts.size())) return 0.;
0257     double norm = std::accumulate(_rho1wgts.begin(),_rho1wgts.end(),0.);
0258     unsigned int imin=0,imax=_rho1wgts.size();
0259     if(ires>0) {
0260       imin=ires;
0261       imax=imin+1;
0262     }
0263     Complex output(0.);
0264     for(unsigned int ix=imin;ix<imax;++ix)
0265       output+=_rho1wgts[ix]*
0266     Resonance::BreitWignerPWave(q2,_rho1mass[ix],_rho1width[ix],_mpi,_mpi);
0267     return output/norm;
0268   }
0269 
0270   /**
0271    *  The \f$\rho\f$ lineshape for the vector terms
0272    * @param q2 The scale \f$q^2\f$ for the lineshape
0273    * @param ires Which \f$\rho\f$ multiplet
0274    */
0275   Complex Trho2(Energy2 q2,int ires) const {
0276     if(ires>=int(_rho2wgts.size())) return 0.;
0277     double norm = std::accumulate(_rho2wgts.begin(),_rho2wgts.end(),0.);
0278     unsigned int imin=0,imax=_rho2wgts.size();
0279     if(ires>0) {
0280       imin=ires;
0281       imax=imin+1;
0282     }
0283     Complex output(0.);
0284     for(unsigned int ix=imin;ix<imax;++ix)
0285       output+=_rho2wgts[ix]*
0286     Resonance::BreitWignerPWave(q2,_rho2mass[ix],_rho2width[ix],_mpi,_mpi);
0287     return output/norm;
0288   }
0289 
0290   /**
0291    *  The \f$K^*\f$ lineshape for the axial-vector terms
0292    * @param q2 The scale \f$q^2\f$ for the lineshape
0293    * @param ires Which \f$K^*\f$ multiplet
0294    */
0295   Complex TKstar1(Energy2 q2,int ires) const {
0296     if(ires>=int(_kstar1wgts.size())) return 0.;
0297     double norm = std::accumulate(_kstar1wgts.begin(),_kstar1wgts.end(),0.);
0298     unsigned int imin=0,imax=_kstar1wgts.size();
0299     if(ires>0) {
0300       imin=ires;
0301       imax=imin+1;
0302     }
0303     Complex output(0.);
0304     for(unsigned int ix=imin;ix<imax;++ix)
0305       output+=_kstar1wgts[ix]*
0306     Resonance::BreitWignerPWave(q2,_kstar1mass[ix],_kstar1width[ix],_mK,_mpi);
0307     return output/norm;
0308   }
0309 
0310   /**
0311    * \f$a_1\f$ Breit-Wigner
0312    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
0313    * @return The Breit-Wigner
0314    */
0315   Complex a1BreitWigner(Energy2 q2) const {
0316     Complex ii(0.,1.);
0317     Energy2 m2(_a1mass*_a1mass);
0318     Energy  q(sqrt(q2));
0319     Energy gamma = !_a1opt ?
0320       _a1mass*_a1width*Resonance::ga1(q2)/Resonance::ga1(_a1mass*_a1mass)/sqrt(q2) : (*_a1runinter)(q2);
0321     return m2/(m2-q2-ii*q*gamma);
0322   }
0323 
0324   /**
0325    * Initialize the \f$a_1\f$ running width
0326    * @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
0327    */
0328   void inita1Width(int iopt);
0329 
0330   /**
0331    *  The \f$T_\omega\f$ function
0332    * @param q2 The scale
0333    * @param ires the resonance
0334    */
0335   Complex Tomega(Energy2 q2, int ires) const;
0336 
0337   /**
0338    *  The \f$\omega\f$ and \f$\phi\f$ Breit-Wigner
0339    * @param q2 The scale
0340    * @param ires the resonance
0341    */
0342   Complex OmegaPhiBreitWigner(Energy2 q2, unsigned int ires) const {
0343     Energy2 m2,mg;
0344     if(ires==0) {
0345       m2=sqr(_omegamass);
0346       mg=_omegamass*_omegawidth;
0347     }
0348     else {
0349       m2=sqr(_phimass);
0350       mg=_phimass*_phiwidth;
0351     }
0352     return (-m2+Complex(0.,1.)*mg)/(q2-m2+Complex(0.,1.)*mg);
0353   }
0354 
0355   /**
0356    * The \f$\omega-\phi\f$ \f$K^*\f$ form-factor for the \f$F_5\f$ form-factor
0357    * @param s1 The scale \f$s_1\f$.
0358    * @param s2 The scale \f$s_2\f$.
0359    * @param ires Which resonances to use
0360    * @return The mixed Breit-Wigner
0361    */
0362   Complex TOmegaKStar(Energy2 s1,Energy2 s2,int ires) const;
0363 
0364 private:
0365 
0366   /**
0367    *  Parameters for the \f$\rho\f$ in the axial-vector terms
0368    */
0369   //@{
0370   /**
0371    *  Weight for the different resonances
0372    */
0373   vector<double> _rho1wgts;
0374 
0375   /**
0376    *  Masses
0377    */
0378   vector<Energy> _rho1mass;
0379 
0380   /**
0381    *  Widths
0382    */
0383   vector<Energy> _rho1width;
0384   //@}
0385 
0386   /**
0387    *  Parameters for the \f$\rho\f$ in the vector terms
0388    */
0389   //@{
0390   /**
0391    *  Weight for the different resonances
0392    */
0393   vector<double> _rho2wgts;
0394 
0395   /**
0396    *  Masses
0397    */
0398   vector<Energy> _rho2mass;
0399 
0400   /**
0401    *  Widths
0402    */
0403   vector<Energy> _rho2width;
0404   //@}
0405 
0406   /**
0407    *  Parameters for the \f$K^*\f$ in the axial-vector terms
0408    */
0409   //@{
0410   /**
0411    *  Weight for the different resonances
0412    */
0413   vector<double> _kstar1wgts;
0414 
0415   /**
0416    *  Masses
0417    */
0418   vector<Energy> _kstar1mass;
0419 
0420   /**
0421    *  Widths
0422    */
0423   vector<Energy> _kstar1width;
0424   //@}
0425 
0426   /**
0427    *  Parameters for the three meson resonances
0428    */
0429   //@{
0430   /**
0431    * The mass of the \f$a_1\f$ resonances.
0432    */
0433   Energy _a1mass;
0434   
0435   /**
0436    * The width of the \f$a_1\f$ resonances.
0437    */
0438   Energy _a1width;
0439 
0440   /**
0441    * The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
0442    */
0443   vector<Energy>  _a1runwidth;
0444 
0445   /**
0446    * The \f$q^2\f$ for the running \f$a_1\f$  width calculation.
0447    */
0448   vector<Energy2> _a1runq2;
0449 
0450   /**
0451    * The interpolator for the running \f$a_1\f$ width calculation.
0452    */
0453   Interpolator<Energy,Energy2>::Ptr _a1runinter;
0454   //@}
0455 
0456   /**
0457    *  Parameters for the \f$T_\omega\f$ function
0458    */
0459   //@{
0460   /**
0461    *  Mixing parameter
0462    */
0463   double _epsomega;
0464 
0465   /**
0466    *  Mass of the \f$\omega\f$
0467    */
0468   Energy _omegamass;
0469 
0470   /**
0471    *  Width of the \f$\omega\f$
0472    */
0473   Energy _omegawidth;
0474 
0475   /**
0476    *  Mass of the \f$\phi\f$
0477    */
0478   Energy _phimass;
0479 
0480   /**
0481    *  Width of the \f$\phi\f$
0482    */
0483   Energy _phiwidth;
0484   //@}
0485 
0486   /**
0487    * The relative weight of the \f$\omega-\phi\f$ and \f$K^*\f$ where needed.
0488    */
0489   double _omegaKstarwgt;
0490 
0491   /**
0492    * The pion decay constant, \f$f_\pi\f$.
0493    */
0494   Energy _fpi;
0495 
0496   /**
0497    * The pion mass
0498    */
0499   Energy _mpi;
0500 
0501   /**
0502    * The kaon mass
0503    */
0504   Energy _mK;
0505 
0506   /**
0507    *  Initialization switches
0508    */
0509   //@{
0510   /**
0511    * Initialize the running \f$a_1\f$ width.
0512    */
0513   bool _initializea1;
0514 
0515   /**
0516    * Option for the \f$a_1\f$ width
0517    */
0518   bool _a1opt;
0519   //@}
0520 
0521   /**
0522    *  The maximum mass of the hadronic system
0523    */
0524   Energy _maxmass;
0525 
0526   /**
0527    *  The maximum mass when the running width was calculated
0528    */
0529   Energy _maxcalc;
0530 };
0531 
0532 }
0533 
0534 #endif /* HERWIG_TwoKaonOnePionCurrent_H */