---

SageDecayer.cxx

Go to the documentation of this file.

//  $Header: /nfs/slac/g/glast/ground/cvs/gismo/src/SageDecayer.cxx,v 1.1 2000/01/18 00:33:35 burnett Exp $
//

#include "gismo/SageDecayer.h"
#include "Random.h"

#include "gismo/GParticle.h"
#include <strstream>
#include "facilities/error.h"

typedef long int integer;
typedef float real;

#define GOGEN gogen_
#define SAGEXX sagexx_
#define SAGERR sagerr_
#define RAND  rand_
#define SAGEWT sagewt_


// Constructor //
SageDecayer::SageDecayer(const char *) : Decayer("Sage")
{
}
// This is an interface routine for sage
extern "C"{

void RAND(real r[], integer* n);
void SAGERR(char what[4], integer* code);
extern void GOGEN(real*, real[], integer[], real[], real[], real[][4]);
extern  struct sagexx_1_ {
         real hend, hhlab, hcm, fac1[51];
  } SAGEXX;
extern struct sagewt_1_{
         real swt[2], wsq, w[2];
  } SAGEWT;
}

void RAND(real r[], integer* n)
{
  int i;
  for(i=0;i<*n;i++)
    r[i]=Random::flat();
}
void SAGERR(char what[4], integer* code)
{
   std::ostrstream serr;
   serr << "Sage Error from "<<what << *code <<'\n';
   WARNING(serr.str());
}

// Member function


int SageDecayer::decay(GParticle* parent)const
{

        int i,n = parent->numChildren();
        // uses the sage decay MC to calculate the four vectors of the decay
        //   products (given by 'pl') of particle p

        real U0, P0[4], M[30], U[30], P[30][4];           // 'FORTRAN' vectors
        U0    = parent->mass();
        P0[0] = SAGEXX.hcm;     // set CM flag so Sage knows not to boost
        integer NP[]={n,1};
        for(i=0; i<n; i++)
        {
                M[i] = parent->child(i)->mass();

        }
        SAGEWT.w[0] = -1.;
        GOGEN(&U0, P0, NP, M, U, P);
  //          WSQ=SQRT(WSQ)
  //    A1(1)=ABS(W(1))*WSQ*SWT(1)

        weight = SAGEWT.w[0]*sqrt(SAGEWT.wsq)*SAGEWT.swt[0];;
        for(i=0; i<n; i++)
        {
                parent->child(i)->setMomentum(Vector(P[i][0],P[i][1],P[i][2]));
        }
   return 1;
}

---