Herwig/Decay/MamboDecayer.h

0001 // -*- C++ -*-
0002 //
0003 // MamboDecayer.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_MamboDecayer_H
0010 #define HERWIG_MamboDecayer_H
0011 //
0012 // This is the declaration of the MamboDecayer class.
0013 //
0014
0015 #include "HwDecayerBase.h"
0016 #include "ThePEG/PDT/DecayMode.h"
0017
0018 namespace Herwig {
0019 using namespace ThePEG;
0020
0021   /**
0022    * The MamboDecayer class inherits from the Decayer class in
0023    * ThePEG and implements the algorithm of R.Kleiss and
0024    * W.J.Stirling NPB 385 (1992) 413-432 for massive multi-particle phase-space
0025    * decays
0026    */
0027 class MamboDecayer: public HwDecayerBase {
0028
0029 public:
0030
0031   /**
0032    * The default constructor.
0033    */
0034   MamboDecayer() : _maxweight(10.), _a0(10,0.), _a1(10,0.) {}
0035
0036   /**
0037    * Check if this decayer can perfom the decay for a particular mode
0038    * @param parent The decaying particle
0039    * @param children The decay products
0040    * @return true If this decayer can handle the given mode, otherwise false.
0041    */
0042   virtual bool accept(tcPDPtr parent, const tPDVector & children) const;
0043
0044   /**
0045    *  Perform the decay of the particle to the specified decay products
0046    * @param parent The decaying particle
0047    * @param children The decay products
0048    * @return a ParticleVector containing the decay products.
0049    */
0050   virtual ParticleVector decay(const Particle & parent,
0051                    const tPDVector & children) const;
0052
0053   /**
0054    * Output the setup information for the particle database
0055    * @param os The stream to output the information to
0056    * @param header Whether or not to output the information for MySQL
0057    */
0058   virtual void dataBaseOutput(ofstream & os,bool header) const;
0059
0060
0061 public:
0062
0063   /** @name Functions used by the persistent I/O system. */
0064   //@{
0065   /**
0066    * Function used to write out object persistently.
0067    * @param os the persistent output stream written to.
0068    */
0069   void persistentOutput(PersistentOStream & os) const;
0070
0071   /**
0072    * Function used to read in object persistently.
0073    * @param is the persistent input stream read from.
0074    * @param version the version number of the object when written.
0075    */
0076   void persistentInput(PersistentIStream & is, int version);
0077   //@}
0078
0079   /**
0080    * The standard Init function used to initialize the interfaces.
0081    * Called exactly once for each class by the class description system
0082    * before the main function starts or
0083    * when this class is dynamically loaded.
0084    */
0085   static void Init();
0086
0087 protected:
0088
0089   /** @name Clone Methods. */
0090   //@{
0091   /**
0092    * Make a simple clone of this object.
0093    * @return a pointer to the new object.
0094    */
0095   virtual IBPtr clone() const {return new_ptr(*this);}
0096
0097   /** Make a clone of this object, possibly modifying the cloned object
0098    * to make it sane.
0099    * @return a pointer to the new object.
0100    */
0101   virtual IBPtr fullclone() const {return new_ptr(*this);}
0102   //@}
0103
0104 protected:
0105
0106   /** @name Standard Interfaced functions. */
0107   //@{
0108   /**
0109    * Initialize this object. Called in the run phase just before
0110    * a run begins.
0111    */
0112   virtual void doinitrun();
0113   //@}
0114
0115 private:
0116
0117   /**
0118    * The assignment operator is private and must never be called.
0119    * In fact, it should not even be implemented.
0120    */
0121   MamboDecayer & operator=(const MamboDecayer &) = delete;
0122
0123 private:
0124
0125    /**
0126      *Set array of mometum to particles
0127      *@param mom   Momentum set to be distributed over phase-space
0128      *@param comEn The mass of the decaying particle
0129      *@return The weight of the configuration
0130      **/
0131   double calculateMomentum(vector<Lorentz5Momentum> & mom,
0132                Energy comEn) const;
0133
0134   /**
0135    * Set up the colour connections for the decay
0136    * @param parent The incoming particle
0137    * @param out The decay products
0138    */
0139   void colourConnections(const Particle & parent,
0140              ParticleVector & out) const;
0141
0142   /** @name Bessel Functions.*/
0143   //@{
0144   /**
0145    * Compute the values \f$K_0(x)/K_1(x)\f$ and it's derivative using
0146    * asymptotic expansion for large x values.
0147    * @param x The argument
0148    * @param f The value of the ratio
0149    * @param fp The value of the derivative ratio
0150    */
0151   void BesselFns(const long double x,
0152          long double & f, long double & fp) const {
0153     assert(x>=0.);
0154     if( x < 10. ) {
0155       f = BesselK0(x)/BesselK1(x);
0156       fp = ( sqr(f)*x + f - x )/x;
0157     }
0158     else
0159       BesselIExpand(-x, f, fp);
0160   }
0161
0162   /**
0163    * Compute the values \f$I_0(x)/I_1(x)\f$ and it's derivative using
0164    * asymptotic expansion.
0165    * @param x The argument
0166    * @param f The value of the ratio
0167    * @param fp The value of the derivative ratio
0168    */
0169   void BesselIExpand(const long double x,
0170              long double & f, long double & fp) const {
0171     long double y = 1./x;
0172     f = 1.+ y*(_a0[0] + y*(_a0[1] + y*(_a0[2] + y*(_a0[3]
0173         + y*(_a0[4] + y*(_a0[5] + y*(_a0[6] + y*(_a0[7]
0174         + y*(_a0[8] + y*_a0[9] )))))))));
0175     fp = -y*y*(_a1[0] + y*(_a1[1] + y*(_a1[2] + y*(_a1[3]
0176         + y*(_a1[4] + y*(_a1[5] + y*(_a1[6] + y*(_a1[7]
0177         + y*(_a1[8] + y*_a1[9] )))))))));
0178   }
0179
0180   /**
0181    * Modified Bessel function of first kind \f$I_0(x)\f$.
0182    *@param x Argument of Bessel Function
0183    **/
0184   long double BesselI0(const long double x) const {
0185     long double y,ans;
0186     if(x < 3.75) {
0187       y = sqr(x/3.75);
0188       ans = 1. + y*(3.5156229 + y*(3.0899424 + y*(1.2067492
0189           + y*(0.2659732 + y*(0.0360768+y*0.0045813)))));
0190     }
0191     else {
0192       y = (3.75/x);
0193       ans = (exp(x)/sqrt(x))*(0.39894228 + y*(0.01328592
0194           + y*(0.00225319 + y*(-0.00157565 + y*(0.00916281
0195           + y*(-0.02057706+y*(0.02635537+y*(-0.01647633+y*0.00392377))))))));
0196     }
0197     return ans;
0198   }
0199
0200   /**
0201    *  Modified Bessel function of first kind \f$I_1(x)\f$.
0202    *@param x Argument of Bessel Function
0203    **/
0204   long double BesselI1(const long double x) const {
0205     long double y,ans;
0206     if(x < 3.75) {
0207       y = sqr(x/3.75);
0208       ans = x*(0.5 + y*(0.87890594 + y*(0.51498869 + y*(0.15084934
0209           + y*(0.02658733 + y*(0.00301532 + y*0.00032411))))));
0210     }
0211     else {
0212       y = 3.75/x;
0213       ans = (0.39894228 + y*(-0.03988024 + y*(-0.00362018
0214           + y*(0.00163801 + y*(-0.01031555 + y*(0.02282967
0215       + y*(-0.02895312 + y*(0.01787654-y*0.00420059))))))))*(exp(x)/sqrt(x));
0216     }
0217     return ans;
0218   }
0219
0220   /**
0221    * Modified Bessel function of second kind \f$K_0(x)\f$.
0222    * @param x Argument of Bessel Function
0223    **/
0224   long double BesselK0(const long double x) const {
0225     long double y,ans;
0226     if(x <= 2.0) {
0227       y = x*x/4.0;
0228       ans = -log(x/2.0)*BesselI0(x) - 0.57721566
0229           + y*(0.42278420 + y*(0.23069756
0230           + y*(0.03488590 + y*(0.00262698 + y*(0.00010750+y*0.00000740)))));
0231     }
0232     else {
0233       y = 2.0/x;
0234       ans = (1.25331414 + y*(-0.07832358 + y*(+0.02189568
0235           + y*(-0.01062446 + y*(0.00587872
0236           + y*(-0.00251540 + y*0.00053208))))))*(exp(-x)/sqrt(x));
0237     }
0238     return ans;
0239   }
0240
0241   /**
0242    * Modified Bessel function of second kind \f$K_1(x)\f$.
0243    * @param x Argument of Bessel Function
0244    **/
0245   long double BesselK1(const long double x) const  {
0246     long double y,ans;
0247     if(x <= 2.0) {
0248       y = x*x/4.;
0249       ans = log(x/2.)*BesselI1(x) + (1./x)*(1. + y*(0.15443144
0250           + y*(-0.67278579 + y*(-0.18156897
0251           + y*(-0.01919402+y*(-0.00110404-(y*0.00004686)))))));
0252     }
0253     else {
0254       y = 2./x;
0255       ans = (exp(-x)/sqrt(x))*(1.25331414 + y*(0.23498619
0256           + y*(-0.03655620 + y*(0.01504268 + y*(-0.00780353
0257           + y*(0.00325614+y*(-0.00068245)))))));
0258     }
0259     return ans;
0260   }
0261   //@}
0262
0263 private:
0264
0265   /**
0266    * Maximum weight
0267    */
0268   double _maxweight;
0269
0270   /**
0271    * Store coefficents for aysymptotic expansion of \f$\frac{I_0}{I_1}\f$
0272    */
0273   vector<double> _a0;
0274
0275   /**
0276    * Store data for aysymptotic expansion of the first derivative
0277    * \f$\frac{I_0}{I_1}\f$.
0278    */
0279   vector<double> _a1;
0280 };
0281
0282 }
0283
0284 #endif /* HERWIG_MamboDecayer_H */