| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Author: Luciano Pandola 27 // Author: Luciano Pandola
28 // 28 //
29 // History: 29 // History:
30 // -------- 30 // --------
31 // 29 Oct 2008 L Pandola Migration from p 31 // 29 Oct 2008 L Pandola Migration from process to model
32 // 15 Apr 2009 V Ivanchenko Cleanup initiali 32 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of
33 // secondaries: 33 // secondaries:
34 // - apply internal high-ener 34 // - apply internal high-energy limit only in constructor
35 // - do not apply low-energy 35 // - do not apply low-energy limit (default is 0)
36 // - do not use G4ElementSele 36 // - do not use G4ElementSelector
37 // 02 Oct 2013 L.Pandola Migration to MT 37 // 02 Oct 2013 L.Pandola Migration to MT
38 38
39 #include "G4PenelopeAnnihilationModel.hh" 39 #include "G4PenelopeAnnihilationModel.hh"
40 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
42 #include "G4ParticleDefinition.hh" 42 #include "G4ParticleDefinition.hh"
43 #include "G4MaterialCutsCouple.hh" 43 #include "G4MaterialCutsCouple.hh"
44 #include "G4ProductionCutsTable.hh" 44 #include "G4ProductionCutsTable.hh"
45 #include "G4DynamicParticle.hh" 45 #include "G4DynamicParticle.hh"
46 #include "G4Gamma.hh" 46 #include "G4Gamma.hh"
47 47
48 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 49
50 G4double G4PenelopeAnnihilationModel::fPielr2 50 G4double G4PenelopeAnnihilationModel::fPielr2 = 0;
51 51
52 G4PenelopeAnnihilationModel::G4PenelopeAnnihil 52 G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition* part,
53 c 53 const G4String& nam)
54 :G4VEmModel(nam),fParticleChange(nullptr),fP << 54 :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false)
55 { 55 {
56 fIntrinsicLowEnergyLimit = 0.0; 56 fIntrinsicLowEnergyLimit = 0.0;
57 fIntrinsicHighEnergyLimit = 100.0*GeV; 57 fIntrinsicHighEnergyLimit = 100.0*GeV;
>> 58 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
58 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 59 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
59 60
60 if (part) 61 if (part)
61 SetParticle(part); 62 SetParticle(part);
62 63
63 //Calculate variable that will be used later 64 //Calculate variable that will be used later on
64 fPielr2 = pi*classic_electr_radius*classic_e 65 fPielr2 = pi*classic_electr_radius*classic_electr_radius;
65 66
66 fVerboseLevel= 0; << 67 verboseLevel= 0;
67 // Verbosity scale: 68 // Verbosity scale:
68 // 0 = nothing 69 // 0 = nothing
69 // 1 = warning for energy non-conservation 70 // 1 = warning for energy non-conservation
70 // 2 = details of energy budget 71 // 2 = details of energy budget
71 // 3 = calculation of cross sections, file o 72 // 3 = calculation of cross sections, file openings, sampling of atoms
72 // 4 = entering in methods 73 // 4 = entering in methods
>> 74
73 } 75 }
74 76
75 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 78
77 G4PenelopeAnnihilationModel::~G4PenelopeAnnihi 79 G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
78 {;} 80 {;}
79 81
80 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 83
82 void G4PenelopeAnnihilationModel::Initialise(c 84 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition* part,
83 const G4DataVector&) 85 const G4DataVector&)
84 { 86 {
85 if (fVerboseLevel > 3) << 87 if (verboseLevel > 3)
86 G4cout << "Calling G4PenelopeAnnihilationM 88 G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
87 SetParticle(part); 89 SetParticle(part);
88 90
89 if (IsMaster() && part == fParticle) 91 if (IsMaster() && part == fParticle)
90 { 92 {
91 93
92 if(fVerboseLevel > 0) { << 94 if(verboseLevel > 0) {
93 G4cout << "Penelope Annihilation model is in 95 G4cout << "Penelope Annihilation model is initialized " << G4endl
94 << "Energy range: " 96 << "Energy range: "
95 << LowEnergyLimit() / keV << " keV - 97 << LowEnergyLimit() / keV << " keV - "
96 << HighEnergyLimit() / GeV << " GeV" 98 << HighEnergyLimit() / GeV << " GeV"
97 << G4endl; 99 << G4endl;
98 } 100 }
99 } 101 }
100 102
101 if(fIsInitialised) return; << 103 if(isInitialised) return;
102 fParticleChange = GetParticleChangeForGamma( 104 fParticleChange = GetParticleChangeForGamma();
103 fIsInitialised = true; << 105 isInitialised = true;
104 } 106 }
105 107
106 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 void G4PenelopeAnnihilationModel::InitialiseLo 109 void G4PenelopeAnnihilationModel::InitialiseLocal(const G4ParticleDefinition* part,
108 G4VEmModel* masterModel) 110 G4VEmModel* masterModel)
109 { 111 {
110 if (fVerboseLevel > 3) << 112 if (verboseLevel > 3)
111 G4cout << "Calling G4PenelopeAnnihilationM 113 G4cout << "Calling G4PenelopeAnnihilationModel::InitialiseLocal()" << G4endl;
112 114
113 // 115 //
114 //Check that particle matches: one might hav 116 //Check that particle matches: one might have multiple master models (e.g.
115 //for e+ and e-). 117 //for e+ and e-).
116 // 118 //
117 if (part == fParticle) 119 if (part == fParticle)
118 { 120 {
119 //Get the const table pointers from the 121 //Get the const table pointers from the master to the workers
120 const G4PenelopeAnnihilationModel* theMo 122 const G4PenelopeAnnihilationModel* theModel =
121 static_cast<G4PenelopeAnnihilationMode 123 static_cast<G4PenelopeAnnihilationModel*> (masterModel);
122 124
123 //Same verbosity for all workers, as the 125 //Same verbosity for all workers, as the master
124 fVerboseLevel = theModel->fVerboseLevel; << 126 verboseLevel = theModel->verboseLevel;
125 } 127 }
126 } 128 }
127 129
128 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 131
130 G4double G4PenelopeAnnihilationModel::ComputeC 132 G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
131 const G 133 const G4ParticleDefinition*,
132 G 134 G4double energy,
133 G 135 G4double Z, G4double,
134 G 136 G4double, G4double)
135 { 137 {
136 if (fVerboseLevel > 3) << 138 if (verboseLevel > 3)
137 G4cout << "Calling ComputeCrossSectionPerA 139 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" <<
138 G4endl; 140 G4endl;
139 141
140 G4double cs = Z*ComputeCrossSectionPerElectr 142 G4double cs = Z*ComputeCrossSectionPerElectron(energy);
141 143
142 if (fVerboseLevel > 2) << 144 if (verboseLevel > 2)
143 G4cout << "Annihilation cross Section at " 145 G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z <<
144 " = " << cs/barn << " barn" << G4endl; 146 " = " << cs/barn << " barn" << G4endl;
145 return cs; 147 return cs;
146 } 148 }
147 149
148 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 151
150 void G4PenelopeAnnihilationModel::SampleSecond 152 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
151 const G4MaterialCutsCouple*, 153 const G4MaterialCutsCouple*,
152 const G4DynamicParticle* aDyna 154 const G4DynamicParticle* aDynamicPositron,
153 G4double, 155 G4double,
154 G4double) 156 G4double)
155 { 157 {
156 // 158 //
157 // Penelope model to sample final state for 159 // Penelope model to sample final state for positron annihilation.
158 // Target eletrons are assumed to be free an 160 // Target eletrons are assumed to be free and at rest. Binding effects enabling
159 // one-photon annihilation are neglected. 161 // one-photon annihilation are neglected.
160 // For annihilation at rest, two back-to-bac 162 // For annihilation at rest, two back-to-back photons are emitted, having energy of 511 keV
161 // and isotropic angular distribution. 163 // and isotropic angular distribution.
162 // For annihilation in flight, it is used th 164 // For annihilation in flight, it is used the theory from
163 // W. Heitler, The quantum theory of radiat 165 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
164 // The two photons can have different energy 166 // The two photons can have different energy. The efficiency of the sampling algorithm
165 // of the photon energy from the dSigma/dE d 167 // of the photon energy from the dSigma/dE distribution is practically 100% for
166 // positrons of kinetic energy < 10 keV. It 168 // positrons of kinetic energy < 10 keV. It reaches a minimum (about 80%) at energy
167 // of about 10 MeV. 169 // of about 10 MeV.
168 // The angle theta is kinematically linked t 170 // The angle theta is kinematically linked to the photon energy, to ensure momentum
169 // conservation. The angle phi is sampled is 171 // conservation. The angle phi is sampled isotropically for the first gamma.
170 // 172 //
171 if (fVerboseLevel > 3) << 173 if (verboseLevel > 3)
172 G4cout << "Calling SamplingSecondaries() o 174 G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
173 175
174 G4double kineticEnergy = aDynamicPositron->G 176 G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
175 177
176 // kill primary 178 // kill primary
177 fParticleChange->SetProposedKineticEnergy(0. 179 fParticleChange->SetProposedKineticEnergy(0.);
178 fParticleChange->ProposeTrackStatus(fStopAnd 180 fParticleChange->ProposeTrackStatus(fStopAndKill);
179 181
180 if (kineticEnergy == 0.0) 182 if (kineticEnergy == 0.0)
181 { 183 {
182 //Old AtRestDoIt 184 //Old AtRestDoIt
183 G4double cosTheta = -1.0+2.0*G4UniformRa 185 G4double cosTheta = -1.0+2.0*G4UniformRand();
184 G4double sinTheta = std::sqrt(1.0-cosThe 186 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
185 G4double phi = twopi*G4UniformRand(); 187 G4double phi = twopi*G4UniformRand();
186 G4ThreeVector direction (sinTheta*std::c 188 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
187 G4DynamicParticle* firstGamma = new G4Dy 189 G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
188 direction, electron_mass_c2 190 direction, electron_mass_c2);
189 G4DynamicParticle* secondGamma = new G4D 191 G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
190 -direction, electron_mass_ 192 -direction, electron_mass_c2);
191 193
192 fvect->push_back(firstGamma); 194 fvect->push_back(firstGamma);
193 fvect->push_back(secondGamma); 195 fvect->push_back(secondGamma);
194 return; 196 return;
195 } 197 }
196 198
197 //This is the "PostStep" case (annihilation 199 //This is the "PostStep" case (annihilation in flight)
198 G4ParticleMomentum positronDirection = 200 G4ParticleMomentum positronDirection =
199 aDynamicPositron->GetMomentumDirection(); 201 aDynamicPositron->GetMomentumDirection();
200 G4double gamma = 1.0 + std::max(kineticEnerg 202 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
201 G4double gamma21 = std::sqrt(gamma*gamma-1); 203 G4double gamma21 = std::sqrt(gamma*gamma-1);
202 G4double ani = 1.0+gamma; 204 G4double ani = 1.0+gamma;
203 G4double chimin = 1.0/(ani+gamma21); 205 G4double chimin = 1.0/(ani+gamma21);
204 G4double rchi = (1.0-chimin)/chimin; 206 G4double rchi = (1.0-chimin)/chimin;
205 G4double gt0 = ani*ani-2.0; 207 G4double gt0 = ani*ani-2.0;
206 G4double test=0.0; 208 G4double test=0.0;
207 G4double epsilon = 0; 209 G4double epsilon = 0;
208 do{ 210 do{
209 epsilon = chimin*std::pow(rchi,G4UniformRa 211 epsilon = chimin*std::pow(rchi,G4UniformRand());
210 G4double reject = ani*ani*(1.0-epsilon)+2. 212 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
211 test = G4UniformRand()*gt0-reject; 213 test = G4UniformRand()*gt0-reject;
212 }while(test>0); 214 }while(test>0);
213 215
214 G4double totalAvailableEnergy = kineticEnerg 216 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
215 G4double photon1Energy = epsilon*totalAvaila 217 G4double photon1Energy = epsilon*totalAvailableEnergy;
216 G4double photon2Energy = (1.0-epsilon)*total 218 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
217 G4double cosTheta1 = (ani-1.0/epsilon)/gamma 219 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
218 G4double cosTheta2 = (ani-1.0/(1.0-epsilon)) 220 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
219 221
>> 222 //G4double localEnergyDeposit = 0.;
>> 223
220 G4double sinTheta1 = std::sqrt(1.-cosTheta1* 224 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
221 G4double phi1 = twopi * G4UniformRand(); 225 G4double phi1 = twopi * G4UniformRand();
222 G4double dirx1 = sinTheta1 * std::cos(phi1); 226 G4double dirx1 = sinTheta1 * std::cos(phi1);
223 G4double diry1 = sinTheta1 * std::sin(phi1); 227 G4double diry1 = sinTheta1 * std::sin(phi1);
224 G4double dirz1 = cosTheta1; 228 G4double dirz1 = cosTheta1;
225 229
226 G4double sinTheta2 = std::sqrt(1.-cosTheta2* 230 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
227 G4double phi2 = phi1+pi; 231 G4double phi2 = phi1+pi;
228 G4double dirx2 = sinTheta2 * std::cos(phi2); 232 G4double dirx2 = sinTheta2 * std::cos(phi2);
229 G4double diry2 = sinTheta2 * std::sin(phi2); 233 G4double diry2 = sinTheta2 * std::sin(phi2);
230 G4double dirz2 = cosTheta2; 234 G4double dirz2 = cosTheta2;
231 235
232 G4ThreeVector photon1Direction (dirx1,diry1, 236 G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
233 photon1Direction.rotateUz(positronDirection) 237 photon1Direction.rotateUz(positronDirection);
234 // create G4DynamicParticle object for the p 238 // create G4DynamicParticle object for the particle1
235 G4DynamicParticle* aParticle1= new G4Dynamic 239 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
236 photon1Direction, 240 photon1Direction,
237 photon1Energy); 241 photon1Energy);
238 fvect->push_back(aParticle1); 242 fvect->push_back(aParticle1);
239 243
240 G4ThreeVector photon2Direction(dirx2,diry2,d 244 G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
241 photon2Direction.rotateUz(positronDirection) 245 photon2Direction.rotateUz(positronDirection);
242 // create G4DynamicParticle object for the p << 246 // create G4DynamicParticle object for the particle2
243 G4DynamicParticle* aParticle2= new G4Dynamic 247 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
244 photon2Direction, 248 photon2Direction,
245 photon2Energy); 249 photon2Energy);
246 fvect->push_back(aParticle2); 250 fvect->push_back(aParticle2);
247 251
248 if (fVerboseLevel > 1) << 252 if (verboseLevel > 1)
249 { 253 {
250 G4cout << "----------------------------- 254 G4cout << "-----------------------------------------------------------" << G4endl;
251 G4cout << "Energy balance from G4Penelop 255 G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
252 G4cout << "Kinetic positron energy: " << 256 G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
253 G4cout << "Total available energy: " << 257 G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
254 G4cout << "----------------------------- 258 G4cout << "-----------------------------------------------------------" << G4endl;
255 G4cout << "Photon energy 1: " << photon1 259 G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
256 G4cout << "Photon energy 2: " << photon2 260 G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
257 G4cout << "Total final state: " << (phot 261 G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV <<
258 " keV" << G4endl; 262 " keV" << G4endl;
259 G4cout << "----------------------------- 263 G4cout << "-----------------------------------------------------------" << G4endl;
260 } 264 }
261 if (fVerboseLevel > 0) << 265 if (verboseLevel > 0)
262 { 266 {
263 G4double energyDiff = std::fabs(totalAva 267 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
264 if (energyDiff > 0.05*keV) 268 if (energyDiff > 0.05*keV)
265 G4cout << "Warning from G4PenelopeAnnihilati 269 G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
266 (photon1Energy+photon2Energy)/keV << 270 (photon1Energy+photon2Energy)/keV <<
267 " keV (final) vs. " << 271 " keV (final) vs. " <<
268 totalAvailableEnergy/keV << " keV (initial 272 totalAvailableEnergy/keV << " keV (initial)" << G4endl;
269 } 273 }
270 return; 274 return;
271 } 275 }
272 276
273 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 278
275 G4double G4PenelopeAnnihilationModel:: Compute 279 G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
276 { 280 {
277 // 281 //
278 // Penelope model to calculate cross section 282 // Penelope model to calculate cross section for positron annihilation.
279 // The annihilation cross section per electr 283 // The annihilation cross section per electron is calculated according
280 // to the Heitler formula 284 // to the Heitler formula
281 // W. Heitler, The quantum theory of radiat 285 // W. Heitler, The quantum theory of radiation, Oxford University Press (1954)
282 // in the assumptions of electrons free and 286 // in the assumptions of electrons free and at rest.
283 // 287 //
284 G4double gamma = 1.0+std::max(energy,1.0*eV) 288 G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
285 G4double gamma2 = gamma*gamma; 289 G4double gamma2 = gamma*gamma;
286 G4double f2 = gamma2-1.0; 290 G4double f2 = gamma2-1.0;
287 G4double f1 = std::sqrt(f2); 291 G4double f1 = std::sqrt(f2);
288 G4double crossSection = fPielr2*((gamma2+4.0 << 292 G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
289 - (gamma+3.0)/f1)/(gamma+1.0); 293 - (gamma+3.0)/f1)/(gamma+1.0);
290 return crossSection; 294 return crossSection;
291 } 295 }
292 296
293 //....oooOO0OOooo........oooOO0OOooo........oo 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
294 298
295 void G4PenelopeAnnihilationModel::SetParticle( 299 void G4PenelopeAnnihilationModel::SetParticle(const G4ParticleDefinition* p)
296 { 300 {
297 if(!fParticle) { 301 if(!fParticle) {
298 fParticle = p; 302 fParticle = p;
299 } 303 }
300 } 304 }
301 305