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(double* Pg){
00009   //
00010   // Pg = in and out gamma 4 vector: e,px,py,pz
00011   //
00012   static double m = MELECTRON;
00013   //
00014   // sample Klein-Nishina cross section
00015   //
00016   double ki = Pg[0];
00017   double w = ki/m;
00018   double kf;
00019   double cosg;
00020   double phig;
00021   double wtmax=2;
00022   bool OK=0;
00023   while(!OK){
00024     const int len=3;
00025     float rv[len];
00026     ranlux_(rv,len);
00027     //
00028     cosg = -1+2*rv[0];
00029     phig = 2*PI*rv[1];
00030     //
00031     kf = ki/(1+w*(1-cosg));
00032     double KNXC = kf*kf/ki/ki*(kf/ki+ki/kf-1+cosg*cosg);
00033     //
00034     double test = wtmax*rv[2];
00035     OK=test<KNXC;
00036   }
00037   //
00038   double sing=sqrt(1-cosg*cosg);
00039   //
00040   double t[3];
00041   t[0] = Pg[1]/ki;
00042   t[1] = Pg[2]/ki;
00043   t[2] = Pg[3]/ki;
00044   double v[3];
00045   double pt = sqrt(Pg[1]*Pg[1]+Pg[2]*Pg[2]);
00046   v[0]=-Pg[2]/pt;
00047   v[1]= Pg[1]/pt;
00048   v[2]=0;
00049   double u[3];
00050   u[0] = v[1]*t[2]-v[2]*t[1];
00051   u[1] = v[2]*t[0]-v[0]*t[2];
00052   u[2] = v[0]*t[1]-v[1]*t[0];
00053   //
00054   for(int i=0;i<3;i++){
00055     Pg[i+1] = kf*t[i]*cosg;
00056     Pg[i+1]+= kf*u[i]*sing*cos(phig);
00057     Pg[i+1]+= kf*v[i]*sing*sin(phig);
00058   }
00059   Pg[0] = kf;
00060   //  cout << " ki, kf, cosg = " << ki << " " << kf << " " << cosg << endl;
00061   return;
00062 }
00063