EIC code displayed by LXR master/​include/​Herwig/​Decay/​WeakCurrents/​TwoKaonOnePionDefaultCurrent.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

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

 // -*- C++ -*-
 //
 // TwoKaonOnePionDefaultCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
 // Copyright (C) 2002-2019 The Herwig Collaboration
 //
 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 #ifndef HERWIG_TwoKaonOnePionDefaultCurrent_H
 #define HERWIG_TwoKaonOnePionDefaultCurrent_H
 //
 // This is the declaration of the TwoKaonOnePionDefaultCurrent class.
 //
 #include "WeakCurrent.h"
 #include "Herwig/Utilities/Interpolator.h"
 #include "Herwig/Utilities/Kinematics.h"
 #include "ThePEG/StandardModel/StandardModelBase.h"
 #include "Herwig/Decay/ResonanceHelpers.h"
 #include <numeric>
 
 namespace Herwig {
 using namespace ThePEG;
 
 /** \ingroup Decay
  *
  * The TwoKaonOnePionDefaultCurrent class implements the currents from Z.Phys.C58:445 (1992),
  * this paper uses the form from Z.Phys.C48:445 (1990) for the \f$a_1\f$ width and
  * is the default model in TAUOLA.
  *
  *  The following three meson modes are implemented.
  *
  * - \f$    K^-   \pi^-    K^+ \f$, (imode=0)
  * - \f$    K^0   \pi^-    \bar{K}^0\f$, (imode=1)
  * - \f$    K^-   \pi^0    K^0 \f$, (imode=2)
  *
  *  using the currents from TAUOLA
  *
  *
  * @see WeakCurrent
  * @see Defaulta1MatrixElement
  * 
  */
 class TwoKaonOnePionDefaultCurrent: public WeakCurrent {
 
   /**
    * The matrix element for the running \f$a_1\f$ width is a friend to 
    * keep some members private.
    */
   friend class Defaulta1MatrixElement;
 
 public:
 
   /**
    * Default constructor
    */
   TwoKaonOnePionDefaultCurrent();
 
   /**
    * Hadronic current. This method is purely virtual and must be implemented in
    * all classes inheriting from this one.
    * @param resonance If specified only include terms with this particle
    * @param flavour Information on the required flavours of the quarks
    * @param imode The mode
    * @param ichan The phase-space channel the current is needed for.
    * @param scale The invariant mass of the particles in the current.
    * @param outgoing The particles produced in the decay
    * @param momenta  The momenta of the particles produced in the decay
    * @param meopt Option for the calculation of the matrix element
    * @return The current. 
    */
   virtual vector<LorentzPolarizationVectorE> 
   current(tcPDPtr resonance,
       FlavourInfo flavour,
       const int imode, const int ichan,Energy & scale,
       const tPDVector & outgoing,
       const vector<Lorentz5Momentum> & momenta,
       DecayIntegrator::MEOption meopt) const;
 
   /**
    * Accept the decay. Checks the mesons against the list.
    * @param id The id's of the particles in the current.
    * @return Can this current have the external particles specified.
    */
   virtual bool accept(vector<int> id);
 
   /**
    * Return the decay mode number for a given set of particles in the current. 
    * Checks the mesons against the list.
    * @param id The id's of the particles in the current.
    * @return The number of the mode
    */
   virtual unsigned int decayMode(vector<int> id);
 
   /**
    * The particles produced by the current. This returns the mesons for the mode.
    * @param icharge The total charge of the particles in the current.
    * @param imode The mode for which the particles are being requested
    * @param iq The PDG code for the quark
    * @param ia The PDG code for the antiquark
    * @return The external particles for the current.
    */
   virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
 
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * Standard Init function used to initialize the interfaces.
    */
   static void Init();
 
 public:
 
   /** @name Methods for the construction of the phase space integrator. */
   //@{
   /**
    * Complete the construction of the decay mode for integration.classes inheriting
    * from this one.
    * This method is purely virtual and must be implemented in the classes inheriting
    * from WeakCurrent.
    * @param icharge   The total charge of the outgoing particles in the current.
    * @param resonance If specified only include terms with this particle
    * @param flavour Information on the required flavours of the quarks
    * @param imode     The mode in the current being asked for.
    * @param mode      The phase space mode for the integration
    * @param iloc      The location of the of the first particle from the current in
    *                  the list of outgoing particles.
    * @param ires      The location of the first intermediate for the current.
    * @param phase     The prototype phase space channel for the integration.
    * @param upp       The maximum possible mass the particles in the current are
    *                  allowed to have.
    * @return Whether the current was sucessfully constructed.
    */
   virtual bool createMode(int icharge, tcPDPtr resonance,
               FlavourInfo flavour,
               unsigned int imode,PhaseSpaceModePtr mode,
               unsigned int iloc,int ires,
               PhaseSpaceChannel phase, Energy upp );
   //@}
 
   /**
    * Output the setup information for the particle database
    * @param os The stream to output the information to
    * @param header Whether or not to output the information for MySQL
    * @param create Whether or not to add a statement creating the object
    */
   virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
   
   /**
    * the matrix element for the \f$a_1\f$ decay to calculate the running width
    * @param imode The mode for which the matrix element is needed.
    * @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
    * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
    * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
    * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
    * @param m1 The mass of the first  outgoing particle.
    * @param m2 The mass of the second outgoing particle.
    * @param m3 The mass of the third  outgoing particle.
    * @return The matrix element squared summed over spins.
    */
   double threeBodyMatrixElement(const int imode,  const Energy2 q2,
                 const Energy2 s3, const Energy2 s2, 
                 const Energy2 s1, const Energy  m1, 
                 const Energy  m2, const Energy  m3) const;
 
 protected:
 
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const {return new_ptr(*this);}
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const {return new_ptr(*this);}
   //@}
 
 protected:
 
   /** @name Standard Interfaced functions. */
   //@{
   /**
    * Initialize this object after the setup phase before saving and
    * EventGenerator to disk.
    * @throws InitException if object could not be initialized properly.
    */
   virtual void doinit();
 
   /**
    * Initialize this object to the begining of the run phase.
    */
   virtual void doinitrun();
 
   /**
    * Check sanity of the object during the setup phase.
    */
   virtual void doupdate();
   //@}
 
 private:
 
   /**
    * Private and non-existent assignment operator.
    */
   TwoKaonOnePionDefaultCurrent & operator=(const TwoKaonOnePionDefaultCurrent &) = delete;
 
 private:
   
   /**
    * The \f$\rho\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
    * @param ires Which \f$\rho\f$ multiplet
    * @return The Breit-Wigner 
    */
   Complex BrhoF123(Energy2 q2,int ires) const {
     if(ires>=int(_rhoF123wgts.size())) return 0.;
     Complex output(0.);
     Complex norm = std::accumulate(_rhoF123wgts.begin(),
                    _rhoF123wgts.end(),Complex(0.));
     unsigned int imin=0,imax=_rhoF123wgts.size();
     if(ires>0) {
       imin=ires;
       imax=imin+1;
     }
     for(unsigned int ix=imin;ix<imax;++ix)
       output+=_rhoF123wgts[ix]*Resonance::BreitWignerPWave(q2,_rhoF123masses[ix],
                                _rhoF123widths[ix],_mpi,_mpi);
     return output/norm;
   }
 
   /**
    * The \f$\rho\f$ Breit-Wigner for the \f$F_5\f$ form factors.
    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
    * @param ires Which \f$\rho\f$ multiplet
    * @return The Breit-Wigner 
    */
   Complex BrhoF5(Energy2 q2,int ires) const {
     if(ires>=int(_rhoF5wgts.size())) return 0.;
     Complex output(0.);
     Complex norm = std::accumulate(_rhoF5wgts.begin(),
                    _rhoF5wgts.end(),Complex(0.));
     unsigned int imin=0,imax=_rhoF5wgts.size();
     if(ires>0) {
       imin=ires;
       imax=imin+1;
     }
     for(unsigned int ix=imin;ix<imax;++ix)
       output+=_rhoF5wgts[ix]*Resonance::BreitWignerPWave(q2,_rhoF5masses[ix],
                                _rhoF5widths[ix],_mpi,_mpi);
     return output/norm;
   }
 
   /**
    * The \f$K^*\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
    * @param ires Which \f$\rho\f$ multiplet
    * @return The Breit-Wigner 
    */
   Complex BKstarF123(Energy2 q2,int ires) const {
     if(ires>=int(_kstarF123wgts.size())) return 0.;
     Complex output(0.);
     Complex norm = std::accumulate(_kstarF123wgts.begin(),
                    _kstarF123wgts.end(),Complex(0.));
     unsigned int imin=0,imax=_kstarF123wgts.size();
     if(ires>0) {
       imin=ires;
       imax=imin+1;
     }
     assert(imax<=_kstarF123wgts.size());
     for(unsigned int ix=imin;ix<imax;++ix)
       output+=_kstarF123wgts[ix]*Resonance::BreitWignerPWave(q2,_kstarF123masses[ix],
                                  _kstarF123widths[ix],_mpi,_mK);
     return output/norm;
   }
   
   /**
    * Mixed Breit Wigner for the \f$F_5\f$ form factor
    * @param si The scale \f$s_1\f$.
    * @param sj The scale \f$s_2\f$.
    * @param ires Which resonances to use
    * @return The mixed Breit-Wigner
    */
   Complex FKrho(Energy2 si,Energy2 sj,int ires) const {
     Complex output;
     if(ires<0)
       output = _rhoKstarwgt*BKstarF123(si,-1)+BrhoF123(sj,-1);
     else if(ires%2==0)
       output= _rhoKstarwgt*BKstarF123(si,ires/2);
     else if(ires%2==1)
       output=BrhoF123(sj,ires/2);
     return output/(1.+_rhoKstarwgt);
   }
   
   /**
    * \f$a_1\f$ Breit-Wigner
    * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
    * @return The Breit-Wigner
    */
   Complex a1BreitWigner(Energy2 q2) const  {
     if(!_a1opt)
       return Resonance::BreitWignera1(q2,_a1mass,_a1width);
     Complex ii(0.,1.);
     Energy2 m2(_a1mass*_a1mass);
     Energy  q(sqrt(q2));
     Energy width = (*_a1runinter)(q2);
     return m2/(m2-q2-ii*q*width);
   }
 
   /**
    * Initialize the \f$a_1\f$ running width
    * @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
    */
   void inita1Width(int iopt);
 
 private:
   
   /**
    * Parameters for the \f$\rho\f$ Breit-Wigner in the
    * \f$F_{1,2,3}\f$ form factors.
    */
   vector<double> _rhoF123wgts;
 
   /**
    * Parameters for the \f$K^*\f$ Breit-Wigner in the
    * \f$F_{1,2,3}\f$ form factors.
    */
   vector<double> _kstarF123wgts;
   
   /**
    * Parameters for the \f$\rho\f$ Breit-Wigner in the
    * \f$F_5\f$ form factors.
    */
   vector<double> _rhoF5wgts;
   
   /**
    * The relative weight of the \f$\rho\f$ and \f$K^*\f$ where needed.
    */
   double _rhoKstarwgt;
   
   /**
    * The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
    */
   vector<Energy>  _a1runwidth;
 
   /**
    * The \f$q^2\f$ for the running \f$a_1\f$  width calculation.
    */
   vector<Energy2> _a1runq2;
 
   /**
    * The interpolator for the running \f$a_1\f$ width calculation.
    */
   Interpolator<Energy,Energy2>::Ptr _a1runinter;
 
   /**
    * Initialize the running \f$a_1\f$ width.
    */
   bool _initializea1;
   
   /**
    * The mass of the \f$a_1\f$ resonances.
    */
   Energy _a1mass;
 
   /**
    * The width of the \f$a_1\f$ resonances.
    */
   Energy _a1width;
 
   /**
    * The pion decay constant, \f$f_\pi\f$.
    */
   Energy _fpi;
 
   /**
    * The pion mass
    */
   Energy _mpi;
 
   /**
    * The kaon mass
    */
   Energy _mK;
 
   /**
    * The \f$\rho\f$ masses for the \f$F_{1,2,3}\f$ form factors.
    */
   vector<Energy> _rhoF123masses;
 
   /**
    * The \f$\rho\f$ masses for the \f$F_5\f$ form factors.
    */
   vector<Energy> _rhoF5masses;
 
   /**
    * The \f$\rho\f$ widths for the \f$F_{1,2,3}\f$ form factors.
    */
   vector<Energy> _rhoF123widths;
 
   /**
    * The \f$\rho\f$ widths for the \f$F_5\f$ form factors.
    */
   vector<Energy> _rhoF5widths;
 
   /**
    * The \f$K^*\f$ masses for the \f$F_{1,2,3}\f$ form factors.
    */
   vector<Energy> _kstarF123masses;
 
   /**
    * The \f$K^*\f$ widths for the \f$F_{1,2,3}\f$ form factors.
    */
   vector<Energy> _kstarF123widths;
 
   /**
    * Option for the \f$a_1\f$ width
    */
   bool _a1opt;
 
   /**
    *  The maximum mass of the hadronic system
    */
   Energy _maxmass;
 
   /**
    *  The maximum mass when the running width was calculated
    */
   Energy _maxcalc;
   
 };
 
 }
 
 #endif /* HERWIG_TwoKaonOnePionDefaultCurrent_H */

This page was automatically generated by the 2.3.7 LXR engine. The LXR team