---

GammaComptonRange.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream.h>
#include <string>
#include "ReactorConstants.hh"

double GammaComptonRange(std::string material,double egamma){

  //cout << "gamma propogating through " << material << endl;
  //
  // Compute the mean free path for Compton scattering.
  //
  static ReactorConstants k;
  static double pi = k.pi;
  static double m = k.Melectron;
  static double alpha = 1/k.alphainv;
  static double Units = k.XcMeVtoCmsqrd;
  //
  // Overall cross section scale
  //
  static double Thomson = 8*pi*alpha*alpha/3/m/m*Units;
  double w = egamma/m;
  //
  // Klein-Nishina cross section
  //
  double XC = 3*Thomson/4*
    ((1+w)/w/w/w*(2*w*(1+w)/(1+2*w)-log(1+2*w))
     +log(1+2*w)/2/w-(1+3*w)/(1+2*w)/(1+2*w));
  //
  // Need electron density to get mfp. Get H and C contributions.
  //
  static double density = k.density;
  static double Na = k.Navogadro;
  static double A[2] = {k.A0,k.A1};
  static double Z[2] = {k.Z0,k.Z1};
  static double f[2];
  if(material == "scintillator"){
    f[0] = k.f0;
    f[1] = k.f1;
  }
  else if(material == "oil"){
    f[0] = k.f4;
    f[1] = k.f5;
  }
  else{
    cout << "GammaComptonRange: " << material 
         << " is an unknown material!  Assuming scintillator." << endl;
    f[0] = k.f0;
    f[1] = k.f1;
  }    
  //for(int i=0;i<2;i++) cout << "f["<<i<<"] " << f[i] << endl;
  static double nel[2];
  static bool first = 1;
  for(int i=0;i<2;i++){
    nel[i]=Na*density*f[i]*Z[i]/A[i];
    first = 0;
  }
  //
  //  Compute mean free path.
  //  
  return 1/((nel[0]+nel[1])*XC);
}
                          

---