Kinematics.h

// -*- C++ -*-
//
// Kinematics.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_Kinematics_H
#define HERWIG_Kinematics_H
 
// This is the declaration of the Kinematics class.
 
#include <ThePEG/Config/ThePEG.h>
#include "ThePEG/Vectors/ThreeVector.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/Repository/UseRandom.h"
#include <ThePEG/Vectors/Lorentz5Vector.h>
 
namespace Herwig {
 
  using namespace ThePEG;
 
  namespace Kinematics {
 
    bool twoBodyDecay(const Lorentz5Momentum & p, 
                             const Energy m1, const Energy m2,
                             const Axis & unitDir1,
                             Lorentz5Momentum & p1, Lorentz5Momentum & p2);
    
    inline Axis unitDirection(const double cosTheta, const double phi) {
      return ( fabs( cosTheta ) <= 1.0  ? 
         Axis( cos(phi)*sqrt(1.0-cosTheta*cosTheta) , 
         sin(phi)*sqrt(1.0-cosTheta*cosTheta) , cosTheta) : Axis() );
    }
 
    inline bool twoBodyDecay(const Lorentz5Momentum & p, 
                             const Energy m1, const Energy m2,
                             const double cosThetaStar1, 
                             const double phiStar1,
                             Lorentz5Momentum & p1, Lorentz5Momentum & p2) {
      return twoBodyDecay(p,m1,m2,unitDirection(cosThetaStar1,phiStar1),p1,p2); 
    }
 
    bool threeBodyDecay(Lorentz5Momentum p0, Lorentz5Momentum &p1, 
                               Lorentz5Momentum &p2, Lorentz5Momentum &p3,
                               double (*fcn)(Energy2,Energy2,Energy2,InvEnergy4) = NULL);
 
    inline Energy pstarTwoBodyDecay(const Energy M, 
                                    const Energy m1, const Energy m2) {
      return ( M > ZERO &&  m1 >=ZERO && m2 >= ZERO  && M > m1+m2 ?  
               Energy(sqrt(( sqr(M) - sqr(m1+m2) )*( sqr(M) - sqr(m1-m2) )) 
                      / (2.0*M) ) : ZERO); 
    }
    
    inline void generateAngles(double & ct, double & az) {
      ct = UseRandom::rnd()*2.0 - 1.0;  // Flat from -1..1
      az = UseRandom::rnd()*Constants::twopi;   
    }
  }
 
}
 
#endif /* HERWIG_Kinematics_H */