source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilationModel.cc@ 1197

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 10.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4PenelopeAnnihilationModel.cc,v 1.4 2009/06/10 13:32:36 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 29 Oct 2008 L Pandola Migration from process to model
34// 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
35// - apply internal high-energy limit only in constructor
36// - do not apply low-energy limit (default is 0)
37// - do not use G4ElementSelector
38
39#include "G4PenelopeAnnihilationModel.hh"
40#include "G4ParticleDefinition.hh"
41#include "G4MaterialCutsCouple.hh"
42#include "G4ProductionCutsTable.hh"
43#include "G4DynamicParticle.hh"
44#include "G4Gamma.hh"
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48
49G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*,
50 const G4String& nam)
51 :G4VEmModel(nam),isInitialised(false)
52{
53 fIntrinsicLowEnergyLimit = 0.0;
54 fIntrinsicHighEnergyLimit = 100.0*GeV;
55 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
56 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
57
58 //Calculate variable that will be used later on
59 fPielr2 = pi*classic_electr_radius*classic_electr_radius;
60
61 verboseLevel= 0;
62 // Verbosity scale:
63 // 0 = nothing
64 // 1 = warning for energy non-conservation
65 // 2 = details of energy budget
66 // 3 = calculation of cross sections, file openings, sampling of atoms
67 // 4 = entering in methods
68
69}
70
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72
73G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
74{;}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*,
79 const G4DataVector&)
80{
81 if (verboseLevel > 3)
82 G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
83
84 if(verboseLevel > 0) {
85 G4cout << "Penelope Annihilation model is initialized " << G4endl
86 << "Energy range: "
87 << LowEnergyLimit() / keV << " keV - "
88 << HighEnergyLimit() / GeV << " GeV"
89 << G4endl;
90 }
91
92 if(isInitialised) return;
93 fParticleChange = GetParticleChangeForGamma();
94 isInitialised = true;
95}
96
97//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98
99G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
100 const G4ParticleDefinition*,
101 G4double energy,
102 G4double Z, G4double,
103 G4double, G4double)
104{
105 if (verboseLevel > 3)
106 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" <<
107 G4endl;
108
109 G4double cs = Z*ComputeCrossSectionPerElectron(energy);
110
111 if (verboseLevel > 2)
112 G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z <<
113 " = " << cs/barn << " barn" << G4endl;
114 return cs;
115}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
120 const G4MaterialCutsCouple*,
121 const G4DynamicParticle* aDynamicPositron,
122 G4double,
123 G4double)
124{
125 //
126 // Penelope model to sample final state for positron annihilation.
127 // Target eletrons are assumed to be free and at rest. Binding effects enabling
128 // one-photon annihilation are neglected.
129 // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV
130 // and isotropic angular distribution.
131 // For annihilation in flight, it is used the theory from
132 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
133 // The two photons can have different energy. The efficiency of the sampling algorithm
134 // of the photon energy from the dSigma/dE distribution is practically 100% for
135 // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy
136 // of about 10 MeV.
137 // The angle theta is kinematically linked to the photon energy, to ensure momentum
138 // conservation. The angle phi is sampled isotropically for the first gamma.
139 //
140 if (verboseLevel > 3)
141 G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
142
143 G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
144
145 // kill primary
146 fParticleChange->SetProposedKineticEnergy(0.);
147 fParticleChange->ProposeTrackStatus(fStopAndKill);
148
149 if (kineticEnergy == 0.0)
150 {
151 //Old AtRestDoIt
152 G4double cosTheta = -1.0+2.0*G4UniformRand();
153 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
154 G4double phi = twopi*G4UniformRand();
155 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
156 G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
157 direction, electron_mass_c2);
158 G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
159 -direction, electron_mass_c2);
160
161 fvect->push_back(firstGamma);
162 fvect->push_back(secondGamma);
163 return;
164 }
165
166 //This is the "PostStep" case (annihilation in flight)
167 G4ParticleMomentum positronDirection =
168 aDynamicPositron->GetMomentumDirection();
169 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
170 G4double gamma21 = std::sqrt(gamma*gamma-1);
171 G4double ani = 1.0+gamma;
172 G4double chimin = 1.0/(ani+gamma21);
173 G4double rchi = (1.0-chimin)/chimin;
174 G4double gt0 = ani*ani-2.0;
175 G4double test=0.0;
176 G4double epsilon = 0;
177 do{
178 epsilon = chimin*std::pow(rchi,G4UniformRand());
179 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
180 test = G4UniformRand()*gt0-reject;
181 }while(test>0);
182
183 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
184 G4double photon1Energy = epsilon*totalAvailableEnergy;
185 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
186 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
187 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
188
189 //G4double localEnergyDeposit = 0.;
190
191 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
192 G4double phi1 = twopi * G4UniformRand();
193 G4double dirx1 = sinTheta1 * std::cos(phi1);
194 G4double diry1 = sinTheta1 * std::sin(phi1);
195 G4double dirz1 = cosTheta1;
196
197 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
198 G4double phi2 = phi1+pi;
199 G4double dirx2 = sinTheta2 * std::cos(phi2);
200 G4double diry2 = sinTheta2 * std::sin(phi2);
201 G4double dirz2 = cosTheta2;
202
203 G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
204 photon1Direction.rotateUz(positronDirection);
205 // create G4DynamicParticle object for the particle1
206 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
207 photon1Direction,
208 photon1Energy);
209 fvect->push_back(aParticle1);
210
211 G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
212 photon2Direction.rotateUz(positronDirection);
213 // create G4DynamicParticle object for the particle2
214 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
215 photon2Direction,
216 photon2Energy);
217 fvect->push_back(aParticle2);
218
219 if (verboseLevel > 1)
220 {
221 G4cout << "-----------------------------------------------------------" << G4endl;
222 G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
223 G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
224 G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
225 G4cout << "-----------------------------------------------------------" << G4endl;
226 G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
227 G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
228 G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV <<
229 " keV" << G4endl;
230 G4cout << "-----------------------------------------------------------" << G4endl;
231 }
232 if (verboseLevel > 0)
233 {
234 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
235 if (energyDiff > 0.05*keV)
236 G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
237 (photon1Energy+photon2Energy)/keV <<
238 " keV (final) vs. " <<
239 totalAvailableEnergy/keV << " keV (initial)" << G4endl;
240 }
241 return;
242}
243
244//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245
246G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
247{
248 //
249 // Penelope model to calculate cross section for positron annihilation.
250 // The annihilation cross section per electron is calculated according
251 // to the Heitler formula
252 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
253 // in the assumptions of electrons free and at rest.
254 //
255 G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
256 G4double gamma2 = gamma*gamma;
257 G4double f2 = gamma2-1.0;
258 G4double f1 = std::sqrt(f2);
259 G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
260 - (gamma+3.0)/f1)/(gamma+1.0);
261 return crossSection;
262}
Note: See TracBrowser for help on using the repository browser.