source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeAnnihilation.cc@ 1036

Last change on this file since 1036 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 10.7 KB
RevLine 
[819]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// Author: Luciano Pandola
27//
28// History:
29// --------
30// 02 Jul 2003 L.Pandola First implementation
31// 16 Mar 2004 L.Pandola Removed unnecessary calls to std::pow(a,b)
32
33#include "G4PenelopeAnnihilation.hh"
34#include "Randomize.hh"
35#include "G4UnitsTable.hh"
36#include "G4PhysicsTable.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4VParticleChange.hh"
39#include "G4Gamma.hh"
40#include "G4Positron.hh"
41#include "G4Track.hh"
42#include "G4Step.hh"
43#include "G4PhysicsLogVector.hh"
44#include "G4ElementTable.hh"
45#include "G4Material.hh"
46#include "G4MaterialCutsCouple.hh"
47
48// constructor
49
50G4PenelopeAnnihilation::G4PenelopeAnnihilation(const G4String& processName)
51 : G4VRestDiscreteProcess (processName),
52 lowEnergyLimit(250*eV),
53 highEnergyLimit(100*GeV),
54 nBins(200),
55 cutForLowEnergySecondaryPhotons(250.0*eV)
56{
57 meanFreePathTable = 0;
58 if (verboseLevel > 0)
59 {
60 G4cout << GetProcessName() << " is created " << G4endl
61 << "Energy range: "
62 << lowEnergyLimit / keV << " keV - "
63 << highEnergyLimit / GeV << " GeV"
64 << G4endl;
65 }
66}
67
68// destructor
69G4PenelopeAnnihilation::~G4PenelopeAnnihilation()
70{
71 if (meanFreePathTable) {
72 meanFreePathTable->clearAndDestroy();
73 delete meanFreePathTable;
74 }
75}
76
77
78void G4PenelopeAnnihilation::BuildPhysicsTable(const G4ParticleDefinition& )
79{
80 G4double lowEdgeEnergy, value;
81 G4PhysicsLogVector* dataVector;
82
83 // Build mean free path table for the e+e- annihilation
84
85 if (meanFreePathTable) {
86 meanFreePathTable->clearAndDestroy(); delete meanFreePathTable;}
87
88 meanFreePathTable =
89 new G4PhysicsTable(G4Material::GetNumberOfMaterials());
90 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
91 G4Material* material;
92
93 for (size_t j=0;j<G4Material::GetNumberOfMaterials();j++)
94 {
95 //create physics vector then fill it ....
96 dataVector = new G4PhysicsLogVector(lowEnergyLimit,
97 highEnergyLimit,nBins);
98 material = (*theMaterialTable)[j];
99 const G4ElementVector* theElementVector = material->GetElementVector();
100 const G4double* NbOfAtomsPerVolume =
101 material->GetVecNbOfAtomsPerVolume();
102 for (G4int i=0;i<nBins;i++)
103 {
104 lowEdgeEnergy = dataVector->GetLowEdgeEnergy(i);
105 G4double sigma=0.0 ;
106 for (size_t elm=0;elm<material->GetNumberOfElements();elm++)
107 {
108 G4double Z = (*theElementVector)[elm]->GetZ();
109 sigma += NbOfAtomsPerVolume[elm] * Z *
110 calculateCrossSectionPerElectron(lowEdgeEnergy);
111 }
112 if (sigma > DBL_MIN)
113 {
114 value = 1.0/sigma;
115 }
116 else
117 {
118 value = DBL_MAX;
119 }
120 dataVector->PutValue(i,value);
121 }
122 meanFreePathTable->insertAt(j,dataVector);
123 }
124 PrintInfoDefinition();
125}
126
127G4double G4PenelopeAnnihilation::calculateCrossSectionPerElectron
128 (G4double ene)
129{
130 //Heitler dcs formula for annihilation with free electrons at rest
131 G4double crossSection=0.0;
132 G4double gamma = 1.0+std::max(ene,1.0*eV)/electron_mass_c2;
133 G4double gamma2 = gamma*gamma;
134 G4double f2 = gamma2-1.0;
135 G4double f1 = std::sqrt(f2);
136 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius;
137 crossSection = pielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
138 - (gamma+3.0)/f1)/(gamma+1.0);
139 return crossSection;
140}
141
142G4VParticleChange* G4PenelopeAnnihilation::PostStepDoIt(const G4Track& aTrack,
143 const G4Step& )
144{
145 aParticleChange.Initialize(aTrack);
146
147 const G4DynamicParticle* incidentPositron = aTrack.GetDynamicParticle();
148 G4double kineticEnergy = incidentPositron->GetKineticEnergy();
149 G4ParticleMomentum positronDirection =
150 incidentPositron->GetMomentumDirection();
151
152 // Do not make anything if particle is stopped, the annihilation then
153 // should be performed by the AtRestDoIt!
154 if (aTrack.GetTrackStatus() == fStopButAlive) return &aParticleChange;
155
156 //G4cout << "Sono nel PostStep" << G4endl;
157
158 //Annihilation in flight
159 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
160 G4double gamma21 = std::sqrt(gamma*gamma-1);
161 G4double ani = 1.0+gamma;
162 G4double chimin = 1.0/(ani+gamma21);
163 G4double rchi = (1.0-chimin)/chimin;
164 G4double gt0 = ani*ani-2.0;
165 G4double epsilon=0.0, reject=0.0, test=0.0;
166 do{
167 epsilon = chimin*std::pow(rchi,G4UniformRand());
168 reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
169 test = G4UniformRand()*gt0-reject;
170 }while(test>0);
171
172 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
173 G4double photon1Energy = epsilon*totalAvailableEnergy;
174 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
175 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
176 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
177
178 aParticleChange.SetNumberOfSecondaries(2);
179 G4double localEnergyDeposit = 0.;
180
181 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
182 G4double phi1 = twopi * G4UniformRand();
183 G4double dirx1 = sinTheta1 * std::cos(phi1);
184 G4double diry1 = sinTheta1 * std::sin(phi1);
185 G4double dirz1 = cosTheta1;
186
187 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
188 G4double phi2 = phi1+pi;
189 G4double dirx2 = sinTheta2 * std::cos(phi2);
190 G4double diry2 = sinTheta2 * std::sin(phi2);
191 G4double dirz2 = cosTheta2;
192
193 if (photon1Energy > cutForLowEnergySecondaryPhotons) {
194 G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
195 photon1Direction.rotateUz(positronDirection);
196 // create G4DynamicParticle object for the particle1
197 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
198 photon1Direction,
199 photon1Energy);
200 aParticleChange.AddSecondary(aParticle1);
201 }
202 else localEnergyDeposit += photon1Energy;
203
204 if (photon2Energy > cutForLowEnergySecondaryPhotons) {
205 G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
206 photon2Direction.rotateUz(positronDirection);
207 // create G4DynamicParticle object for the particle2
208 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
209 photon2Direction,
210 photon2Energy);
211 aParticleChange.AddSecondary(aParticle2);
212 }
213 else localEnergyDeposit += photon2Energy;
214
215 aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit);
216
217 aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
218 aParticleChange.ProposeEnergy(0.);
219 aParticleChange.ProposeTrackStatus(fStopAndKill);
220
221 return &aParticleChange;
222}
223
224
225G4VParticleChange* G4PenelopeAnnihilation::AtRestDoIt(const G4Track& aTrack,
226 const G4Step& )
227{
228 aParticleChange.Initialize(aTrack);
229 aParticleChange.SetNumberOfSecondaries(2);
230 G4double cosTheta = -1.0+2.0*G4UniformRand();
231 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
232 G4double phi = twopi*G4UniformRand();
233 //G4cout << "cosTheta: " << cosTheta << " sinTheta: " << sinTheta << G4endl;
234 //G4cout << "phi: " << phi << G4endl;
235 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
236 aParticleChange.AddSecondary(new G4DynamicParticle (G4Gamma::Gamma(),
237 direction, electron_mass_c2) );
238 aParticleChange.AddSecondary(new G4DynamicParticle (G4Gamma::Gamma(),
239 -direction, electron_mass_c2) );
240
241 aParticleChange.ProposeLocalEnergyDeposit(0.);
242
243 // Kill the incident positron
244 //
245 aParticleChange.ProposeTrackStatus(fStopAndKill);
246
247 return &aParticleChange;
248}
249
250void G4PenelopeAnnihilation::PrintInfoDefinition()
251{
252 G4String comments = "Total cross section from Heilter formula"
253 "(annihilation into 2 photons).";
254 comments += "\n Gamma energies sampled according Heitler";
255 comments += "\n It can be used for positrons";
256 comments += " in the energy range [250eV,100GeV].";
257 G4cout << G4endl << GetProcessName() << ": " << comments
258 << G4endl;
259}
260
261G4bool G4PenelopeAnnihilation::IsApplicable(
262 const G4ParticleDefinition& particle)
263{
264 return ( &particle == G4Positron::Positron() );
265}
266
267
268G4double G4PenelopeAnnihilation::GetMeanFreePath(const G4Track& aTrack,
269 G4double,
270 G4ForceCondition*)
271{
272 const G4DynamicParticle* incidentPositron = aTrack.GetDynamicParticle();
273 G4double kineticEnergy = incidentPositron->GetKineticEnergy();
274 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
275 const G4Material* material = couple->GetMaterial();
276 G4int materialIndex = material->GetIndex();
277 if (kineticEnergy<lowEnergyLimit)
278 {
279 kineticEnergy = lowEnergyLimit;
280 }
281
282 G4double meanFreePath;
283 G4bool isOutRange ;
284
285 if (kineticEnergy>highEnergyLimit)
286 {
287 meanFreePath = DBL_MAX;
288 }
289 else
290 {
291 meanFreePath = (*meanFreePathTable)(materialIndex)->
292 GetValue(kineticEnergy,isOutRange);
293 }
294
295 return meanFreePath;
296}
297
298G4double G4PenelopeAnnihilation::GetMeanLifeTime(const G4Track&,
299 G4ForceCondition*)
300
301{
302 return 0.0;
303}
Note: See TracBrowser for help on using the repository browser.