Back to home page

EIC code displayed by LXR

 
 

    


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

0001 // -*- C++ -*-
0002 //
0003 // ThreePionDefaultCurrent.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_ThreePionDefaultCurrent_H
0010 #define HERWIG_ThreePionDefaultCurrent_H
0011 //
0012 // This is the declaration of the ThreePionDefaultCurrent class.
0013 //
0014 #include "WeakCurrent.h"
0015 #include "Herwig/Utilities/Interpolator.h"
0016 #include "Herwig/Utilities/Kinematics.h"
0017 #include "ThePEG/StandardModel/StandardModelBase.h"
0018 #include "Herwig/Decay/ResonanceHelpers.h"
0019 #include <numeric>
0020 
0021 namespace Herwig {
0022 using namespace ThePEG;
0023 
0024 /** \ingroup Decay
0025  *
0026  * The ThreePionDefaultCurrent class implements the currents from Z.Phys.C58:445 (1992),
0027  * this paper uses the form from Z.Phys.C48:445 (1990) for the \f$a_1\f$ width and
0028  * is the default model in TAUOLA.
0029  *
0030  *  The following three meson modes are implemented.
0031  *
0032  * - \f$    \pi^-  \pi^-    \pi^+ \f$, (imode=0)
0033  * - \f$    \pi^0  \pi^0    \pi^- \f$, (imode=1)
0034  * - \f$    \pi^+  \pi^-    \pi^0 \f$, (imode=2)
0035  *
0036  *  using the currents from TAUOLA
0037  *
0038  *
0039  * @see WeakCurrent
0040  * @see Defaulta1MatrixElement
0041  * 
0042  */
0043 class ThreePionDefaultCurrent: public WeakCurrent {
0044 
0045   /**
0046    * The matrix element for the running \f$a_1\f$ width is a friend to 
0047    * keep some members private.
0048    */
0049   friend class Defaulta1MatrixElement;
0050 
0051 public:
0052 
0053   /**
0054    * Default constructor
0055    */
0056   ThreePionDefaultCurrent();
0057 
0058   /**
0059    * Hadronic current. This method is purely virtual and must be implemented in
0060    * all classes inheriting from this one.
0061    * @param resonance If specified only include terms with this particle
0062    * @param flavour Information on the required flavours of the quarks
0063    * @param imode The mode
0064    * @param ichan The phase-space channel the current is needed for.
0065    * @param scale The invariant mass of the particles in the current.
0066    * @param outgoing The particles produced in the decay
0067    * @param momenta  The momenta of the particles produced in the decay
0068    * @param meopt Option for the calculation of the matrix element
0069    * @return The current. 
0070    */
0071   virtual vector<LorentzPolarizationVectorE> 
0072   current(tcPDPtr resonance,
0073       FlavourInfo flavour,
0074       const int imode, const int ichan,Energy & scale,
0075       const tPDVector & outgoing,
0076       const vector<Lorentz5Momentum> & momenta,
0077       DecayIntegrator::MEOption meopt) const;
0078 
0079   /**
0080    * Accept the decay. Checks the mesons against the list.
0081    * @param id The id's of the particles in the current.
0082    * @return Can this current have the external particles specified.
0083    */
0084   virtual bool accept(vector<int> id);
0085 
0086   /**
0087    * Return the decay mode number for a given set of particles in the current. 
0088    * Checks the mesons against the list.
0089    * @param id The id's of the particles in the current.
0090    * @return The number of the mode
0091    */
0092   virtual unsigned int decayMode(vector<int> id);
0093 
0094   /**
0095    * The particles produced by the current. This returns the mesons for the mode.
0096    * @param icharge The total charge of the particles in the current.
0097    * @param imode The mode for which the particles are being requested
0098    * @param iq The PDG code for the quark
0099    * @param ia The PDG code for the antiquark
0100    * @return The external particles for the current.
0101    */
0102   virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
0103 
0104 public:
0105 
0106   /** @name Functions used by the persistent I/O system. */
0107   //@{
0108   /**
0109    * Function used to write out object persistently.
0110    * @param os the persistent output stream written to.
0111    */
0112   void persistentOutput(PersistentOStream & os) const;
0113 
0114   /**
0115    * Function used to read in object persistently.
0116    * @param is the persistent input stream read from.
0117    * @param version the version number of the object when written.
0118    */
0119   void persistentInput(PersistentIStream & is, int version);
0120   //@}
0121 
0122   /**
0123    * Standard Init function used to initialize the interfaces.
0124    */
0125   static void Init();
0126 
0127 public:
0128 
0129   /** @name Methods for the construction of the phase space integrator. */
0130   //@{
0131   /**
0132    * Complete the construction of the decay mode for integration.classes inheriting
0133    * from this one.
0134    * This method is purely virtual and must be implemented in the classes inheriting
0135    * from WeakCurrent.
0136    * @param icharge   The total charge of the outgoing particles in the current.
0137    * @param resonance If specified only include terms with this particle
0138    * @param flavour Information on the required flavours of the quarks
0139    * @param imode     The mode in the current being asked for.
0140    * @param mode      The phase space mode for the integration
0141    * @param iloc      The location of the of the first particle from the current in
0142    *                  the list of outgoing particles.
0143    * @param ires      The location of the first intermediate for the current.
0144    * @param phase     The prototype phase space channel for the integration.
0145    * @param upp       The maximum possible mass the particles in the current are
0146    *                  allowed to have.
0147    * @return Whether the current was sucessfully constructed.
0148    */
0149   virtual bool createMode(int icharge, tcPDPtr resonance,
0150               FlavourInfo flavour,
0151               unsigned int imode,PhaseSpaceModePtr mode,
0152               unsigned int iloc,int ires,
0153               PhaseSpaceChannel phase, Energy upp );
0154   //@}
0155 
0156   /**
0157    * Output the setup information for the particle database
0158    * @param os The stream to output the information to
0159    * @param header Whether or not to output the information for MySQL
0160    * @param create Whether or not to add a statement creating the object
0161    */
0162   virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
0163   
0164   /**
0165    * the matrix element for the \f$a_1\f$ decay to calculate the running width
0166    * @param imode The mode for which the matrix element is needed.
0167    * @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
0168    * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
0169    * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
0170    * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
0171    * @param m1 The mass of the first  outgoing particle.
0172    * @param m2 The mass of the second outgoing particle.
0173    * @param m3 The mass of the third  outgoing particle.
0174    * @return The matrix element squared summed over spins.
0175    */
0176   double threeBodyMatrixElement(const int imode,  const Energy2 q2,
0177                 const Energy2 s3, const Energy2 s2, 
0178                 const Energy2 s1, const Energy  m1, 
0179                 const Energy  m2, const Energy  m3) const;
0180 protected:
0181 
0182   /** @name Clone Methods. */
0183   //@{
0184   /**
0185    * Make a simple clone of this object.
0186    * @return a pointer to the new object.
0187    */
0188   virtual IBPtr clone() const {return new_ptr(*this);}
0189 
0190   /** Make a clone of this object, possibly modifying the cloned object
0191    * to make it sane.
0192    * @return a pointer to the new object.
0193    */
0194   virtual IBPtr fullclone() const {return new_ptr(*this);}
0195   //@}
0196 
0197 protected:
0198 
0199   /** @name Standard Interfaced functions. */
0200   //@{
0201   /**
0202    * Initialize this object after the setup phase before saving and
0203    * EventGenerator to disk.
0204    * @throws InitException if object could not be initialized properly.
0205    */
0206   virtual void doinit();
0207 
0208   /**
0209    * Initialize this object to the begining of the run phase.
0210    */
0211   virtual void doinitrun();
0212 
0213   /**
0214    * Check sanity of the object during the setup phase.
0215    */
0216   virtual void doupdate();
0217   //@}
0218 
0219 private:
0220 
0221   /**
0222    * Private and non-existent assignment operator.
0223    */
0224   ThreePionDefaultCurrent & operator=(const ThreePionDefaultCurrent &) = delete;
0225 
0226 private:
0227   
0228   /**
0229    * The \f$\rho\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
0230    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
0231    * @param ires Which \f$\rho\f$ multiplet
0232    * @return The Breit-Wigner 
0233    */
0234   Complex BrhoF123(Energy2 q2,int ires) const {
0235     if(ires>=int(_rhoF123wgts.size())) return 0.;
0236     Complex output(0.);
0237     Complex norm = std::accumulate(_rhoF123wgts.begin(),
0238                    _rhoF123wgts.end(),Complex(0.));
0239     unsigned int imin=0,imax=_rhoF123wgts.size();
0240     if(ires>0) {
0241       imin=ires;
0242       imax=imin+1;
0243     }
0244     for(unsigned int ix=imin;ix<imax;++ix)
0245       output+=_rhoF123wgts[ix]*Resonance::BreitWignerPWave(q2,_rhoF123masses[ix],
0246                                _rhoF123widths[ix],_mpi,_mpi);
0247     return output/norm;
0248   }
0249   
0250   /**
0251    * \f$a_1\f$ Breit-Wigner
0252    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
0253    * @return The Breit-Wigner
0254    */
0255   Complex a1BreitWigner(Energy2 q2) const  {
0256     if(!_a1opt)
0257       return Resonance::BreitWignera1(q2,_a1mass,_a1width);
0258     Complex ii(0.,1.);
0259     Energy2 m2(_a1mass*_a1mass);
0260     Energy  q(sqrt(q2));
0261     Energy width = (*_a1runinter)(q2);
0262     return m2/(m2-q2-ii*q*width);
0263   }
0264 
0265   /**
0266    * Initialize the \f$a_1\f$ running width
0267    * @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
0268    */
0269   void inita1Width(int iopt);
0270 
0271 private:
0272   
0273   /**
0274    * Parameters for the \f$\rho\f$ Breit-Wigner in the
0275    * \f$F_{1,2,3}\f$ form factors.
0276    */
0277   vector<double> _rhoF123wgts;
0278   
0279   /**
0280    * The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
0281    */
0282   vector<Energy>  _a1runwidth;
0283 
0284   /**
0285    * The \f$q^2\f$ for the running \f$a_1\f$  width calculation.
0286    */
0287   vector<Energy2> _a1runq2;
0288 
0289   /**
0290    * The interpolator for the running \f$a_1\f$ width calculation.
0291    */
0292   Interpolator<Energy,Energy2>::Ptr _a1runinter;
0293 
0294   /**
0295    * Initialize the running \f$a_1\f$ width.
0296    */
0297   bool _initializea1;
0298   
0299   /**
0300    * The mass of the \f$a_1\f$ resonances.
0301    */
0302   Energy _a1mass;
0303 
0304   /**
0305    * The width of the \f$a_1\f$ resonances.
0306    */
0307   Energy _a1width;
0308 
0309   /**
0310    * The pion decay constant, \f$f_\pi\f$.
0311    */
0312   Energy _fpi;
0313 
0314   /**
0315    * The pion mass
0316    */
0317   Energy _mpi;
0318 
0319   /**
0320    * The \f$\rho\f$ masses for the \f$F_{1,2,3}\f$ form factors.
0321    */
0322   vector<Energy> _rhoF123masses;
0323 
0324   /**
0325    * The \f$\rho\f$ widths for the \f$F_{1,2,3}\f$ form factors.
0326    */
0327   vector<Energy> _rhoF123widths;
0328 
0329   /**
0330    * Option for the \f$a_1\f$ width
0331    */
0332   bool _a1opt;
0333 
0334   /**
0335    *  The maximum mass of the hadronic system
0336    */
0337   Energy _maxmass;
0338 
0339   /**
0340    *  The maximum mass when the running width was calculated
0341    */
0342   Energy _maxcalc;
0343   
0344 };
0345 
0346 }
0347 
0348 #endif /* HERWIG_ThreePionDefaultCurrent_H */