bhabha.cxx Source File

00001 // -*- C++ -*-  $Id: bhabha.cxx,v 1.1 2000/09/01 20:20:18 burnett Exp $
00002 //
00003 // This file is part of Gismo 2
00004 //
00005 // file: bhabah.cpp
00006 // created 1-Dec-1991  T. Burnett
00007 
00008 // Contents -------------------------------------------------------
00009 //
00010 //  EGS::bhabha
00011 //
00012 // End -----------------------------------------------------------
00013 
00014 #include "EGS.h"
00015 #include "PEGSData.h"
00016 #include "Random.h"
00017 
00018 
00019 void EGS::bhabha(const PEGSData& mat)
00020 {
00021     //                                version 4.00  --  26 jan 1986/1900
00022     //******************************************************************
00023     //   discrete bhabha scattering (a call to this routine) has been
00024     //   arbitrarily defined and calculated to mean bhabha scatterings
00025     //   which impart to the secondary electron sufficient energy that
00026     //   it be transported discretely, i.e. e=ae or t=te.  it is not
00027     //   guaranteed that the final positron will have this much energy
00028     //   however.  the exact bhabha differential cross section is used.
00029     //******************************************************************
00030 
00031     double peip = this->e,         //energy of incident positron
00032            prm  = rm,              // electron mass
00033            pekin=peip-prm;         //precise k.e. of incident positron
00034     float  ekin =pekin,
00035            br,
00036 
00037           t0 = ekin/rm,
00038           e0 = t0+1.,
00039           yy = 1./(t0+2.),
00040           e02= e0*e0,
00041           betai2=e02/(e02-1.),
00042           ep0 = mat.te/ekin,   // te is a data member
00043           ep0c=1.-ep0,
00044           y2=yy*yy,
00045           yp=1.-2.*yy;
00046     float yp2=yp*yp,
00047           b4 =yp2*yp,
00048           b3 =b4+yp2,
00049           b2 =yp*(3.+y2),
00050           b1 =2.-y2;
00051     //   sample br from minimum(ep0) to 1.
00052 
00053     for(;;)
00054     {   float rnno03 = Random::flat();
00055         br=ep0/(1.-ep0c*rnno03);
00056         //   apply rejection function
00057         float rnno04 = Random::flat();
00058         float rejf2=ep0c*(betai2-br*(b1-br*(b2-br*(b3-br*b4))));
00059         if ( rnno04<=rejf2) break ;
00060     }
00061 
00062 
00063 
00064     double pekse2=br*ekin; //precise kinetic energy of secondary 'electron' 2
00065     double pese1=peip-pekse2; //precise energy of secondary 'electron' 1
00066     double pese2=pekse2+prm;  //precise energy of secondary 'electron' 2
00067 
00068     //   bhabha angles are uniquely determined by kinematics
00069 
00070     double h1=(peip+prm)/pekin;
00071 
00072     //   direction cosine change for 'old' electron
00073     float phi = Random::flat(2.*M_PI);   // random phi
00074     double dcosth=h1*(pese1-prm)/(pese1+prm);
00075     float costhe = sqrt(dcosth);
00076     rotate(costhe, phi);
00077     this->e = pese1;
00078 
00079     EGS& el2 = *new EGS(this);   // create new electron
00080     el2.iq=-1;
00081     el2.e = pese2;
00082     dcosth=h1*(pese2-prm)/(pese2+prm);
00083     costhe= sqrt(dcosth);
00084     el2.rotate(costhe, phi+M_PI);
00085 
00086 }  //end of subroutine bhabha
00087