---

PhaseSpace.cxx

Go to the documentation of this file.

// -*- C++ -*-  
// $Id: PhaseSpace.cxx,v 1.4 2001/01/25 22:24:32 burnett Exp $
//
// This file is part of Gismo 2
//

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "gismo/PhaseSpace.h"

#include "gismo/GParticle.h"
#include <cmath>
#include "Random.h"
#include "facilities/error.h"
#include <algorithm>
#include <cmath>
static inline double sqr(double x) { return x*x;}

int PhaseSpace::decay(GParticle* parent)const
{
        weight = 1;

        int n =parent->numChildren();
        if( n==2 ) return decay2( parent);
        else if(n==3) return decay3( parent );
        else
                FATAL("PhaseSpace: not set up to decay other than 2 or 3 secondaries");
        return 0;
}
int PhaseSpace::decay2(GParticle* parent)const
{

        GParticle& b = *parent->child(0);
        GParticle& c = *parent->child(1);
        float m1 = b.mass();
        float m2 = c.mass();

        float M = parent->mass();

        float pstar = (sqr(M)-sqr(m1+m2))*(sqr(M)-sqr(m1-m2));
        if (pstar<0.)
        {  WARNING( "*** decay kinematically forbidden");
                return 0;
        }
        pstar = sqrt(pstar)/(2.0*M);
        float costh = Random::flat(-1.0, 1.0);
        float phi   = Random::flat(2*M_PI);
        Vector pcm(0,0,pstar);
        pcm.Hep3Vector::rotateY(acos(costh));  // rotate about y axis by acos(costh)
        pcm.Hep3Vector::rotateZ(phi);         // rotate about z axis by phi

        b.setMomentum(pcm);
        c.setMomentum(-pcm);

        return 1;
}
int PhaseSpace::decay3(GParticle* parent)const
{
        GParticle& a = *parent->child(0);
        GParticle& b = *parent->child(1);
        GParticle& c = *parent->child(2);
        double m1 = a.mass();
        double m2 = b.mass();
        double m3 = c.mass();

        double M = parent->mass();


         /* Function Body */
         double t0 = M - m1 - m2 - m3;
         if( t0 <= 0 ) return 0;   // not kinematically allowed
         double te, tm, tn;
         double pa1, pa2, pa3, pmax;
         do
         {
                double g1 = Random::flat();
                double g2 = Random::flat();
                double gmi = std::min(g1,g2);
                double gma = std::max(g1,g2);
                te = gmi * t0;
                tm = (1. - gma) * t0;
                tn = (gma - gmi) * t0;
                pa1 = sqrt(sqr(te) + m1 * 2. * te);
                pa2 = sqrt(sqr(tm) + m2 * 2. * tm);
                pa3 = sqrt(sqr(tn) + m3 * 2. * tn);
/* -------------------------------------------------------- */
/*     CHECK OF THE MOMENTUM CONSERVATION */
/* -------------------------------------------------------- */
                pmax = std::max(pa3,std::max(pa1,pa2));
         }
         while( pmax >= pa1 + pa2 + pa3 - pmax );

/* -------------------------------------------------------- */
/*     CALCULATIONS OF MOMENTA */
/* -------------------------------------------------------- */
         double cte = Random::flat(-1,1);
         double ste = sqrt( abs(1. - sqr(cte)) );
         double fe = Random::flat(0, 2*M_PI);
         double cfe = cos(fe);
         double sfe = sin(fe);
         Hep3Vector pa(pa1 * ste * cfe,pa1 * ste * sfe,pa1 * cte);
         a.setMomentum(pa);
         double cten = (sqr(pa2) - sqr(pa1) - sqr(pa3)) / (pa1 * 2. * pa3);
         double sten = sqrt(abs(1. - sqr(cten)) );
         double fen = Random::flat(0,2*M_PI);
         double cfen = cos(fen);
         double sfen = sin(fen);
         Hep3Vector pc(
                pa3 * ( sten * cfen * cte * cfe - sten * sfen * sfe + cten * ste * cfe)
          ,pa3 * ( sten * cfen * cte * sfe + sten * sfen * cfe + cten * ste * sfe)
          ,pa3 * (-sten * cfen * ste + cten * cte)
          );
         c.setMomentum( pc);
         b.setMomentum( -pa-pc );
         return 1;
}

---