G4eeToTwoGammaModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:   G4eeToTwoGammaModel
00034 //
00035 // Author:        Vladimir Ivanchenko on base of Michel Maire code
00036 //
00037 // Creation date: 02.08.2004
00038 //
00039 // Modifications:
00040 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
00041 // 18-04-05 Compute CrossSectionPerVolume (V.Ivanchenko)
00042 // 06-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
00043 // 29-06-06 Fix problem for zero energy incident positron (V.Ivanchenko) 
00044 // 20-10-06 Add theGamma as a member (V.Ivanchenko)
00045 //
00046 //
00047 // Class Description:
00048 //
00049 // Implementation of e+ annihilation into 2 gamma
00050 //
00051 // The secondaries Gamma energies are sampled using the Heitler cross section.
00052 //
00053 // A modified version of the random number techniques of Butcher & Messel
00054 // is used (Nuc Phys 20(1960),15).
00055 //
00056 // GEANT4 internal units.
00057 //
00058 // Note 1: The initial electron is assumed free and at rest.
00059 //
00060 // Note 2: The annihilation processes producing one or more than two photons are
00061 //         ignored, as negligible compared to the two photons process.
00062 
00063 
00064 
00065 //
00066 // -------------------------------------------------------------------
00067 //
00068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00070 
00071 #include "G4eeToTwoGammaModel.hh"
00072 #include "G4PhysicalConstants.hh"
00073 #include "G4SystemOfUnits.hh"
00074 #include "G4TrackStatus.hh"
00075 #include "G4Electron.hh"
00076 #include "G4Positron.hh"
00077 #include "G4Gamma.hh"
00078 #include "Randomize.hh"
00079 #include "G4ParticleChangeForGamma.hh"
00080 
00081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00082 
00083 using namespace std;
00084 
00085 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
00086                                          const G4String& nam)
00087   : G4VEmModel(nam),
00088     pi_rcl2(pi*classic_electr_radius*classic_electr_radius),
00089     isInitialised(false)
00090 {
00091   theGamma = G4Gamma::Gamma();
00092   fParticleChange = 0;
00093 }
00094 
00095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00096 
00097 G4eeToTwoGammaModel::~G4eeToTwoGammaModel()
00098 {}
00099 
00100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00101 
00102 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
00103                                      const G4DataVector&)
00104 {
00105   if(isInitialised) { return; }
00106   fParticleChange = GetParticleChangeForGamma();
00107   isInitialised = true;
00108 }
00109 
00110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00111 
00112 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(
00113                                        const G4ParticleDefinition*,
00114                                        G4double kineticEnergy,
00115                                        G4double, G4double)
00116 {
00117   // Calculates the cross section per electron of annihilation into two photons
00118   // from the Heilter formula.
00119 
00120   G4double ekin  = std::max(eV,kineticEnergy);   
00121 
00122   G4double tau   = ekin/electron_mass_c2;
00123   G4double gam   = tau + 1.0;
00124   G4double gamma2= gam*gam;
00125   G4double bg2   = tau * (tau+2.0);
00126   G4double bg    = sqrt(bg2);
00127 
00128   G4double cross = pi_rcl2*((gamma2+4*gam+1.)*log(gam+bg) - (gam+3.)*bg)
00129                  / (bg2*(gam+1.));
00130   return cross;  
00131 }
00132 
00133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00134 
00135 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
00136                                     const G4ParticleDefinition* p,
00137                                     G4double kineticEnergy, G4double Z,
00138                                     G4double, G4double, G4double)
00139 {
00140   // Calculates the cross section per atom of annihilation into two photons
00141   
00142   G4double cross = Z*ComputeCrossSectionPerElectron(p,kineticEnergy);
00143   return cross;  
00144 }
00145 
00146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00147 
00148 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
00149                                         const G4Material* material,
00150                                         const G4ParticleDefinition* p,
00151                                               G4double kineticEnergy,
00152                                               G4double, G4double)
00153 {
00154   // Calculates the cross section per volume of annihilation into two photons
00155   
00156   G4double eDensity = material->GetElectronDensity();
00157   G4double cross = eDensity*ComputeCrossSectionPerElectron(p,kineticEnergy);
00158   return cross;
00159 }
00160 
00161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00162 
00163 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00164                                             const G4MaterialCutsCouple*,
00165                                             const G4DynamicParticle* dp,
00166                                             G4double,
00167                                             G4double)
00168 {
00169   G4double PositKinEnergy = dp->GetKineticEnergy();
00170   G4DynamicParticle *aGamma1, *aGamma2;
00171    
00172   // Case at rest
00173   if(PositKinEnergy == 0.0) {
00174     G4double cost = 2.*G4UniformRand()-1.;
00175     G4double sint = sqrt((1. - cost)*(1. + cost));
00176     G4double phi  = twopi * G4UniformRand();
00177     G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
00178     phi = twopi * G4UniformRand();
00179     G4ThreeVector pol(cos(phi), sin(phi), 0.0);
00180     pol.rotateUz(dir);
00181     aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
00182     aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
00183     aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
00184     aGamma1->SetPolarization(-pol.x(),-pol.y(),-pol.z());
00185 
00186   } else {
00187 
00188     G4ThreeVector PositDirection = dp->GetMomentumDirection();
00189 
00190     G4double tau     = PositKinEnergy/electron_mass_c2;
00191     G4double gam     = tau + 1.0;
00192     G4double tau2    = tau + 2.0;
00193     G4double sqgrate = sqrt(tau/tau2)*0.5;
00194     G4double sqg2m1  = sqrt(tau*tau2);
00195 
00196     // limits of the energy sampling
00197     G4double epsilmin = 0.5 - sqgrate;
00198     G4double epsilmax = 0.5 + sqgrate;
00199     G4double epsilqot = epsilmax/epsilmin;
00200 
00201     //
00202     // sample the energy rate of the created gammas
00203     //
00204     G4double epsil, greject;
00205 
00206     do {
00207       epsil = epsilmin*pow(epsilqot,G4UniformRand());
00208       greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
00209     } while( greject < G4UniformRand() );
00210 
00211     //
00212     // scattered Gamma angles. ( Z - axis along the parent positron)
00213     //
00214 
00215     G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
00216     if(std::abs(cost) > 1.0) {
00217       G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost
00218              << " positron Ekin(MeV)= " << PositKinEnergy
00219              << " gamma epsil= " << epsil
00220              << G4endl;
00221       if(cost > 1.0) cost = 1.0;
00222       else cost = -1.0; 
00223     }
00224     G4double sint = sqrt((1.+cost)*(1.-cost));
00225     G4double phi  = twopi * G4UniformRand();
00226 
00227     //
00228     // kinematic of the created pair
00229     //
00230 
00231     G4double TotalAvailableEnergy = PositKinEnergy + 2.0*electron_mass_c2;
00232     G4double Phot1Energy = epsil*TotalAvailableEnergy;
00233 
00234     G4ThreeVector Phot1Direction(sint*cos(phi), sint*sin(phi), cost);
00235     Phot1Direction.rotateUz(PositDirection);
00236     aGamma1 = new G4DynamicParticle (theGamma,Phot1Direction, Phot1Energy);
00237     phi = twopi * G4UniformRand();
00238     G4ThreeVector pol1(cos(phi), sin(phi), 0.0);
00239     pol1.rotateUz(Phot1Direction);
00240     aGamma1->SetPolarization(pol1.x(),pol1.y(),pol1.z());
00241 
00242     G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
00243     G4double PositP= sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
00244     G4ThreeVector dir = PositDirection*PositP - Phot1Direction*Phot1Energy;
00245     G4ThreeVector Phot2Direction = dir.unit();
00246 
00247     // create G4DynamicParticle object for the particle2
00248     aGamma2 = new G4DynamicParticle (theGamma,Phot2Direction, Phot2Energy);
00249 
00251     aGamma2->SetPolarization(-pol1.x(),-pol1.y(),-pol1.z());
00252   }
00253   /*
00254     G4cout << "Annihilation in fly: e0= " << PositKinEnergy
00255     << " m= " << electron_mass_c2
00256     << " e1= " << Phot1Energy 
00257     << " e2= " << Phot2Energy << " dir= " <<  dir 
00258     << " -> " << Phot1Direction << " " 
00259     << Phot2Direction << G4endl;
00260   */
00261  
00262   vdp->push_back(aGamma1);
00263   vdp->push_back(aGamma2);
00264   fParticleChange->SetProposedKineticEnergy(0.);
00265   fParticleChange->ProposeTrackStatus(fStopAndKill);
00266 }
00267 
00268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

Generated on Mon May 27 17:48:06 2013 for Geant4 by  doxygen 1.4.7