---

EGS.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: EGS.cxx,v 1.2 2001/01/25 22:24:33 burnett Exp $
//
// This file is part of Gismo 2
//
// file: egs.cpp
// created 1-Dec-1991, T. Burnett


// Contents ----------------------------------------------------------------
//
//      EGS::EGS
//      EGS::~EGS
//      EGS::rotate
//      EGS::interact
//      EGS::meanFreePath
//      EGS::dedx
//
// Description
//
//      Implementation of EGS particle routines
//
// End --------------------------------------------------------------------

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

#include "EGS.h"

#include "PEGSData.h"

#include "Random.h"
#include "facilities/error.h"
#include <math.h>
#include <iomanip>

// ----------------constructor EGS::EGS ----------------------------------

EGS::EGS(EGS* old)
{
    // if arg present, assume that it is the top of the particle stack:
    // add this entry to it, and copy all data

    if (old )
    {
        *this = *old;   // copy data
        old->next = this; // link in
        prev = old;
        next = 0;

    }
    else
    {   // this is the TOS: set up defaults
        next=prev=0;
        x=y=z=0;
        u=v=0; w=1.;
        iq = 0;
    }
}

// ----------------destructor   EGS::~EGS ------------------------------

EGS::~EGS()
{
   // remove this particle and all following it in the stack

   EGS* pt = next;
   if (prev ) prev->next = 0;

   if( next )
        delete pt;
 }

// ----------------member function EGS::printOn ------------------------

void EGS::printOn(std::ostream& cout)
{   //cout << setprecision(3);

#if 1
    std::cout << std::setw(5)  << iq
         << std::setw(10) << e
         << std::setw(10) << u
         << std::setw(10) << v
         << std::setw(10) << w
         << '\n';
#else
    std::cout << iq << e << u << v << w << '\n';
#endif
    if (next != NULL) next->printOn(cout);
}


// ----------------member function EGS::rotate ------------------------

void EGS::rotate(double costhe, double phi)
{
    double sinthsq = 1.-costhe*costhe;
    rotate(costhe, phi, sinthsq<0? 0: sqrt(sinthsq)  );
}
void EGS::rotate(double costhe, double phi, double sinthe)
{
    // rotate the direction cosines by (theta, phi)
    // --------------------------------------------
    //  see h.h. nagel dissertation for coordinate system description.
    //  a rotation is performed to transform direction cosines of the
    //  particle back to the physical frame (from the transport frame)

    static float phiold=0, sinphi=0, cosphi=1;

    if (phi!=phiold)
    {   if (phi==phiold+M_PI )
        {   // this happpens for the 2nd of correlated particles
            sinphi = -sinphi; cosphi = -cosphi;
        }
        else
        {   sinphi = sin(phi); cosphi = cos(phi);
        }
        phiold = phi;
    }

    float  sinps2=u*u+v*v;

    //   if sinps2 is small, no rotation is needed

    if (sinps2 < 1.0e-20)
    {
        //small polar angle case
        u = sinthe*cosphi;
        v = sinthe*sinphi;
        w = w*costhe;
    }
    else
    {
        //large polar angle case
        float a=u, b=v, c=w,
              sinpsi=sqrt(sinps2),
              us=sinthe*cosphi,
              vs=sinthe*sinphi,
              sindel=b/sinpsi,
              cosdel=a/sinpsi;
        u = c*cosdel*us - sindel*vs + a*costhe;
        v = c*sindel*us + cosdel*vs + b*costhe;
        w =  -sinpsi*us             + c*costhe;
    }
}

// ----------------member function EGS::interact ------------------------

void EGS::interact(const PEGSData &m)
{
    // interact the particle in the material. The particle is assumed
    // to be at the top of the stack (but no check is made) It is
    // replaced by its daughter(s)

  //    int qIn = iq;

    PEGSData::channel which = m.chooseChannel(Random::flat(),e,iq);
    switch (which)    {
        case PEGSData::STOP:    e=rm;     break;
        case PEGSData::BREMS:   brems(m); break;  //e- -> e- + gamma
        case PEGSData::MOLLER:  moller(m);break;  //e- -> e- + e-
        case PEGSData::ANNIH:   annih();  break;  //e+ => gamma+gamma
        case PEGSData::BHABHA:  bhabha(m);break;  //e+ -> e+ + e-
        case PEGSData::PAIR:  pairProd(m);break;  //gamma-> e+ + e-
        case PEGSData::COMPT:   compt();  break;  //gamma-> gamma + e-
        case PEGSData::PHOTO:   photo(m); break;  //gamma-> e- (+gamma)
        case PEGSData::RAYLEIGH: rayleigh(m); break; // gamma->gamma
        

    }
//    if(!next && qIn)
//         WARNING("EGS::interact -- next == 0!! ");
}
//--------------------------------------------------------------------
//------------------member function EGS::rayleigh --------------------
//--------------------------------------------------------------------

/* should duplicate the following EGS4 2.0 mortran:
    :sampling-loop: loop [

        $randomset xxx;
        $set interval xxx,rco;
        $evaluate x2 using rsct(xxx);

        q2=x2*rmsq/(20.60744*20.60744);
        costhe=1.-q2/(2.*e(np)*e(np));
        if (abs(costhe).gt.1.0) go to :sampling-loop:;

        csqthe=costhe*costhe;
        rejf=(1.0+csqthe)/2.0;
        $randomset rnnorj;

    ] until rnnorj.le.rejf;

    sinthe=sqrt(1.0-csqthe);
    call uphi(2,1);

*/
void EGS::rayleigh(const PEGSData &m )
{
    float costhe, rejf;
    do {
        float x2 = m.ph.rsct(Random::flat());
        float q2 = x2*rmsq/(20.60744*20.60744);
        costhe=1.-q2/(2.*e*e);
        if( abs(costhe) > 1.0 ) continue;
        float csqthe=costhe*costhe;
        rejf=(1.0+csqthe)/2.0;
        if( Random::flat() <= rejf ) break;
    } while(1);

    rotate( costhe, Random::flat(2.*M_PI) );
}

// ----------------- static variables ------------------------------
double EGS::rm2 = 0.51100340; // electron mass in MeV
double EGS::rm = rm2;

double EGS::rmsq = rm2*rm2;

---