---

PositronRange.cpp

#include <math.h>
#include <iostream.h>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
//
// Compute range of positron from dE/dx (Berger-Seltzer formula) 
// and annihilation.
// This function is part of dE/dx.
//
static double Fp(double t,double tup){
  if(t<=0||tup<=0||t<=tup)return 0;
  double y = 1/(2+t);
  return log(t*tup)
    -tup*tup*(t+2*tup-3*tup*tup*y/2
              -(tup-exp(3*log(tup))/3)*y*y
              -(tup*tup/2
                -t*exp(3*log(tup))/2
                +exp(4*log(tup))/4)*exp(3*log(y)));
}
//
// This function is the positron annihilation cross section minus an overall
// scale.
//
static double AnniFun(double g){
  return ((g*g+4*g+1)/(g*g-1)*log(g+sqrt(g*g-1))
          -(g+3)/sqrt(g*g-1))/(1+g);
}

void PositronRange(double KE,double& range,double& time){

  static ReactorConstants k;

  double step=0.01;
  double Emin=0.01;

  range = 0;
  time = 0;
  if(KE<=Emin)return;

  static double re = k.Relectron;
  static double pi = k.pi;
  static double me = k.Melectron;
  static double Na = k.Navogadro;
  static double c = k.c;
  static double density = k.density;
  //
  // For electron density.
  //
  static double nel[2];
  static double A[2] = {k.A0,k.A1};
  static double Z[2] = {k.Z0,k.Z1};
  static double f[2] = {k.f0,k.f1};
  static double I[2] = {k.I0,k.I1};
  static bool first = 1;
  for(int i=0;i<2;i++){
    nel[i]=Na*density*f[i]*Z[i]/A[i];
    first = 0;
  }
  //
  // Begin energy loss loop
  //
  double dE=0;
  while(KE-dE>=Emin){
    range+=step;
    double gamma = 1+(KE-dE)/me;
    double beta = sqrt(1-1/gamma/gamma);
    time+=step/beta/c;
    //
    // check for annihilation during step
    //
    double AnniMfp = 1/((nel[0]+nel[1])*pi*re*re*AnniFun(gamma));
    int len=1;
    float rv[len];
    ranlux_(rv,len);
    if(rv[0]<1-exp(-step/AnniMfp))return;
    //
    // increment energy loss
    //
    double t = gamma-1;
    double tc = Emin/me;
    double tmax = t;
    double tup = tc;
    if(tc>tmax)tup=tmax;
    for (int k=0;k<2;k++){
      dE+=2*pi*re*re*me*nel[k]*step/beta/beta*
        (log(2*(gamma+1)/(I[k]*I[k]/me/me))
         +Fp(t,tup));
    }
  }
  return;
}
         
                                     
        

---