---

Ranbo.cxx

Go to the documentation of this file.

//
// This file is part of Gismo 2
// Implementation of the phase space generator Ranbo

#include "Ranbo.h"

#include "Random.h"  //TODO fix
static inline double sqr(double x) { return x*x;}

HepRandom* ran;  // file scope to pass to Ran4Vec constructor

Ranbo::Ranbo(double W, unsigned n, double mu[], HepRandom& _ran)
{

        // create a set of randomized 4-momenta (see constructor)
        ran = &_ran;
        q =  new Ran4Vec[n] ;

        HepLorentzVector qsum;
        unsigned i;
        for(  i=0; i<n; i++)
          qsum += q[i];

        double x = W/qsum.mag();
        Hep3Vector beta = -qsum.boostVector();

        // now boost to CM and rescale to correct total energy
        for(i=0; i<n; i++)
        {
                q[i].boost(beta);
                q[i] *= x;
        }

        // iterate to find factor to multiply all 3-momenta so that
        // total energy is OK
        double xisq = 1;  // start from here
        double* musq = new double[n];
        double* psq = new double[n];
        double* energy = new double[n];
        for( i=0; i<n; i++)
        {
                musq[i]= sqr(mu[i]);
                psq[i] = sqr(q[i].t());
        }
        if( n==2 ) // this case is exact
        {
                xisq = (1-sqr((mu[0]+mu[1])/W))
                                *(1-sqr((mu[0]-mu[1])/W));
        }
        else   // otherwise iterate
        {
          double f,df;   // E-conservation function & deriviate
          do
          {
                 f = -W; df=0;
                 for( i=0; i<n; i++)
                 {  double e = energy[i] = sqrt(musq[i]+xisq*psq[i]);
                        f += e;
                        df += psq[i]/e;
                 }
                 df *= 0.5;
                 xisq -= f/df;
          }while( fabs(f)>1e-6*W ); // this should always converge
        }
        if( xisq != 1 )
        {
          double xi = sqrt(xisq);
          for( i=0; i<n; i++)
                 (Hep3Vector&)q[i] *= xi;
        }
        // finally calculate weight
        double numerator=W,denominator=0;
        for(i=0; i<n; i++)
        {
                double betasq = xisq*psq[i]/sqr(energy[i]);
                numerator *= sqrt(betasq);
                denominator += betasq*energy[i];
        }
        _weight = numerator/denominator*pow(sqrt(xisq),int(2*n-3));

        delete [] musq;
        delete [] psq;
        delete [] energy;
}

Ranbo::~Ranbo()
{
        delete []q;
}


Ranbo::Ran4Vec::Ran4Vec()
{
  double c = -1 + 2*ran->flat();
  double phi = 2*M_PI*ran->flat();
  double r = -log(ran->flat()*ran->flat());
  double s = sqrt(1-c*c);
  setX(r*s*cos(phi));
  setY(r*s*sin(phi));
  setZ(r*c);
  setT(r);
}

Ranbo::Ran4Vec &
Ranbo::Ran4Vec::operator *= (float f)
{
          setX(f*x());
          setY(f*y());
          setZ(f*z());
          setT(f*t());
          return *this;
}

---