---

brems.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: brems.cxx,v 1.1 2000/09/01 20:20:18 burnett Exp $
//
// This file is part of Gismo 2
//
// created: 1-Dec-1991 T. Burnett

// Contents -------------------------------------------------------------
//
//   void EGS::brems(PEGSData&)
//
// End ------------------------------------------------------------------

#include "EGS.h"
#include "PEGSData.h"
#include "Random.h"

inline int  ilog2(float x)
{
  return (int)(1.44269*log(x));  //  "1.44269=1/LN(2)"
}

static float ai2ln2=0.7213475; //1./(2.*alog(2.))

#define MAXPWR2I 50
float pwr2i[MAXPWR2I]={0};

// used below in calculation of brem angle
static double set_brem_rejection_function(double input, double RJARG1, 
                                          double RJARG2, double RJARG3, double ZTARG, double ESEDEI);


void EGS::brems(const PEGSData & mat)
{
   if (pwr2i[0]==0){float p=2.; for(int i=0;i<MAXPWR2I;i++)pwr2i[i]=(p/=2.);}

   float eie = this->e;  // incoming electron/positron energy
   //******************************************************************
   //   for electron energy greater than 5.0 mev, the bethe-heitler
   //   cross section is employed.
   //******************************************************************

   //   decide which distribution to use (b-h coulomb corrected is
   //   used from 50 to 20000 mev, b-h is used 1.5 to 50 mev)

   int lvx=1, lvl0=3, lvl;
   if (eie < 50.0){lvx=0;lvl0=0;}

   //   the method of butcher and messel for sampling a class of
   //   factorizable frequency distributions is used.
   //   our 'br' variable is the same as their 'epsilon' variable.
   //   see butcher and messel,nucl.phys.,vol.20,pp23,24.
   //   compute number of subdistributions needed to produce photons
   //   of minimum discrete transport enrgy ap, in case the (1-br)/br
   //   part of the distribution is used.

   float abrems=float(ilog2( eie/mat.ap));

   //   ilog2(x) is that integer function of x such that . . .
   //   2**(ilog2(x)-1) .le. x .lt. 2**(ilog2(x))

   float esg, ese, br;

   for(;;)
   {
       //various rejections can cause resample
       //   decide whether to sample from the (1-br)/br or the 2*br part
       //   of the distribution.

       float rnno06 = Random::flat();
       if (0.5 < ((abrems*mat.alphi(lvx)+0.5)*rnno06))
       {
           //   use the (1-br)/br part.  which subdistribution?
           float rnno07 = Random::flat();
           int idistr =(int)(abrems*rnno07);

           //   this chooses idistr at random from the set
           //   (0,1,2, . . . , nbrems - 1 )
           float p= pwr2i[idistr];
           //   select screening rejection function
           //   lvl=0    uncoulomb corrected     a(delta)
           //   lvl=1    uncoulomb corrected     b(delta)
           //   lvl=2    uncoulomb corrected     c(delta)
           //   lvl=3      coulomb corrected     a(delta)
           //   lvl=4      coulomb corrected     b(delta)
           //   lvl=5      coulomb corrected     c(delta)
           lvl=lvl0;

           //   use a(delta), either born or coulomb corrected, depending on
           //   whether lvl has been previously set to 0 or 3.
           //   all subdistributions are sampled by first sampling from
           //            (1./log(2.))*(1.-br)/br     if 0.5 .le. br .le. 1.
           //            1./log(2.)                  if   br.lt. 0.5
           //   and then taking br = br*p
           //   ai2ln2 is actually 1./(2.*log(2.)), which is the probability
           //   that br is less than 0.5 in the elementary distribution above.
           float rnno08 = Random::flat();

           if (rnno08>=ai2ln2)
           {
               for(;;)
               {  float rnno09 = Random::flat(), //if rejected for subdistribution
                        rnno10 = Random::flat(),
                        rnno11 = Random::flat(),
                        h=std::max(rnno10,rnno11);
                  br=1.0-0.5*h;

                  if ( rnno09 <= 0.5/br) break;   //rejection condition
               }

           }   //end br.ge.0.5 part
           else
           {
               //sample br.lt.0.5 part
               float rnno12 = Random::flat();
               br=rnno12*0.5;
           } //end of br.lt.0.5 part

           //   product energy fraction chosen
           br=br*p;
       } //end (1-br)/br part
       else
       {
          //use the 2*br part
           float rnno13 = Random::flat();
           float rnno14 = Random::flat();
           br=std::max(rnno13,rnno14);
           lvl=lvl0+1; //use b(delta) for screening function
       } //end of 2*br part of distribution

       //   energy of new photon

       esg=eie*br; //energy of secondary gamma

       //   ap=0.256 mev  ---  rm=0.511 mev
       //   ap is selected in the routine pegs.
       //   minimum hardness requirement, corresponding to lower bound
       //   choice for total cross section integral

       if (esg < mat.ap )continue;

       ese=(double)eie-(double)esg; //precise energy of secondary 'electron'

       //   the electron must have a minimum energy equal to 0.511 mev

       if (ese < rm ) continue;

       //   delta=136.0*exp(zg)*rm*ee/(e*(1.0-ee))
       //        =delcm*ee/(e*(1.0-ee))=delcm*del
       //   where e=electron incident energy(MeV), and ee=(photon energy)/e
       //   zg is defined in the program shinp, and is a weighted average
       //   of log(z**(-1./3.))  over the various types of atoms in the
       //   molecule (Butcher and Messel, op.cit., p.17-19,22-24).

       float del = br/ese;

       float rejf = mat.screening(del, lvx, lvl);
       if (rejf == 0) continue;

       float rnscrn = Random::flat(); //random number for screening rejection
       if (rnscrn <= rejf) break;
    }  //loop until value accepted

    this->e = ese;    // change electron or positron energy, not angle

    //now, get photon's direction 

    double theta;
    
    //set IBRDST to 0 for EGS4 default rm/eie;
    //set IBRDST to 1 for Koch & Motz angle distribution
    //Adapted from Alex F. Bielajew's paper
    //"Improved bremsstrahlung photon angular sampleing
    //in the EGS4 code system"  Nov, 14, 1989
    
    int IBRDST = 1; //should change to #ifdef KOCHnMOTZ 
    if (IBRDST == 0){
        theta = rm / eie;
    }
    if (IBRDST == 1){
        double rm2,test,HPI,ZTARG,TTEIE,TTESE,ESEDEI,Y2MAX;
        double RJARG1,RJARG2,RJARG3,REJMIN,REJMID,REJMAX,REJTOP;
        double Y2TST,REJTST,RTEST;
                
        rm2 = 0.51100340; //ELECTRON REST MASS MEV/CC
        HPI = 3.141593;
        ZTARG = (double) mat.hzbrang(); //((1/111)*zeff**(1/3))**2
        TTEIE = (double) eie / rm2; //incoming e- energy in rest mass units
        TTESE = (double) ese / rm2; //outgoing e- energy in rest mass units
        ESEDEI = TTESE/TTEIE;
        Y2MAX = pow((HPI*TTEIE),2);
        RJARG1 = (1.0+pow(ESEDEI,2));
        RJARG2 = 3.0 * RJARG1 - 2.0*ESEDEI;
        RJARG3 = pow(((1.0-ESEDEI)/(2.0*TTEIE*ESEDEI)),2);
        REJMIN = set_brem_rejection_function(0.0E0,RJARG1,RJARG2,RJARG3,ZTARG,ESEDEI);
        REJMID = set_brem_rejection_function(1.0E0,RJARG1,RJARG2,RJARG3,ZTARG,ESEDEI);
        REJMAX = set_brem_rejection_function(Y2MAX,RJARG1,RJARG2,RJARG3,ZTARG,ESEDEI);
        test = REJMIN;
        if (test <= REJMID) test = REJMID;
        if (test <= REJMAX) test = REJMAX;
        REJTOP = test;
        do {
            Y2TST = (double) Random::flat();
            Y2TST = Y2TST/(1.0-Y2TST+1.0/Y2MAX);
            REJTST = set_brem_rejection_function(Y2TST,RJARG1,RJARG2,RJARG3,ZTARG,ESEDEI);
            RTEST = (double) Random::flat();
        } while(RTEST > (REJTST/REJTOP));
        if (Y2TST < 0){ //error protection
            Y2TST = 0;
        }
        
        theta = sqrt(Y2TST)/TTEIE;
        //      if (theta > 3.1415926) theta = 3.141592; //to be safe
        //      if (theta < 0) theta = 0; //to be safe
    }
    
    EGS *photon = new EGS(this);  // create the photon
    photon->iq = 0;
    photon->e = esg;
    photon->rotate( cos(theta), Random::flat(2.*M_PI) );

} // end of PDT::brems

static double set_brem_rejection_function(double input, double RJARG1, 
                                          double RJARG2, double RJARG3, double ZTARG, double ESEDEI)
{ //harrison
    double Y2TST1, result;
    Y2TST1 = pow((1.0 + input),2);
    result = (4.0 + log(RJARG3+ZTARG/Y2TST1))*(4.0*ESEDEI*input/Y2TST1-RJARG1)+RJARG2;
    return result;
}

---