GammaCompton.cpp Source File

00001 #include "ReactorConstants.hh"
00002 #include "ReactorFortran.hh"
00003 #include "ReactorXC.hh"
00004 #include <cmath>
00005 #include <iostream>
00006 using namespace std;
00007 
00008 void GammaCompton(ReactorConstants& RC,double* Pg){
00009   //
00010   // RC = reactor constants
00011   // Pg = in and out gamma 4 vector: e,px,py,pz
00012   //
00013   static double pi = RC.pi;
00014   static double m = RC.Melectron;
00015   //
00016   // sample Klein-Nishina cross section
00017   //
00018   double ki = Pg[0];
00019   double w = ki/m;
00020   double kf;
00021   double cosg;
00022   double phig;
00023   double wtmax=2;
00024   bool OK=0;
00025   while(!OK){
00026     int len=3;
00027     float rv[len];
00028     ranlux_(rv,len);
00029     //
00030     cosg = -1+2*rv[0];
00031     phig = 2*pi*rv[1];
00032     //
00033     kf = ki/(1+w*(1-cosg));
00034     double KNXC = kf*kf/ki/ki*(kf/ki+ki/kf-1+cosg*cosg);
00035     //
00036     double test = wtmax*rv[2];
00037     OK=test<KNXC;
00038   }
00039   //
00040   double sing=sqrt(1-cosg*cosg);
00041   //
00042   double t[3];
00043   t[0] = Pg[1]/ki;
00044   t[1] = Pg[2]/ki;
00045   t[2] = Pg[3]/ki;
00046   double v[3];
00047   double pt = sqrt(Pg[1]*Pg[1]+Pg[2]*Pg[2]);
00048   v[0]=-Pg[2]/pt;
00049   v[1]= Pg[1]/pt;
00050   v[2]=0;
00051   double u[3];
00052   u[0] = v[1]*t[2]-v[2]*t[1];
00053   u[1] = v[2]*t[0]-v[0]*t[2];
00054   u[2] = v[0]*t[1]-v[1]*t[0];
00055   //
00056   for(int i=0;i<3;i++){
00057     Pg[i+1] = kf*t[i]*cosg;
00058     Pg[i+1]+= kf*u[i]*sing*cos(phig);
00059     Pg[i+1]+= kf*v[i]*sing*sin(phig);
00060   }
00061   Pg[0] = kf;
00062   //  cout << " ki, kf, cosg = " << ki << " " << kf << " " << cosg << endl;
00063   return;
00064 }
00065