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