Main Page   Compound List   File List   Compound Members   File Members  

PositronRange.cpp

00001 #include <math.h>
00002 #include <iostream.h>
00003 #include "ReactorConstants.hh"
00004 #include "ReactorFortran.hh"
00005 //
00006 // Compute range of positron from dE/dx (Berger-Seltzer formula) 
00007 // and annihilation.
00008 // This function is part of dE/dx.
00009 //
00010 static double Fp(double t,double tup){
00011   if(t<=0||tup<=0||t<=tup)return 0;
00012   double y = 1/(2+t);
00013   return log(t*tup)
00014     -tup*tup*(t+2*tup-3*tup*tup*y/2
00015               -(tup-exp(3*log(tup))/3)*y*y
00016               -(tup*tup/2
00017                 -t*exp(3*log(tup))/2
00018                 +exp(4*log(tup))/4)*exp(3*log(y)));
00019 }
00020 //
00021 // This function is the positron annihilation cross section minus an overall
00022 // scale.
00023 //
00024 static double AnniFun(double g){
00025   return ((g*g+4*g+1)/(g*g-1)*log(g+sqrt(g*g-1))
00026           -(g+3)/sqrt(g*g-1))/(1+g);
00027 }
00028 
00029 void PositronRange(double KE,double& range,double& time){
00030 
00031   double step=0.01;
00032   double EminPos=0.01;
00033 
00034   range = 0;
00035   time = 0;
00036   if(KE<=EminPos)return;
00037 
00038   static double re = RELECTRON;
00039   static double me = MELECTRON;
00040   static double Na = NAVOGADRO;
00041   //
00042   // For electron density.
00043   //
00044   static double nel[2];
00045   static double A[2] = {A0,A1};
00046   static double Z[2] = {Z0,Z1};
00047   static double f[2] = {f0,f1};
00048   static double I[2] = {I0,I1};
00049   static bool first = 1;
00050   for(int i=0;i<2;i++){
00051     nel[i]=Na*density*f[i]*Z[i]/A[i];
00052     first = 0;
00053   }
00054   //
00055   // Begin energy loss loop
00056   //
00057   double dE=0;
00058   while(KE-dE>=EminPos){
00059     range+=step;
00060     double gamma = 1+(KE-dE)/me;
00061     double beta = sqrt(1-1/gamma/gamma);
00062     time+=step/beta/c;
00063     //
00064     // check for annihilation during step
00065     //
00066     double AnniMfp = 1/((nel[0]+nel[1])*PI*re*re*AnniFun(gamma));
00067     const int len=1;
00068     float rv[len];
00069     ranlux_(rv,len);
00070     if(rv[0]<1-exp(-step/AnniMfp))return;
00071     //
00072     // increment energy loss
00073     //
00074     double t = gamma-1;
00075     double tc = EminPos/me;
00076     double tmax = t;
00077     double tup = tc;
00078     if(tc>tmax)tup=tmax;
00079     for (int k=0;k<2;k++){
00080       dE+=2*PI*re*re*me*nel[k]*step/beta/beta*
00081         (log(2*(gamma+1)/(I[k]*I[k]/me/me))
00082          +Fp(t,tup));
00083     }
00084   }
00085   return;
00086 }
00087          
00088                                      
00089         

Generated on Mon Feb 21 16:11:19 2005 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002