Herwig/Utilities/Kinematics.h

0001 // -*- C++ -*-
0002 //
0003 // Kinematics.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 //
0010 #ifndef HERWIG_Kinematics_H
0011 #define HERWIG_Kinematics_H
0012
0013 // This is the declaration of the Kinematics class.
0014
0015 #include <ThePEG/Config/ThePEG.h>
0016 #include "ThePEG/Vectors/ThreeVector.h"
0017 #include "ThePEG/Vectors/LorentzRotation.h"
0018 #include "ThePEG/Repository/UseRandom.h"
0019 #include <ThePEG/Vectors/Lorentz5Vector.h>
0020
0021 namespace Herwig {
0022
0023   using namespace ThePEG;
0024
0025   /** \ingroup Utilities
0026    *  This is a namespace which provides some useful methods
0027    *  for kinematics computation, as the two body decays.
0028    *
0029    *  NB) For other useful kinematical methods (and probably even those
0030    *      implemented in Kinematics class!):
0031    *          @see UtilityBase
0032    */
0033   namespace Kinematics {
0034
0035     /**
0036      *  Calculate the momenta for a two body decay
0037      * The return value indicates success or failure.
0038      * @param p The momentum of the decaying particle
0039      * @param m1 The mass of the first decay product
0040      * @param m2 The mass of the second decay product
0041      * @param unitDir1 Direction for the products in the rest frame of
0042      * the decaying particle
0043      * @param p1 The momentum of the first decay product
0044      * @param p2 The momentum of the second decay product
0045      */
0046     bool twoBodyDecay(const Lorentz5Momentum & p,
0047                  const Energy m1, const Energy m2,
0048                  const Axis & unitDir1,
0049                  Lorentz5Momentum & p1, Lorentz5Momentum & p2);
0050
0051     /**
0052      * It returns the unit 3-vector with the given  cosTheta  and  phi.
0053      */
0054     inline Axis unitDirection(const double cosTheta, const double phi) {
0055       return ( fabs( cosTheta ) <= 1.0  ?
0056          Axis( cos(phi)*sqrt(1.0-cosTheta*cosTheta) ,
0057          sin(phi)*sqrt(1.0-cosTheta*cosTheta) , cosTheta) : Axis() );
0058     }
0059
0060     /**
0061      *  Calculate the momenta for a two body decay
0062      * The return value indicates success or failure.
0063      * @param p The momentum of the decaying particle
0064      * @param m1 The mass of the first decay product
0065      * @param m2 The mass of the second decay product
0066      * @param cosThetaStar1 Polar angle in rest frame
0067      * @param phiStar1 Azimuthal angle in rest frame
0068      * @param p1 The momentum of the first decay product
0069      * @param p2 The momentum of the second decay product
0070      */
0071     inline bool twoBodyDecay(const Lorentz5Momentum & p,
0072                  const Energy m1, const Energy m2,
0073                  const double cosThetaStar1,
0074                  const double phiStar1,
0075                  Lorentz5Momentum & p1, Lorentz5Momentum & p2) {
0076       return twoBodyDecay(p,m1,m2,unitDirection(cosThetaStar1,phiStar1),p1,p2);
0077     }
0078
0079     /**
0080      * As the name implies, this takes the momentum p0 and does a flat three
0081      * body decay into p1..p3. The argument fcn is used to add additional
0082      * weights. If it is not used, the default is just flat in phasespace.
0083      * The return value indicates success or failure.
0084      */
0085     bool threeBodyDecay(Lorentz5Momentum p0, Lorentz5Momentum &p1,
0086                    Lorentz5Momentum &p2, Lorentz5Momentum &p3,
0087                    double (*fcn)(Energy2,Energy2,Energy2,InvEnergy4) = NULL);
0088
0089     /**
0090      * For the two body decay  M -> m1 + m2  it gives the module of the
0091      * 3-momentum of the decay product in the rest frame of M.
0092      */
0093     inline Energy pstarTwoBodyDecay(const Energy M,
0094                     const Energy m1, const Energy m2) {
0095       return ( M > ZERO &&  m1 >=ZERO && m2 >= ZERO  && M > m1+m2 ?
0096            Energy(sqrt(( sqr(M) - sqr(m1+m2) )*( sqr(M) - sqr(m1-m2) ))
0097               / (2.0*M) ) : ZERO);
0098     }
0099
0100     /**
0101      * This just generates angles. First flat -1..1, second flat 0..2Pi
0102      */
0103     inline void generateAngles(double & ct, double & az) {
0104       ct = UseRandom::rnd()*2.0 - 1.0;  // Flat from -1..1
0105       az = UseRandom::rnd()*Constants::twopi;
0106     }
0107   }
0108
0109 }
0110
0111 #endif /* HERWIG_Kinematics_H */
0112