Back to home page

EIC code displayed by LXR

 
 

    


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

0001 // -*- C++ -*-
0002 //
0003 // VertexBase.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
0005 //
0006 // ThePEG 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 ThePEG_VertexBase_H
0010 #define ThePEG_VertexBase_H
0011 //
0012 // This is the declaration of the VertexBase class.
0013 
0014 #include <ThePEG/Interface/Interfaced.h>
0015 #include <ThePEG/PDT/ParticleData.h>
0016 #include <ThePEG/Helicity/HelicityDefinitions.h>
0017 #include <ThePEG/Repository/EventGenerator.h>
0018 #include "ThePEG/StandardModel/StandardModelBase.h"
0019 #include "VertexBase.fh"
0020 #include <array>
0021 
0022 
0023 namespace ThePEG {
0024 namespace Helicity {
0025 
0026 /**
0027  * Namespace for naming of vertices. Each child class should extend this
0028  * with its own spin configuration.
0029  */
0030 namespace VertexType {
0031   typedef unsigned T;
0032   /**
0033    *  Undefined Enum for the Lorentz structures
0034    */
0035   const T UNDEFINED = 0;
0036 }
0037 
0038 /**
0039  * Namespace for naming types of colour structures to allow models to define new type
0040  */
0041 namespace ColourStructure {
0042   typedef unsigned T;
0043   const T UNDEFINED  = 0;
0044   const T SINGLET    = 1;
0045   const T SU3TFUND   = 2;
0046   const T SU3F       = 3;
0047   const T SU3T6      = 4;
0048   const T SU3K6      = 5;
0049   const T EPS        = 6;
0050   const T DELTA      = 7;
0051   const T SU3FF      = 8;
0052   const T SU3TTFUNDS = 9;
0053   const T SU3TTFUNDD = 10;
0054   const T SU3TT6     = 11;
0055   const T SU3I12I34  = 12;
0056   const T SU3I14I23  = 13;
0057   const T SU3T21T43  = 14;
0058   const T SU3T23T41  = 15;
0059 }
0060 
0061 /**
0062  * Namespace for naming types of couplings to allow models to define new type
0063  */
0064 namespace CouplingType {
0065   typedef unsigned T;
0066   const T UNDEFINED = 0;
0067   const T QED = 1;
0068   const T QCD = 2;
0069 }
0070 /** \ingroup Helicity
0071  * 
0072  *  The VertexBase class is the base class for all helicity amplitude
0073  *  vertices. In implements the storage of the particles 
0074  *  which are allowed to interact at the vertex and some simple functions 
0075  *  which are often needed by the classes which implement the specific 
0076  *  vertices.
0077  *
0078  *  In practice little use is made of this information and it is mainly
0079  *  included for future extensions. It can also be used at the development
0080  *  and debugging stage.
0081  *
0082  */
0083 class VertexBase  : public Interfaced {
0084 /**
0085  *  The output operator is a friend to avoid the data being public.
0086  */
0087 friend ostream & operator<<(ostream &, const VertexBase &);
0088 
0089 public:
0090   
0091   /** @name Standard constructors and destructors. */
0092   //@{
0093   /**
0094    * Constructor for \f$n\f$-point vertices.
0095    * @param name The type of vertex
0096    * @param kine Whether the kinematic invariants should be calculated.
0097    */
0098   VertexBase(VertexType::T name, bool kine=false);
0099   //@}
0100 
0101 public:
0102  
0103   /** @name Functions used by the persistent I/O system. */
0104   //@{
0105   /**
0106    * Function used to write out object persistently.
0107    * @param os the persistent output stream written to.
0108    */
0109   void persistentOutput(PersistentOStream & os) const;
0110 
0111   /**
0112    * Function used to read in object persistently.
0113    * @param is the persistent input stream read from.
0114    * @param version the version number of the object when written.
0115    */
0116   void persistentInput(PersistentIStream & is, int version);
0117   //@}
0118   
0119   /**
0120    * Standard Init function used to initialize the interfaces.
0121    */
0122   static void Init();
0123   
0124 public:
0125 
0126   /**
0127    *  Access to the particle information
0128    */
0129   //@{
0130   /**
0131    * Number of different particle combinations allowed.
0132    */
0133   unsigned int size() const { return _particles.size(); }
0134 
0135 public:
0136   /**
0137    * Is a particle allowed as an incoming particle?
0138    * @param p The ParticleData pointer
0139    */
0140   bool isIncoming(tPDPtr p) const {
0141     return _inpart.find(p) != _inpart.end();
0142   }
0143 
0144   /**
0145    * Is a particle allowed as an outgoing particle?
0146    * @param p The ParticleData pointer
0147    */
0148   bool isOutgoing(tPDPtr p) const {
0149     return _outpart.find(p) != _outpart.end();
0150   }
0151 
0152   /**
0153    * Get the list of incoming particles.
0154    */
0155   const set<tPDPtr> & incoming() const { return _inpart; }
0156 
0157   /**
0158    * Get the list of outgoing particles.
0159    */
0160   const set<tPDPtr> & outgoing() const { return _outpart; }
0161 
0162   /**
0163    * Get the coupling.
0164    */
0165   Complex norm() const { return _norm; }
0166 
0167   /**
0168    * Function to search the list.
0169    * @param ilist Which list to search
0170    * @param id The PDG code to look for.
0171    */
0172   vector<long> search(unsigned int ilist,long id) const;
0173 
0174   /**
0175    * Function to search the list.
0176    * @param ilist Which list to search
0177    * @param id The particle to look for.
0178    */
0179   vector<tPDPtr> search(unsigned int ilist,tcPDPtr id) const;
0180 
0181   /**
0182    * Is a given combination allowed.
0183    * @param id1 PDG code of the first particle.
0184    * @param id2 PDG code of the second particle.
0185    * @param id3 PDG code of the third particle.
0186    * @param id4 PDG code of the fourth particle.
0187    */
0188   bool allowed(long id1, long id2, long id3, long id4 = 0) const;
0189 
0190   /**
0191    * Get name of Vertex
0192    */
0193   VertexType::T getName() const { return _theName; }
0194 
0195   /**
0196    * Get number of lines on Vertex
0197    */
0198   unsigned int getNpoint() const { return _npoint; }
0199 
0200   /**
0201    * Get the order in \f$g_EM\f$
0202    */
0203   int orderInGem() const { return couplingOrders_.at(CouplingType::QED); }
0204 
0205   /**
0206    * Get the order in \f$g_s\f$
0207    */
0208   int orderInGs() const { return couplingOrders_.at(CouplingType::QCD); }
0209 
0210   /**
0211    *  Get the order in a specific coupling
0212    */
0213   int orderInCoupling(CouplingType::T cType) const {
0214     if(couplingOrders_.find(cType) !=couplingOrders_.end())
0215       return couplingOrders_.at(cType);
0216     else
0217       return 0;
0218   }
0219 
0220   /**
0221    *  Get the total order of the vertex
0222    */
0223   int orderInAllCouplings() const {
0224     int output(0);
0225     for(auto & p : couplingOrders_)
0226       output += p.second;
0227     return output;
0228   }
0229 
0230   /**
0231    *  Get the colour structure
0232    */
0233   ColourStructure::T colourStructure() const {return colourStructure_;}
0234   //@}
0235 
0236 public:
0237 
0238   /**
0239    * @name Calculation of the strong, electromagnetic and weak couplings
0240    */
0241   //@{
0242   /**
0243    *  Strong coupling
0244    */
0245   double strongCoupling(Energy2 q2) const {
0246     if(_coupopt==0) {
0247       double val = 4.0*Constants::pi*generator()->standardModel()->alphaS(q2);
0248       assert(val>=0.);
0249       return sqrt(val);
0250     }
0251     else if(_coupopt==1)
0252       return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaS());
0253     else
0254       return _gs;
0255   }
0256 
0257   /**
0258    *  Electromagnetic coupling
0259    */
0260   double electroMagneticCoupling(Energy2 q2) const {
0261     if(_coupopt==0)
0262       return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2));
0263     else if(_coupopt==1)
0264       return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ());
0265     else
0266       return _ee;
0267   }
0268 
0269   /**
0270    *  Weak coupling
0271    */
0272   double weakCoupling(Energy2 q2) const {
0273     if( _coupopt == 0 )
0274       return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2)/
0275           generator()->standardModel()->sin2ThetaW());
0276     else if( _coupopt == 1 )
0277       return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ()/
0278           generator()->standardModel()->sin2ThetaW());
0279     else
0280       return _ee/_sw;
0281   }
0282 
0283   double sin2ThetaW() const {
0284     if( _coupopt == 0 || _coupopt  == 1)
0285       return generator()->standardModel()->sin2ThetaW();
0286     else
0287       return sqr(_sw);
0288   }
0289   //@}
0290 
0291 public:
0292 
0293   /**
0294    *   Set coupling methods
0295    */
0296   //@{
0297   /**
0298    * Calculate the couplings for a three point interaction.
0299    * This method is virtual and must be implemented in 
0300    * classes inheriting from this.
0301    * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
0302    * @param part1 The ParticleData pointer for the first  particle.
0303    * @param part2 The ParticleData pointer for the second particle.
0304    * @param part3 The ParticleData pointer for the third  particle.
0305    */
0306   virtual void setCoupling(Energy2 q2,tcPDPtr part1,
0307                tcPDPtr part2,tcPDPtr part3)=0;
0308 
0309   /**
0310    * Calculate the couplings for a four point interaction.
0311    * This method is virtual and must be implemented in 
0312    * classes inheriting from this.
0313    * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
0314    * @param part1 The ParticleData pointer for the first  particle.
0315    * @param part2 The ParticleData pointer for the second particle.
0316    * @param part3 The ParticleData pointer for the third  particle.
0317    * @param part4 The ParticleData pointer for the fourth particle.
0318    */
0319   virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3,
0320                tcPDPtr part4)=0;
0321   //@}
0322 
0323 protected:
0324 
0325   /** @name Standard Interfaced functions. */
0326   //@{
0327   /**
0328    * Initialize this object after the setup phase before saving an
0329    * EventGenerator to disk.
0330    * @throws InitException if object could not be initialized properly.
0331    */
0332   virtual void doinit();
0333 
0334   /**
0335    * Rebind pointer to other Interfaced objects. Called in the setup phase
0336    * after all objects used in an EventGenerator has been cloned so that
0337    * the pointers will refer to the cloned objects afterwards.
0338    * @param trans a TranslationMap relating the original objects to
0339    * their respective clones.
0340    * @throws RebindException if no cloned object was found for a given
0341    * pointer.
0342    */
0343   virtual void rebind(const TranslationMap & trans);
0344 
0345   /**
0346    * Return a vector of all pointers to Interfaced objects used in this
0347    * object.
0348    * @return a vector of pointers.
0349    */
0350   virtual IVector getReferences();
0351   //@}
0352 
0353 protected:
0354   /**
0355    *  Members to set-up the particles
0356    */
0357   //@{
0358   /**
0359    * Set up the lists of outer particles for the vertex.
0360    * @param ids A vector of PDG codes for the particles.
0361    */
0362   void addToList(const vector<long> & ids);
0363 
0364   /**
0365    * Set up the lists of outer particles for the three-/four-point vertex.
0366    * For small vertices, this form is much easier to use.
0367    * @param ida The PDG codes for the first  set of particles.
0368    * @param idb The PDG codes for the second set of particles.
0369    * @param idc The PDG codes for the third  set of particles.
0370    * @param idd The PDG codes for the fourth set of particles.
0371    */
0372   void addToList(long ida, long idb, long idc, long idd = 0);
0373   //@}
0374 
0375 protected:
0376   /**
0377    *  Members for the amplitude calculations
0378    */
0379   //@{
0380   /**
0381    * Set the coupling.
0382    * @param coup The coupling.
0383    */
0384   void norm(const Complex & coup) { _norm = coup; }
0385 
0386   /**
0387    * Calculate the propagator for a diagram.
0388    * @param iopt The option for the Breit-Wigner shape
0389    * @param q2 The scale
0390    * @param part The ParticleData pointer for the off-shell particle.
0391    * @param mass The mass if not to be taken from the ParticleData object
0392    * @param width The width if not to be taken from the ParticleData object
0393    */
0394   virtual Complex propagator(int iopt, Energy2 q2,tcPDPtr part,
0395                  complex<Energy> mass=-GeV,
0396                  complex<Energy> width=-GeV);
0397 
0398   /**
0399    * Calculate propagator multiplied by coupling.
0400    * @param iopt The option for the Breit-Wigner shape
0401    * @param q2 The scale
0402    * @param part The ParticleData pointer for the off-shell particle.
0403    * @param mass The mass if not to be taken from the ParticleData object
0404    * @param width The width if not to be taken from the ParticleData object
0405    */
0406   Complex normPropagator(int iopt, Energy2 q2,tcPDPtr part,
0407              complex<Energy> mass=-GeV, 
0408              complex<Energy> width=-GeV) {
0409     return _norm*propagator(iopt,q2,part,mass,width);
0410   }
0411   //@}    
0412 
0413 public:
0414   /** @name Kinematic invariants for loop diagrams */
0415   //@{
0416 
0417   /**
0418    * Whether or not to calculate the kinematics invariants
0419    */
0420   bool kinematics() const { return _calckinematics; }
0421 
0422   /**
0423    * Set whether or not to calculate the kinematics invariants
0424    */
0425   void kinematics(bool kine ) { _calckinematics=kine; }
0426 
0427   /**
0428    *  Calculate the kinematics for a 3-point vertex
0429    */
0430   void calculateKinematics(const Lorentz5Momentum & p0,
0431                const Lorentz5Momentum & p1,
0432                const Lorentz5Momentum & p2) {
0433     _kine[0][0]=p0*p0;
0434     _kine[1][1]=p1*p1;
0435     _kine[2][2]=p2*p2;
0436     _kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
0437     _kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
0438     _kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
0439   }
0440   
0441   /**
0442    *  Calculate the kinematics for a 4-point vertex
0443    */
0444   void calculateKinematics(const Lorentz5Momentum & p0,
0445                const Lorentz5Momentum & p1,
0446                const Lorentz5Momentum & p2,
0447                const Lorentz5Momentum & p3) {
0448     _kine[0][0]=p0*p0;
0449     _kine[1][1]=p1*p1;
0450     _kine[2][2]=p2*p2;
0451     _kine[3][3]=p3*p3;
0452     _kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
0453     _kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
0454     _kine[0][3]=p0*p3;_kine[3][0]=_kine[0][3];
0455     _kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
0456     _kine[1][3]=p1*p3;_kine[3][1]=_kine[1][3];
0457     _kine[2][3]=p2*p3;_kine[3][2]=_kine[2][3];
0458   }
0459   
0460   /**
0461    *  Calculate the kinematics for a n-point vertex
0462    */
0463   void calculateKinematics(const vector<Lorentz5Momentum> & p) {
0464     for(size_t ix=0;ix<p.size();++ix) {
0465       for(size_t iy=0;iy<=ix;++ix) {
0466     _kine[ix][iy]=p[ix]*p[iy];
0467     _kine[iy][ix]=_kine[ix][iy];
0468       }
0469     }
0470   }
0471 
0472   /**
0473    * Get one of the kinematic invariants
0474    */
0475   Energy2 invariant(unsigned int ix ,unsigned int iy) const {
0476     assert ( ix < _npoint && iy < _npoint );
0477     return _kine[ix][iy];
0478   }
0479   //@}
0480   
0481 protected:
0482 
0483   /**
0484    * Set the order in \f$g_EM\f$
0485    * @param order The order of the vertex in \f$g_EM\f$
0486    */
0487   void orderInGem(int order) { couplingOrders_[CouplingType::QED] = order; }
0488 
0489   /**
0490    * Set the order in \f$g_s\f$
0491    * @param order The order of the vertex in \f$g_s\f$
0492    */
0493   void orderInGs (int order) { couplingOrders_[CouplingType::QCD] = order; }
0494 
0495   /**
0496    *  Set the order in a specifc type of coupling
0497    */
0498   void orderInCoupling(CouplingType::T cType, int order) {
0499     couplingOrders_[cType] = order;
0500   }
0501 
0502   /**
0503    *  Set the colour structure
0504    */
0505   void colourStructure(ColourStructure::T structure) {
0506     colourStructure_ = structure;
0507   }
0508   
0509 private:
0510   
0511   /**
0512    * Private and non-existent assignment operator.
0513    */
0514   VertexBase & operator=(const VertexBase &) = delete;
0515   
0516 private:
0517 
0518   /**
0519    * Storage of the particles.   
0520    */
0521   //@{
0522   /**
0523    *  Particles interacting at the vertex
0524    */
0525   vector<vector<PDPtr> > _particles;
0526 
0527   /**
0528    *  Number of particles at the vertex
0529    */
0530   unsigned int _npoint;
0531 
0532   /**
0533    * ParticleData pointers for the allowed incoming particles.
0534    */
0535   set<tPDPtr> _inpart;
0536 
0537   /**
0538    * ParticleData pointers for the allowed outgoing particles.
0539    */
0540   set<tPDPtr> _outpart;
0541   //@}
0542 
0543   /**
0544    * The overall coupling.
0545    */
0546   Complex _norm;
0547 
0548   /**
0549    * Whether or not to calculate the kinematic invariants for the vertex
0550    */
0551   bool _calckinematics;
0552 
0553   /**
0554    * Kinematica quantities needed for loop vertices
0555    */
0556   std::array<std::array<Energy2,5>,5> _kine;
0557 
0558   /**
0559    * Name of vertex
0560    */
0561   VertexType::T _theName;
0562 
0563   /**
0564    * Colour structure of the vertex
0565    */
0566   ColourStructure::T colourStructure_;
0567 
0568   /**
0569    *  The order of the vertex in specific couplings
0570    */
0571   map<CouplingType::T,int> couplingOrders_;
0572 
0573   /**
0574    *  option for the coupling
0575    */
0576   unsigned int _coupopt;
0577 
0578   /**
0579    *  Fixed value of strong coupling to use
0580    */
0581   double _gs;
0582 
0583   /**
0584    *  Fixed value of the electromagentic coupling to use
0585    */
0586   double _ee;
0587 
0588   /**
0589    *  Fixed value of \f$\sin\theta_W\f$ to use
0590    */
0591   double _sw;
0592 };
0593   
0594 /**
0595  * Output the information on the vertex.
0596  */
0597 ostream & operator<<(ostream &, const VertexBase &);
0598   
0599 }
0600 }
0601 
0602 #endif /* ThePEG_VertexBase_H */