Back to home page

EIC code displayed by LXR

 
 

    


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

0001 // -*- C++ -*-
0002 //
0003 // FourPionNovosibirskCurrent.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_FourPionNovosibirskCurrent_H
0010 #define HERWIG_FourPionNovosibirskCurrent_H
0011 //
0012 // This is the declaration of the FourPionNovosibirskCurrent class.
0013 //
0014 #include "WeakCurrent.h"
0015 #include "Herwig/Utilities/Interpolator.h"
0016 #include "Herwig/Utilities/Kinematics.h"
0017 
0018 namespace Herwig {
0019 using namespace ThePEG;
0020 
0021 /** \ingroup Decay
0022  *
0023  * The <code>FourPionNovosibirskCurrent</code> class implements the decay of the weak 
0024  * current to 4 pions using the hadronic currents of
0025  * Comput. Phys. Commun. 146: 139-153, 2002,
0026  * which is a model based on the \f$e^+e^-\to4\pi\f$ data from Novosibirsk.
0027  *
0028  * It should be noted that there were a large number of mistakes in this paper which 
0029  * were corrected in hep-ph/0312240.
0030  *
0031  * @see WeakCurrent
0032  * @see FourPionDefaultMatrixElement
0033  * 
0034  * \author Peter Richardson
0035  * 
0036  */
0037 class FourPionNovosibirskCurrent: public WeakCurrent {
0038 
0039   /**
0040    * The FourPionDefaultMatrixElement class is a friend so it can perform the
0041    * integration.
0042    */
0043   friend class FourPionDefaultMatrixElement;
0044 
0045 public:
0046 
0047   /**
0048    * Default constructor
0049    */
0050   FourPionNovosibirskCurrent();
0051 
0052   /** @name Functions used by the persistent I/O system. */
0053   //@{
0054   /**
0055    * Function used to write out object persistently.
0056    * @param os the persistent output stream written to.
0057    */
0058   void persistentOutput(PersistentOStream & os) const;
0059 
0060   /**
0061    * Function used to read in object persistently.
0062    * @param is the persistent input stream read from.
0063    * @param version the version number of the object when written.
0064    */
0065   void persistentInput(PersistentIStream & is, int version);
0066   //@}
0067 
0068   /**
0069    * Standard Init function used to initialize the interfaces.
0070    */
0071   static void Init();
0072 
0073 public:
0074 
0075   /** @name Methods for the construction of the phase space integrator. */
0076   //@{
0077   /**
0078    * Complete the construction of the decay mode for integration.classes inheriting
0079    * from this one.
0080    * This method is purely virtual and must be implemented in the classes inheriting
0081    * from WeakCurrent.
0082    * @param icharge   The total charge of the outgoing particles in the current.
0083    * @param resonance If specified only include terms with this particle
0084    * @param flavour Information on the required flavours of the quarks
0085    * @param imode     The mode in the current being asked for.
0086    * @param mode      The phase space mode for the integration
0087    * @param iloc      The location of the of the first particle from the current in
0088    *                  the list of outgoing particles.
0089    * @param ires      The location of the first intermediate for the current.
0090    * @param phase     The prototype phase space channel for the integration.
0091    * @param upp       The maximum possible mass the particles in the current are
0092    *                  allowed to have.
0093    * @return Whether the current was sucessfully constructed.
0094    */
0095   virtual bool createMode(int icharge, tcPDPtr resonance,
0096               FlavourInfo flavour,
0097               unsigned int imode,PhaseSpaceModePtr mode,
0098               unsigned int iloc,int ires,
0099               PhaseSpaceChannel phase, Energy upp );
0100 
0101   /**
0102    * The particles produced by the current. This returns the four pions for the
0103    * current.
0104    * @param icharge The total charge of the particles in the current.
0105    * @param imode The mode for which the particles are being requested
0106    * @param iq The PDG code for the quark
0107    * @param ia The PDG code for the antiquark
0108    * @return The external particles for the current.
0109    */
0110   virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
0111   //@}
0112 
0113 
0114   /**
0115    * Hadronic current. This method is purely virtual and must be implemented in
0116    * all classes inheriting from this one.
0117    * @param resonance If specified only include terms with this particle
0118    * @param flavour Information on the required flavours of the quarks
0119    * @param imode The mode
0120    * @param ichan The phase-space channel the current is needed for.
0121    * @param scale The invariant mass of the particles in the current.
0122    * @param outgoing The particles produced in the decay
0123    * @param momenta  The momenta of the particles produced in the decay
0124    * @param meopt Option for the calculation of the matrix element
0125    * @return The current. 
0126    */
0127   virtual vector<LorentzPolarizationVectorE> 
0128   current(tcPDPtr resonance,
0129       FlavourInfo flavour,
0130       const int imode, const int ichan,Energy & scale,
0131       const tPDVector & outgoing,
0132       const vector<Lorentz5Momentum> & momenta,
0133       DecayIntegrator::MEOption meopt) const;
0134 
0135   /**
0136    * Accept the decay. Checks this is one of the four pion modes.
0137    * @param id The id's of the particles in the current.
0138    * @return Can this current have the external particles specified.
0139    */
0140   virtual bool accept(vector<int> id);
0141 
0142   /**
0143    * Return the decay mode number for a given set of particles in the current. 
0144    * Works out which four pion mode this is.
0145    * @param id The id's of the particles in the current.
0146    * @return The number of the mode
0147    */
0148   virtual unsigned int decayMode(vector<int> id);
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 to evaluate the \f$a_1\f$ running width.
0160    * @param iopt The mode
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 iopt, 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 protected:
0176 
0177   /** @name Clone Methods. */
0178   //@{
0179   /**
0180    * Make a simple clone of this object.
0181    * @return a pointer to the new object.
0182    */
0183   virtual IBPtr clone() const;
0184 
0185   /** Make a clone of this object, possibly modifying the cloned object
0186    * to make it sane.
0187    * @return a pointer to the new object.
0188    */
0189   virtual IBPtr fullclone() const;
0190   //@}
0191   
0192 protected:
0193 
0194   /** @name Standard Interfaced functions. */
0195   //@{
0196   /**
0197    * Initialize this object after the setup phase before saving and
0198    * EventGenerator to disk.
0199    * @throws InitException if object could not be initialized properly.
0200    */
0201   virtual void doinit();
0202 
0203   /**
0204    * Initialize this object to the begining of the run phase.
0205    */
0206   virtual void doinitrun();
0207 
0208   /**
0209    * Check sanity of the object during the setup phase.
0210    */
0211   virtual void doupdate();
0212   //@}
0213 
0214 private:
0215 
0216   /**
0217    * Private and non-existent assignment operator.
0218    */
0219   FourPionNovosibirskCurrent & operator=(const FourPionNovosibirskCurrent &) = delete;
0220       
0221 protected:
0222   
0223 
0224   /**
0225    * Initialize the \f$a_1\f$ width.
0226    * @param iopt Initialization option
0227    *  (-1 is full initialization and 0 sets up the interpolator for the running width)
0228    */
0229   void inita1width(int iopt);
0230   
0231   /**
0232    * Form foactor for the \f$a_1\f$ vertex.
0233    * @param q2 The scale \f$q^2\f$.
0234    * @return The \f$a_1\f$ form factor.
0235    */
0236   double a1FormFactor(Energy2 q2) const {
0237     return sqr((1.+_a1massolam2)/(1.+q2*_onedlam2));   
0238   }
0239 
0240   /**
0241    * Breit-Wigner for the \f$\sigma\f$ meson
0242    * @param q2 The scale \f$q^2\f$.
0243    * @param iopt The pion masses to used (0=\f$\pi^0\f$, 1=\f$\pi^+\f$)
0244    * @return The Breit-Wigner for the \f$\sigma\f$ meson
0245    */
0246   Complex sigmaBreitWigner(Energy2 q2,unsigned int iopt) const;
0247 
0248   /**
0249    * The \f$a_1\f$ breit wigner.
0250    * @param q2 The scale \f$q^2\f$.
0251    * @return The Breit-Wigner for the \f$a_1\f$.
0252    */
0253   Complex a1BreitWigner(Energy2 q2) const;
0254 
0255   /**
0256    * The Breit-Wigner for the \f$\omega\f$.
0257    * @param q2 The scale \f$q^2\f$.
0258    * @return The Breit-Wigner for the \f$\omega\f$.
0259    */
0260   Complex omegaBreitWigner(Energy2 q2) const;
0261 
0262   /**
0263    * The Breit-Wigner for the \f$\rho\f$.
0264    * @param q2 The scale \f$q^2\f$.
0265    * @return The Breit-Wigner for the \f$\rho\f$.
0266    */
0267   Complex rhoBreitWigner(Energy2 q2) const;
0268 
0269   /**
0270    * Return the \f$a_1\f$ running width.
0271    * @param q2 The scale \f$q^2\f$.
0272    * @return The running width.
0273    */
0274   Energy a1width(Energy2 q2) const {return (*_a1runinter)(q2);}
0275 
0276   /**
0277    * The \f$t_1\f$ current used in calculating the current.
0278    * @param q1 The first momentum.
0279    * @param q2 The first momentum.
0280    * @param q3 The first momentum.
0281    * @param q4 The first momentum.
0282    * @return The current \f$t_1\f$.
0283    */
0284   LorentzVector<complex<Energy5> > 
0285   t1(Lorentz5Momentum & q1,Lorentz5Momentum & q2,
0286      Lorentz5Momentum & q3,Lorentz5Momentum & q4) const;
0287 
0288   /**
0289    * The \f$t_2\f$ current used in calculating the current.
0290    * @param q1 The first momentum.
0291    * @param q2 The first momentum.
0292    * @param q3 The first momentum.
0293    * @param q4 The first momentum.
0294    * @param iopt 0 for \f$\sigma\to\pi^+\pi^-\f$ and 1 for \f$\sigma\to\pi^0\pi^0\f$
0295    * @return The current \f$t_2\f$.
0296    */
0297   LorentzVector<complex<Energy5> > 
0298   t2(Lorentz5Momentum & q1,Lorentz5Momentum & q2,
0299      Lorentz5Momentum & q3,Lorentz5Momentum & q4,
0300      unsigned int iopt) const;
0301 
0302   /**
0303    * The \f$t_3\f$ current used in calculating the current.
0304    * @param q1 The first momentum.
0305    * @param q2 The first momentum.
0306    * @param q3 The first momentum.
0307    * @param q4 The first momentum.
0308    * @return The current \f$t_3\f$.
0309    */
0310   LorentzVector<complex<Energy5> > 
0311   t3(Lorentz5Momentum & q1,Lorentz5Momentum & q2,
0312      Lorentz5Momentum & q3,Lorentz5Momentum & q4) const;
0313 
0314   /**
0315    * The G functions of hep-ph/0201149
0316    * @param q2 The scale \f$q^2\f$.
0317    * @param ichan Which of the four pion channels this is for.
0318    * @return The G function.
0319    */
0320   InvEnergy6 gFunction(Energy2 q2, int ichan) const;
0321 
0322   /**
0323    * The d parameter in \f$\rho\f$ the propagator.
0324    */
0325   Energy2 DParameter() const;
0326 
0327   /**
0328    * The \f$\frac{dh}{dq^2}\f$ function in the rho propagator evaluated at \f$q^2=m^2\f$.
0329    */
0330   double dhdq2Parameter() const;
0331 
0332   /**
0333    * The h function in the \f$\rho\f$ propagator.
0334    * @param q The scale.
0335    * @return The h function.
0336    */
0337   Energy2 hFunction(const Energy q) const;
0338 
0339 private:
0340   
0341   /**
0342    * Interpolating functions for the G functions of hep-ph/0201149
0343    */
0344   //@{
0345   /**
0346    * The interpolator for the \f$\omega\f$ current.
0347    */
0348   Interpolator<double,Energy>::Ptr _Fomega;
0349 
0350   /**
0351    * The interpolator for the three charged pion \f$a_1\f$ current. 
0352    */
0353   Interpolator<double,Energy>::Ptr _Fthreec;
0354 
0355   /**
0356    * The interpolator for the one   charged pion \f$a_1\f$ current.
0357    */
0358   Interpolator<double,Energy>::Ptr _Fonec;
0359 
0360   /**
0361    * The interpolator for the \f$\sigma\f$ current.
0362    */
0363   Interpolator<double,Energy2>::Ptr _Fsigma;
0364   //@}
0365 
0366   /**
0367    * The charged pion mass
0368    */
0369   Energy _mpic;
0370 
0371   /**
0372    * The neutral pion mass
0373    */
0374   Energy _mpi0;
0375 
0376   /**
0377    * The mass of the \f$\rho\f$ for the current.
0378    */
0379   Energy _rhomass;
0380 
0381   /**
0382    * The mass of the \f$a_1\f$ for the current.
0383    */
0384   Energy _a1mass;
0385 
0386   /**
0387    * The mass of the \f$\omega\f$ for the current.
0388    */
0389   Energy _omegamass;
0390 
0391   /**
0392    * The mass of the \f$\sigma\f$ for the current.
0393    */
0394   Energy _sigmamass;
0395 
0396   /**
0397    * The width for the \f$\rho\f$.
0398    */
0399   Energy _rhowidth;
0400 
0401   /**
0402    *  The \f$a_1\f$ width
0403    */
0404   Energy _a1width;
0405 
0406   /**
0407    *  The \f$\omega\f$ width.
0408    */
0409   Energy _omegawidth;
0410 
0411   /**
0412    *  The \f$\sigma\f$ width.
0413    */
0414   Energy _sigmawidth;
0415 
0416   /**
0417    * Mass for the intermediate in the phase-space, this is a technical parameter to
0418    * improve the phase-space integration efficiency.
0419    */
0420   Energy _intmass;
0421 
0422   /**
0423    * Width for the intermediate in the phase-space, this is a technical parameter to
0424    * improve the phase-space integration efficiency.
0425    */
0426   Energy _intwidth; 
0427 
0428   /**
0429    * The \f$z\f$ \f$\sigma\f$ coupling.
0430    */
0431   Complex _zsigma;
0432 
0433   /**
0434    * The magnitude of the \f$z\f$ \f$\sigma\f$ coupling.
0435    */
0436   double _zmag;
0437 
0438   /**
0439    * The phase of the \f$z\f$ \f$\sigma\f$ coupling.
0440    */
0441   double _zphase;
0442 
0443   /**
0444    * The mass parameter for the \f$a_1\f$ form-factor.
0445    */
0446   Energy2 _lambda2;
0447 
0448   /**
0449    * The inverse of the mass parameter for the \f$a_1\f$ form-factor.
0450    */
0451   InvEnergy2 _onedlam2;
0452 
0453   /**
0454    *  The physical \f$a_1\f$ mass divided by the mass parameter in the 
0455    *  \f$a_1\f$ form-factor.
0456    */
0457   double _a1massolam2;
0458 
0459   /**
0460    * The momentum of the  pions in on-shell \f$\sigma\f$ decay which is used
0461    * in the calculation of the running \f$\sigma\f$ width.
0462    */
0463   vector<Energy> _psigma;
0464 
0465   /**
0466    *  The charged pion mass squared.
0467    */
0468   Energy2 _mpic2;
0469 
0470   /**
0471    * The neutral pion mass squared
0472    */
0473   Energy2 _mpi02;
0474 
0475   /**
0476    *  The h function evaluated at the \f$\rho\f$ mass.
0477    */
0478   Energy2 _hm2;
0479 
0480   /**
0481    * The d parameter for the \f$\rho\f$ width.
0482    */
0483   Energy2 _rhoD;
0484 
0485   /**
0486    * The momentum of the pions produced in on-shell \f$rho\f$ decay.
0487    */
0488   Energy _prho;
0489 
0490   /**
0491    * \f$\frac{dh}{dq^2}\f$ evaluates at \f$q^2=m^2\f$ for the \f$\rho\f$.
0492    */
0493   double _dhdq2m2;
0494 
0495   /**
0496    * Magic number for the omega current.
0497    */
0498   InvEnergy _aomega;
0499 
0500   /**
0501    * Magic number for the three charged pion current.
0502    */
0503   InvEnergy _athreec;
0504 
0505   /**
0506    * Magic number for the one charged pion current
0507    */
0508   InvEnergy _aonec;
0509 
0510   /**
0511    * Magic number for the omega current.
0512    */
0513   double _bomega;
0514 
0515   /**
0516    * Magic number for the three charged pion current.
0517    */
0518   double _bthreec;
0519 
0520   /**
0521    * Magic number for the one charged pion current
0522    */
0523   double _bonec;
0524 
0525   /**
0526    * Magic number for the omega current.
0527    */
0528   double _comega;
0529 
0530   /**
0531    * Magic number for the three charged pion current.
0532    */
0533   double _cthreec;
0534 
0535   /**
0536    * Magic number for the one charged pion current
0537    */
0538   double _conec;
0539 
0540   /**
0541    * magic numbers for the running omega width
0542    */
0543   vector<double> _omegaparam;
0544 
0545   /**
0546    * whether or not to initialize the calculation of the \f$a_1\f$ width
0547    */
0548   bool _initializea1;
0549 
0550   /**
0551    * use local values of the particle masses
0552    */
0553   bool _localparameters;
0554 
0555   /**
0556    * The widths for the interpolation table for the running \f$a_1\f$ width.
0557    */
0558   vector<Energy> _a1runwidth;
0559 
0560   /**
0561    * The \f$q^2\f$ values for the interpolation table for the running \f$a_1\f$ width.
0562    */
0563   vector<Energy2> _a1runq2;
0564 
0565   /**
0566    * The interpolator for the running \f$a_1\f$ width.
0567    */
0568   Interpolator<Energy,Energy2>::Ptr _a1runinter;
0569 
0570   /**
0571    *  The maximum mass of the hadronic system
0572    */
0573   Energy _maxmass;
0574 
0575   /**
0576    *  The maximum mass when the running width was calculated
0577    */
0578   Energy _maxcalc;
0579 };
0580 
0581 }
0582 
0583 #endif /* HERWIG_FourPionNovosibirskCurrent_H */
0584