---

compt.cxx

Go to the documentation of this file.

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

// Contents -------------------------------------------------------
//
//  EGS::compt
//
// End -----------------------------------------------------------


#include "EGS.h"
#include "Random.h"
#include <algorithm>


void EGS::compt()
{
         //                                version 4.00  --  26 jan 1986/1900
         //******************************************************************
         //   Butcher and Messel's cross section expression is used
         //   (Butcher and Messel, op.cit., p. 17-19,25), but the
         //   1/epsilon part is not sampled in the way that they do.
         //   this routine calls their 'epsilon' variable by the name 'br'.
         //   br=final photon energy /initial photon energy.
         //   br0 = minimum br = 1./(1.+2.*(e(np)/rm))
         //   maximum br is 1.
         //   Butcher and Messel's expression for the differential cross
         //   section is proportional to
         //        (1./br+br)*(1.-br*sinthe**2/(1.+br*br))
         //   we shall sample from the first factor from the interval (br0,1)
         //   and use the second factor as a rejection function.
         //******************************************************************

         double prm = rm;
         float  eig = this->e;
         double peig=eig; //precise energy of incident gamma
         float egp=eig/rm,
          br0i=1.+2.*egp,   //    br0i is the inverse of br0
          alph1=log(br0i),
          alph2=egp*(br0i+1.)/(br0i*br0i),
          sumalp = alph1+alph2,
          sinthe,temp,
          a1mibr,
          esg,
          br;
         EGS& newp = *new EGS(this);  // create the electron

         //   due to rejection, some expressions which only appear once in
         //   the code below, may actually be needed more than once.
         for(;;)
         {   //retry if rejected
        //   which part of  1./br + br to sample from ?
        float rnno15 = Random::flat();
        if (alph1 >= sumalp*rnno15)
        {   //use 1/br part of distribution
            float rnno16 = Random::flat();
            br=exp(alph1*rnno16)/br0i;
        }   //end of 1/br part
        else
        {  //use linear ( br ) part of distribution
            float brp = Random::flat();
            float rnno18 = Random::flat();
            if ( egp >= (egp+1.)*rnno18 )
            {  float rnno19 = Random::flat();
               brp=std::max(brp,rnno19);
            }
            br=((br0i-1.)*brp+1.)/br0i;
        } //end sampling of linear part

        //   br=final photon energy fraction
        esg=br*eig; //energy of secondary gamma
        //   the compton angles for photon and recoil electron are uniquely
        //   determined by the conservation laws
        a1mibr = 1.-br;
        temp=rm*a1mibr/esg;
        //   the amax1 in the ff. is to prevent sinthe.lt.0 from trunc error
        sinthe=std::max(0.0,temp*(2.0-temp));
        //   compare rejection function with random number.
        float rnno20 = Random::flat();
        float rejf3=1.0-br*sinthe/(1.0+br*br);
        if ( rnno20 <= rejf3) break; //loop until accepted
    }
    sinthe=sqrt(sinthe);
    float costhe=1.0-temp;

    //   the recoil electron is added to the shower memory.  the extra
    //   rest mass energy will be discarded when the
    //   electron is thrown away.
    //   (how? --THB)

    double pesg=esg;   //precise energy of secondary gamma
    double pese=peig-pesg+prm; //precise energy of secondary electron
    float  ese=pese;   //energy of secondary electron

    float phi = Random::flat(2.*M_PI);       // generate uniform phi
    rotate(costhe, phi);          // and rotate old photon

         //   related change and (x,y,z) setup for new elctron
    float psq=ese*ese-rmsq;

    //   psq here is the momentum squared.
    costhe = (psq<0.)? 0. : (ese+esg)*a1mibr/sqrt(psq);
    newp.rotate(costhe, phi+M_PI);   // and rotate its direction

    // now put the lowest energy particle on top of stack
    if (ese<= esg)
    {  newp.iq = -1;    // electron is on top
       newp.e = pese;
       this->e = pesg;
    }else
    {  newp.iq = 0; this->iq = -1;
       newp.e = pesg;
       this->e = pese;
       float t;
                 t=newp.u; newp.u=u; u=t;
                 t=newp.v; newp.v=v; v=t;
       t=newp.w; newp.w=w; w=t;
    }
}//end of subroutine compt

---