Herwig/Decay/GeneralDecayMatrixElement.h

0001 // -*- C++ -*-
0002 //
0003 // GeneralDecayMatrixElement.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_GeneralDecayMatrixElement_H
0010 #define HERWIG_GeneralDecayMatrixElement_H
0011 //
0012 // This is the declaration of the GeneralDecayMatrixElement class.
0013
0014 #include "DecayMatrixElement.h"
0015 #include "GeneralDecayMatrixElement.fh"
0016
0017 namespace Herwig {
0018
0019 using namespace ThePEG;
0020
0021
0022 /** \ingroup Helicity
0023  *  \author Peter Richardson
0024  *
0025  *  Implementation of the complex matrix element for a decay.
0026  *  An arbitary number of external particles are supported.
0027  *
0028  *  @see RhoDMatrix
0029  *  @see ProductionMatrixElement
0030  *  @see DecayVertex
0031  */
0032 class GeneralDecayMatrixElement: public DecayMatrixElement {
0033
0034 public:
0035
0036   /** @name Standard constructors and destructors. */
0037   //@{
0038   /**
0039    * Default constructor.
0040    */
0041   GeneralDecayMatrixElement() : DecayMatrixElement(999) {}
0042
0043   /**
0044    * Constructor for two body decay.
0045    * @param spinin \f$2S+1\f$ for the decaying particle
0046    * @param outspin1 \f$2S+1\f$ for the first  decay product.
0047    * @param outspin2 \f$2S+1\f$ for the second decay product.
0048    */
0049   GeneralDecayMatrixElement(PDT::Spin spinin,
0050                 PDT::Spin outspin1,PDT::Spin outspin2)
0051     : DecayMatrixElement(2,spinin) {
0052     outspin().push_back(outspin1);
0053     outspin().push_back(outspin2);
0054     setMESize();
0055   }
0056
0057   /**
0058    * Constructor for three body decay.
0059    * @param spinin \f$2S+1\f$ for the decaying particle
0060    * @param outspin1 \f$2S+1\f$ for the first  decay product.
0061    * @param outspin2 \f$2S+1\f$ for the second decay product.
0062    * @param outspin3 \f$2S+1\f$ for the third  decay product.
0063    */
0064   GeneralDecayMatrixElement(PDT::Spin spinin,PDT::Spin outspin1,
0065                 PDT::Spin outspin2,PDT::Spin outspin3)
0066     : DecayMatrixElement(3,spinin) {
0067     outspin().push_back(outspin1);
0068     outspin().push_back(outspin2);
0069     outspin().push_back(outspin3);
0070     setMESize();
0071   }
0072
0073   /**
0074    * Constructor for four body decay.
0075    * @param spinin \f$2S+1\f$ for the decaying particle
0076    * @param outspin1 \f$2S+1\f$ for the first  decay product.
0077    * @param outspin2 \f$2S+1\f$ for the second decay product.
0078    * @param outspin3 \f$2S+1\f$ for the third  decay product.
0079    * @param outspin4 \f$2S+1\f$ for the fourth decay product.
0080    */
0081   GeneralDecayMatrixElement(PDT::Spin spinin,PDT::Spin outspin1,PDT::Spin outspin2,
0082                 PDT::Spin outspin3,PDT::Spin outspin4)
0083     : DecayMatrixElement(4,spinin) {
0084     outspin().push_back(outspin1);
0085     outspin().push_back(outspin2);
0086     outspin().push_back(outspin3);
0087     outspin().push_back(outspin4);
0088     setMESize();
0089   }
0090
0091   /**
0092    * Constructor for five body decay.
0093    * @param spinin \f$2S+1\f$ for the decaying particle
0094    * @param outspin1 \f$2S+1\f$ for the first  decay product.
0095    * @param outspin2 \f$2S+1\f$ for the second decay product.
0096    * @param outspin3 \f$2S+1\f$ for the third  decay product.
0097    * @param outspin4 \f$2S+1\f$ for the fourth decay product.
0098    * @param outspin5 \f$2S+1\f$ for the fifth  decay product.
0099    */
0100   GeneralDecayMatrixElement(PDT::Spin spinin,PDT::Spin outspin1,PDT::Spin outspin2,
0101                 PDT::Spin outspin3,PDT::Spin outspin4,PDT::Spin outspin5)
0102     : DecayMatrixElement(5,spinin) {
0103     outspin().push_back(outspin1);
0104     outspin().push_back(outspin2);
0105     outspin().push_back(outspin3);
0106     outspin().push_back(outspin4);
0107     outspin().push_back(outspin5);
0108     setMESize();
0109   }
0110
0111   /**
0112    * Constructor for six body decay.
0113    * @param spinin \f$2S+1\f$ for the decaying particle
0114    * @param outspin1 \f$2S+1\f$ for the first  decay product.
0115    * @param outspin2 \f$2S+1\f$ for the second decay product.
0116    * @param outspin3 \f$2S+1\f$ for the third  decay product.
0117    * @param outspin4 \f$2S+1\f$ for the fourth decay product.
0118    * @param outspin5 \f$2S+1\f$ for the fifth  decay product.
0119    * @param outspin6 \f$2S+1\f$ for the sixth  decay product.
0120    */
0121   GeneralDecayMatrixElement(PDT::Spin spinin,PDT::Spin outspin1,PDT::Spin outspin2,
0122                 PDT::Spin outspin3,PDT::Spin outspin4,PDT::Spin outspin5,
0123                 PDT::Spin outspin6)
0124     : DecayMatrixElement(6,spinin)  {
0125     outspin().push_back(outspin1);
0126     outspin().push_back(outspin2);
0127     outspin().push_back(outspin3);
0128     outspin().push_back(outspin4);
0129     outspin().push_back(outspin5);
0130     outspin().push_back(outspin6);
0131     setMESize();
0132   }
0133
0134   /**
0135    * Constructor for arbitary body decay.
0136    * @param spinin \f$2S+1\f$ for the decaying particle
0137    * @param spinout \f$2S+1\f$ for the decay products.
0138    */
0139   GeneralDecayMatrixElement(PDT::Spin spinin,vector<PDT::Spin> spinout)
0140     : DecayMatrixElement(spinout.size(),spinin) {
0141     outspin() = spinout;
0142     setMESize();
0143   }
0144
0145   /**
0146    * Constructor for arbitary body decay.
0147    * @param extspin  \f$2S+1\f$ external particles.
0148    */
0149   GeneralDecayMatrixElement(vector<PDT::Spin> extspin)
0150     : DecayMatrixElement(int(extspin.size())-1,extspin[0]) {
0151     outspin() = vector<PDT::Spin>(extspin.begin()+1,extspin.end());
0152     setMESize();
0153   }
0154   //@}
0155
0156 public:
0157
0158   /**
0159    * Spin Density matrices
0160    */
0161   //@{
0162   /**
0163    * Calculate the decay matrix for this decay.
0164    * @param rhoout The \f$D\f$ matrix for this decay.
0165    */
0166   RhoDMatrix calculateDMatrix(const vector<RhoDMatrix> & rhoout) const;
0167
0168   /**
0169    * Calculate the \f$\rho\f$ matrix for a given outgoing particle.
0170    * @param ipart The outgoing particle the \f$\rho\f$ matrix is needed for
0171    * @param rhoin The \f$\rho\f$ matrix for the decaying particle.
0172    * @param rhoout he \f$D\f$ matrices for the other decay products.
0173    */
0174   RhoDMatrix calculateRhoMatrix(int ipart,const RhoDMatrix & rhoin,
0175                 const vector<RhoDMatrix> & rhoout) const;
0176
0177   /**
0178    * Contract the matrix element with the \f$\rho\f$ matrix of the
0179    * incoming particle. The spins of the decay products are summed over.
0180    * @param rhoin The \f$\rho\f$ matrix for the decaying particle.
0181    */
0182   Complex contract(const RhoDMatrix & rhoin) const;
0183
0184   /**
0185    * Contract the matrix element with the \f$\rho\f$ matrix of the
0186    * incoming particle. The spins of the decay products are summed over.
0187    * @param con The conjugate matrix elemetn for the contraction
0188    * @param rhoin The \f$\rho\f$ matrix for the decaying particle.
0189    */
0190   Complex contract(const GeneralDecayMatrixElement & con,
0191            const RhoDMatrix & rhoin);
0192   //@}
0193
0194 public:
0195
0196   /**
0197    * Access to the individual helicity components.
0198    */
0199   //@{
0200   /**
0201    * Get the helicity components for a two body decay
0202    * @param inhel The helicity of the decaying particle.
0203    * @param outhel1 The helicity of the first  decay product.
0204    * @param outhel2 The helicity of the second decay product.
0205    */
0206   Complex   operator () (unsigned int inhel,unsigned int outhel1,
0207              unsigned int outhel2) const {
0208     unsigned int iloc = inhel*constants_[1]+
0209       outhel1*constants_[2]+outhel2*constants_[3];
0210     assert(outspin().size()==2&&iloc<matrixElement_.size());
0211     return matrixElement_[iloc];
0212   }
0213
0214   /**
0215    * Set the helicity components for a two body decay
0216    * @param inhel The helicity of the decaying particle.
0217    * @param outhel1 The helicity of the first  decay product.
0218    * @param outhel2 The helicity of the second decay product.
0219    */
0220   Complex & operator () (unsigned int inhel,unsigned int outhel1,
0221              unsigned int outhel2) {
0222     unsigned int iloc = inhel*constants_[1]+
0223       outhel1*constants_[2]+outhel2*constants_[3];
0224     assert(outspin().size()==2&&iloc<matrixElement_.size());
0225     return matrixElement_[iloc];
0226   }
0227
0228   /**
0229    * Get the helicity components for a three body decay
0230    * @param inhel The helicity of the decaying particle.
0231    * @param outhel1 The helicity of the first  decay product.
0232    * @param outhel2 The helicity of the second decay product.
0233    * @param outhel3 The helicity of the third  decay product.
0234    */
0235   Complex   operator () (unsigned int inhel,unsigned int outhel1,
0236              unsigned int outhel2,unsigned int outhel3) const {
0237     unsigned int iloc = inhel*constants_[1]+outhel1*constants_[2]+
0238       outhel2*constants_[3]+outhel3*constants_[4];
0239     assert(outspin().size()==3&&iloc<matrixElement_.size());
0240     return matrixElement_[iloc];
0241   }
0242
0243   /**
0244    * Set the helicity components for a three body decay
0245    * @param inhel The helicity of the decaying particle.
0246    * @param outhel1 The helicity of the first  decay product.
0247    * @param outhel2 The helicity of the second decay product.
0248    * @param outhel3 The helicity of the third  decay product.
0249    */
0250   Complex & operator () (unsigned int inhel,unsigned int outhel1,
0251              unsigned int outhel2,unsigned int outhel3) {
0252     unsigned int iloc = inhel*constants_[1]+outhel1*constants_[2]+
0253       outhel2*constants_[3]+outhel3*constants_[4];
0254     assert(outspin().size()==3&&iloc<matrixElement_.size());
0255     return matrixElement_[iloc];
0256   }
0257
0258   /**
0259    * Get the helicity components for a four body decay
0260    * @param inhel The helicity of the decaying particle.
0261    * @param outhel1 The helicity of the first  decay product.
0262    * @param outhel2 The helicity of the second decay product.
0263    * @param outhel3 The helicity of the third  decay product.
0264    * @param outhel4 The helicity of the fourth decay product.
0265    */
0266   Complex   operator () (unsigned int inhel,unsigned int outhel1,
0267              unsigned int outhel2,unsigned int outhel3,
0268              unsigned int outhel4) const {
0269     vector<unsigned int> itemp(5);
0270     itemp[0]=inhel  ; itemp[1]=outhel1;
0271     itemp[2]=outhel2; itemp[3]=outhel3;
0272     itemp[4]=outhel4 ;
0273     return (*this)(itemp);
0274   }
0275
0276   /**
0277    * Set the helicity components for a four body decay
0278    * @param inhel The helicity of the decaying particle.
0279    * @param outhel1 The helicity of the first  decay product.
0280    * @param outhel2 The helicity of the second decay product.
0281    * @param outhel3 The helicity of the third  decay product.
0282    * @param outhel4 The helicity of the fourth decay product.
0283    */
0284   Complex & operator () (unsigned int inhel,unsigned int outhel1,
0285              unsigned int outhel2,unsigned int outhel3,
0286              unsigned int outhel4) {
0287     vector<unsigned int> itemp(5);
0288     itemp[0]=inhel  ; itemp[1]=outhel1;
0289     itemp[2]=outhel2; itemp[3]=outhel3;
0290     itemp[4]=outhel4;
0291     return (*this)(itemp);
0292   }
0293
0294   /**
0295    * Get the helicity components for a five body decay
0296    * @param inhel The helicity of the decaying particle.
0297    * @param outhel1 The helicity of the first  decay product.
0298    * @param outhel2 The helicity of the second decay product.
0299    * @param outhel3 The helicity of the third  decay product.
0300    * @param outhel4 The helicity of the fourth decay product.
0301    * @param outhel5 The helicity of the fifth  decay product.
0302    */
0303   Complex   operator () (unsigned int inhel,unsigned int outhel1,
0304              unsigned int outhel2,unsigned int outhel3,
0305              unsigned int outhel4,unsigned int outhel5) const {
0306     vector<unsigned int> itemp(6);
0307     itemp[0]=inhel  ; itemp[1]=outhel1;
0308     itemp[2]=outhel2; itemp[3]=outhel3;
0309     itemp[4]=outhel4 ;itemp[5]=outhel5;
0310     return (*this)(itemp);
0311   }
0312
0313   /**
0314    * Set the helicity components for a five body decay
0315    * @param inhel The helicity of the decaying particle.
0316    * @param outhel1 The helicity of the first  decay product.
0317    * @param outhel2 The helicity of the second decay product.
0318    * @param outhel3 The helicity of the third  decay product.
0319    * @param outhel4 The helicity of the fourth decay product.
0320    * @param outhel5 The helicity of the fifth  decay product.
0321    */
0322   Complex & operator () (unsigned int inhel,unsigned int outhel1,
0323              unsigned int outhel2,unsigned int outhel3,
0324              unsigned int outhel4,unsigned int outhel5) {
0325     vector<unsigned int> itemp(6);
0326     itemp[0]=inhel  ; itemp[1]=outhel1;
0327     itemp[2]=outhel2; itemp[3]=outhel3;
0328     itemp[4]=outhel4; itemp[5]=outhel5;
0329     return (*this)(itemp);
0330   }
0331
0332   /**
0333    * Get the helicity components for a six body decay
0334    * @param inhel The helicity of the decaying particle.
0335    * @param outhel1 The helicity of the first  decay product.
0336    * @param outhel2 The helicity of the second decay product.
0337    * @param outhel3 The helicity of the third  decay product.
0338    * @param outhel4 The helicity of the fourth decay product.
0339    * @param outhel5 The helicity of the fifth  decay product.
0340    * @param outhel6 The helicity of the sixth  decay product.
0341    */
0342   Complex   operator () (unsigned int inhel,unsigned int outhel1,
0343              unsigned int outhel2,unsigned int outhel3,
0344              unsigned int outhel4,unsigned int outhel5,
0345              unsigned int outhel6) const {
0346     vector<unsigned int> itemp(7);
0347     itemp[0]=inhel  ; itemp[1]=outhel1;
0348     itemp[2]=outhel2; itemp[3]=outhel3;
0349     itemp[4]=outhel4 ;itemp[5]=outhel5;
0350     itemp[6]=outhel6;
0351     return (*this)(itemp);
0352   }
0353
0354   /**
0355    * Set the helicity components for a six body decay
0356    * @param inhel The helicity of the decaying particle.
0357    * @param outhel1 The helicity of the first  decay product.
0358    * @param outhel2 The helicity of the second decay product.
0359    * @param outhel3 The helicity of the third  decay product.
0360    * @param outhel4 The helicity of the fourth decay product.
0361    * @param outhel5 The helicity of the fifth  decay product.
0362    * @param outhel6 The helicity of the sixth  decay product.
0363    */
0364   Complex & operator () (unsigned int inhel,unsigned int outhel1,
0365              unsigned int outhel2,unsigned int outhel3,
0366              unsigned int outhel4,unsigned int outhel5,
0367              unsigned int outhel6) {
0368     vector<unsigned int> itemp(7);
0369     itemp[0]=inhel  ; itemp[1]=outhel1;
0370     itemp[2]=outhel2; itemp[3]=outhel3;
0371     itemp[4]=outhel4; itemp[5]=outhel5;
0372     itemp[6]=outhel6;
0373     return (*this)(itemp);
0374   }
0375
0376   /**
0377    * Get the helicity components for an \f$n\f$-body decay.
0378    * @param in The helicities of the external particles.
0379    */
0380   Complex   operator () (const vector<unsigned int> & in) const {
0381     assert(in.size()==outspin().size()+1);
0382     // contribution  for the incoming particle
0383     unsigned int iloc(in[0]*constants_[1]);
0384     // contributions for the outgoing particles
0385     for(unsigned int ix=1;ix<in.size();++ix){iloc+=in[ix]*constants_[ix+1];}
0386     assert(iloc<matrixElement_.size());
0387     return matrixElement_[iloc];
0388   }
0389
0390   /**
0391    * Set the helicity components for an \f$n\f$-body decay.
0392    * @param in The helicities of the external particles.
0393    */
0394   Complex & operator () (const vector<unsigned int> & in) {
0395     assert(in.size()==outspin().size()+1);
0396     // contribution  for the incoming particle
0397     unsigned int iloc(in[0]*constants_[1]);
0398     // contributions for the outgoing particles
0399     for(unsigned int ix=1;ix<in.size();++ix){iloc+=in[ix]*constants_[ix+1];}
0400     assert(iloc<matrixElement_.size());
0401     return matrixElement_[iloc];
0402   }
0403   //@}
0404
0405   /**
0406    *  Member to zero all the elements for the matrix element
0407    */
0408   void zero() {
0409     for(unsigned int ix=0;ix<matrixElement_.size();++ix)
0410       matrixElement_[ix]=0.;
0411   }
0412
0413   /**
0414    * Set the size of the vector containing the matrix element.
0415    */
0416   void setMESize() {
0417     int isize = inspin();
0418     for(unsigned int ix=0;ix<outspin().size();++ix) isize*=outspin()[ix];
0419     matrixElement_.resize(isize,0.);
0420     // set up the constants for the mapping of helicity to vector index
0421     constants_.resize(outspin().size()+2);
0422     int temp=1;
0423     for(unsigned int ix=outspin().size();ix>0;--ix) {
0424       temp*=outspin()[ix-1];constants_[ix]=temp;
0425     }
0426     temp *= inspin();
0427     constants_[0]=temp;
0428     constants_[outspin().size()+1]=1;
0429   }
0430
0431 private:
0432
0433   /**
0434    * Storage of the matrix element, a vector is better for memory usage.
0435    */
0436   mutable vector<Complex> matrixElement_;
0437
0438   /**
0439    * Constants needed to map the index of the vector to a helicity structure.
0440    */
0441   mutable vector<unsigned int> constants_;
0442
0443 };
0444
0445 }
0446
0447 #endif /* HERWIG_GeneralDecayMatrixElement_H */