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

Last change on this file since 1201 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 11.4 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 }
[1055]66
67 G4cout << G4endl;
68 G4cout << "*******************************************************************************" << G4endl;
69 G4cout << "*******************************************************************************" << G4endl;
70 G4cout << " The class G4PenelopeAnnihilation is NOT SUPPORTED ANYMORE. " << G4endl;
71 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
72 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
73 G4cout << "*******************************************************************************" << G4endl;
74 G4cout << "*******************************************************************************" << G4endl;
75 G4cout << G4endl;
[819]76}
77
78// destructor
79G4PenelopeAnnihilation::~G4PenelopeAnnihilation()
80{
81 if (meanFreePathTable) {
82 meanFreePathTable->clearAndDestroy();
83 delete meanFreePathTable;
84 }
85}
86
87
88void G4PenelopeAnnihilation::BuildPhysicsTable(const G4ParticleDefinition& )
89{
90 G4double lowEdgeEnergy, value;
91 G4PhysicsLogVector* dataVector;
92
93 // Build mean free path table for the e+e- annihilation
94
95 if (meanFreePathTable) {
96 meanFreePathTable->clearAndDestroy(); delete meanFreePathTable;}
97
98 meanFreePathTable =
99 new G4PhysicsTable(G4Material::GetNumberOfMaterials());
100 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
101 G4Material* material;
102
103 for (size_t j=0;j<G4Material::GetNumberOfMaterials();j++)
104 {
105 //create physics vector then fill it ....
106 dataVector = new G4PhysicsLogVector(lowEnergyLimit,
107 highEnergyLimit,nBins);
108 material = (*theMaterialTable)[j];
109 const G4ElementVector* theElementVector = material->GetElementVector();
110 const G4double* NbOfAtomsPerVolume =
111 material->GetVecNbOfAtomsPerVolume();
112 for (G4int i=0;i<nBins;i++)
113 {
114 lowEdgeEnergy = dataVector->GetLowEdgeEnergy(i);
115 G4double sigma=0.0 ;
116 for (size_t elm=0;elm<material->GetNumberOfElements();elm++)
117 {
118 G4double Z = (*theElementVector)[elm]->GetZ();
119 sigma += NbOfAtomsPerVolume[elm] * Z *
120 calculateCrossSectionPerElectron(lowEdgeEnergy);
121 }
122 if (sigma > DBL_MIN)
123 {
124 value = 1.0/sigma;
125 }
126 else
127 {
128 value = DBL_MAX;
129 }
130 dataVector->PutValue(i,value);
131 }
132 meanFreePathTable->insertAt(j,dataVector);
133 }
134 PrintInfoDefinition();
135}
136
137G4double G4PenelopeAnnihilation::calculateCrossSectionPerElectron
138 (G4double ene)
139{
140 //Heitler dcs formula for annihilation with free electrons at rest
141 G4double crossSection=0.0;
142 G4double gamma = 1.0+std::max(ene,1.0*eV)/electron_mass_c2;
143 G4double gamma2 = gamma*gamma;
144 G4double f2 = gamma2-1.0;
145 G4double f1 = std::sqrt(f2);
146 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius;
147 crossSection = pielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
148 - (gamma+3.0)/f1)/(gamma+1.0);
149 return crossSection;
150}
151
152G4VParticleChange* G4PenelopeAnnihilation::PostStepDoIt(const G4Track& aTrack,
153 const G4Step& )
154{
155 aParticleChange.Initialize(aTrack);
156
157 const G4DynamicParticle* incidentPositron = aTrack.GetDynamicParticle();
158 G4double kineticEnergy = incidentPositron->GetKineticEnergy();
159 G4ParticleMomentum positronDirection =
160 incidentPositron->GetMomentumDirection();
161
162 // Do not make anything if particle is stopped, the annihilation then
163 // should be performed by the AtRestDoIt!
164 if (aTrack.GetTrackStatus() == fStopButAlive) return &aParticleChange;
165
166 //G4cout << "Sono nel PostStep" << G4endl;
167
168 //Annihilation in flight
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 epsilon=0.0, reject=0.0, test=0.0;
176 do{
177 epsilon = chimin*std::pow(rchi,G4UniformRand());
178 reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
179 test = G4UniformRand()*gt0-reject;
180 }while(test>0);
181
182 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
183 G4double photon1Energy = epsilon*totalAvailableEnergy;
184 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
185 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
186 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
187
188 aParticleChange.SetNumberOfSecondaries(2);
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 if (photon1Energy > cutForLowEnergySecondaryPhotons) {
204 G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
205 photon1Direction.rotateUz(positronDirection);
206 // create G4DynamicParticle object for the particle1
207 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
208 photon1Direction,
209 photon1Energy);
210 aParticleChange.AddSecondary(aParticle1);
211 }
212 else localEnergyDeposit += photon1Energy;
213
214 if (photon2Energy > cutForLowEnergySecondaryPhotons) {
215 G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
216 photon2Direction.rotateUz(positronDirection);
217 // create G4DynamicParticle object for the particle2
218 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
219 photon2Direction,
220 photon2Energy);
221 aParticleChange.AddSecondary(aParticle2);
222 }
223 else localEnergyDeposit += photon2Energy;
224
225 aParticleChange.ProposeLocalEnergyDeposit(localEnergyDeposit);
226
227 aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
228 aParticleChange.ProposeEnergy(0.);
229 aParticleChange.ProposeTrackStatus(fStopAndKill);
230
231 return &aParticleChange;
232}
233
234
235G4VParticleChange* G4PenelopeAnnihilation::AtRestDoIt(const G4Track& aTrack,
236 const G4Step& )
237{
238 aParticleChange.Initialize(aTrack);
239 aParticleChange.SetNumberOfSecondaries(2);
240 G4double cosTheta = -1.0+2.0*G4UniformRand();
241 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
242 G4double phi = twopi*G4UniformRand();
243 //G4cout << "cosTheta: " << cosTheta << " sinTheta: " << sinTheta << G4endl;
244 //G4cout << "phi: " << phi << G4endl;
245 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
246 aParticleChange.AddSecondary(new G4DynamicParticle (G4Gamma::Gamma(),
247 direction, electron_mass_c2) );
248 aParticleChange.AddSecondary(new G4DynamicParticle (G4Gamma::Gamma(),
249 -direction, electron_mass_c2) );
250
251 aParticleChange.ProposeLocalEnergyDeposit(0.);
252
253 // Kill the incident positron
254 //
255 aParticleChange.ProposeTrackStatus(fStopAndKill);
256
257 return &aParticleChange;
258}
259
260void G4PenelopeAnnihilation::PrintInfoDefinition()
261{
262 G4String comments = "Total cross section from Heilter formula"
263 "(annihilation into 2 photons).";
264 comments += "\n Gamma energies sampled according Heitler";
265 comments += "\n It can be used for positrons";
266 comments += " in the energy range [250eV,100GeV].";
267 G4cout << G4endl << GetProcessName() << ": " << comments
268 << G4endl;
269}
270
271G4bool G4PenelopeAnnihilation::IsApplicable(
272 const G4ParticleDefinition& particle)
273{
274 return ( &particle == G4Positron::Positron() );
275}
276
277
278G4double G4PenelopeAnnihilation::GetMeanFreePath(const G4Track& aTrack,
279 G4double,
280 G4ForceCondition*)
281{
282 const G4DynamicParticle* incidentPositron = aTrack.GetDynamicParticle();
283 G4double kineticEnergy = incidentPositron->GetKineticEnergy();
284 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
285 const G4Material* material = couple->GetMaterial();
286 G4int materialIndex = material->GetIndex();
287 if (kineticEnergy<lowEnergyLimit)
288 {
289 kineticEnergy = lowEnergyLimit;
290 }
291
292 G4double meanFreePath;
293 G4bool isOutRange ;
294
295 if (kineticEnergy>highEnergyLimit)
296 {
297 meanFreePath = DBL_MAX;
298 }
299 else
300 {
301 meanFreePath = (*meanFreePathTable)(materialIndex)->
302 GetValue(kineticEnergy,isOutRange);
303 }
304
305 return meanFreePath;
306}
307
308G4double G4PenelopeAnnihilation::GetMeanLifeTime(const G4Track&,
309 G4ForceCondition*)
310
311{
312 return 0.0;
313}
Note: See TracBrowser for help on using the repository browser.