| 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
28 #include "G4AdjointComptonModel.hh"
29
30 #include "G4AdjointCSManager.hh"
31 #include "G4AdjointElectron.hh"
32 #include "G4AdjointGamma.hh"
33 #include "G4Gamma.hh"
34 #include "G4KleinNishinaCompton.hh"
35 #include "G4ParticleChange.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4TrackStatus.hh"
39 #include "G4VEmProcess.hh"
40
41 //////////////////////////////////////////////
42 G4AdjointComptonModel::G4AdjointComptonModel()
43 : G4VEmAdjointModel("AdjointCompton")
44 {
45 SetApplyCutInRange(false);
46 SetUseMatrix(false);
47 SetUseMatrixPerElement(true);
48 SetUseOnlyOneMatrixForAllElements(true);
49 fAdjEquivDirectPrimPart = G4AdjointGamma::
50 fAdjEquivDirectSecondPart = G4AdjointElectro
51 fDirectPrimaryPart = G4Gamma::Gamma()
52 fSecondPartSameType = false;
53 fDirectModel =
54 new G4KleinNishinaCompton(G4Gamma::Gamma()
55 }
56
57 //////////////////////////////////////////////
58 G4AdjointComptonModel::~G4AdjointComptonModel(
59
60 //////////////////////////////////////////////
61 void G4AdjointComptonModel::SampleSecondaries(
62
63
64 {
65 if(!fUseMatrix)
66 return RapidSampleSecondaries(aTrack, isSc
67
68 // A recall of the compton scattering law:
69 // Egamma2=Egamma1/(1+(Egamma1/E0_electron)(
70 // Therefore Egamma2_max= Egamma2(cos_th=1)
71 // and Egamma2_min= Egamma2(cos_th=-1)
72 // Egamma1/(1+2.(Egamma1/E0_elec
73 const G4DynamicParticle* theAdjointPrimary =
74 G4double adjointPrimKinEnergy = theAdjointPr
75 if(adjointPrimKinEnergy > GetHighEnergyLimit
76 {
77 return;
78 }
79
80 // Sample secondary energy
81 G4double gammaE1;
82 gammaE1 =
83 SampleAdjSecEnergyFromCSMatrix(adjointPrim
84
85 // gammaE2
86 G4double gammaE2 = adjointPrimKinEnergy;
87 if(!isScatProjToProj)
88 gammaE2 = gammaE1 - adjointPrimKinEnergy;
89
90 // Cos th
91 G4double cos_th = 1. + electron_mass_c2 * (1
92 if(!isScatProjToProj)
93 {
94 cos_th =
95 (gammaE1 - gammaE2 * cos_th) / theAdjoin
96 }
97 G4double sin_th = 0.;
98 if(std::abs(cos_th) > 1.)
99 {
100 if(cos_th > 0.)
101 {
102 cos_th = 1.;
103 }
104 else
105 cos_th = -1.;
106 sin_th = 0.;
107 }
108 else
109 sin_th = std::sqrt(1. - cos_th * cos_th);
110
111 // gamma0 momentum
112 G4ThreeVector dir_parallel = theAdjointPrima
113 G4double phi = G4UniformRand()
114 G4ThreeVector gammaMomentum1 =
115 gammaE1 *
116 G4ThreeVector(std::cos(phi) * sin_th, std:
117 gammaMomentum1.rotateUz(dir_parallel);
118
119 // correct the weight of particles before ad
120 CorrectPostStepWeight(fParticleChange, aTrac
121 adjointPrimKinEnergy,
122
123 if(!isScatProjToProj)
124 { // kill the primary and add a secondary
125 fParticleChange->ProposeTrackStatus(fStopA
126 fParticleChange->AddSecondary(
127 new G4DynamicParticle(fAdjEquivDirectPri
128 }
129 else
130 {
131 fParticleChange->ProposeEnergy(gammaE1);
132 fParticleChange->ProposeMomentumDirection(
133 }
134 }
135
136 //////////////////////////////////////////////
137 void G4AdjointComptonModel::RapidSampleSeconda
138 const G4Track& aTrack, G4bool isScatProjToPr
139 G4ParticleChange* fParticleChange)
140 {
141 const G4DynamicParticle* theAdjointPrimary =
142 DefineCurrentMaterial(aTrack.GetMaterialCuts
143
144 G4double adjointPrimKinEnergy = theAdjointPr
145
146 if(adjointPrimKinEnergy > GetHighEnergyLimit
147 {
148 return;
149 }
150
151 G4double diffCSUsed =
152 0.1 * fCurrentMaterial->GetElectronDensity
153 G4double gammaE1 = 0.;
154 G4double gammaE2 = 0.;
155 if(!isScatProjToProj)
156 {
157 G4double Emax = GetSecondAdjEnergyMaxForPr
158 G4double Emin = GetSecondAdjEnergyMinForPr
159 if(Emin >= Emax)
160 return;
161 G4double f1 = (Emin - adjointPrimKinEnergy
162 G4double f2 = (Emax - adjointPrimKinEnergy
163 gammaE1 = adjointPrimKinEnergy / (1. - f1
164 gammaE2 = gammaE1 - adjointPrimKinEnergy;
165 diffCSUsed =
166 diffCSUsed *
167 (1. + 2. * std::log(1. + electron_mass_c
168 adjointPrimKinEnergy / gammaE1 / gammaE2
169 }
170 else
171 {
172 G4double Emax =
173 GetSecondAdjEnergyMaxForScatProjToProj(a
174 G4double Emin =
175 GetSecondAdjEnergyMinForScatProjToProj(a
176 if(Emin >= Emax)
177 return;
178 gammaE2 = adjointPrimKinEnergy;
179 gammaE1 = Emin * std::pow(Emax / Emin,
180 diffCSUsed = diffCSUsed / gammaE1;
181 }
182
183 // Weight correction
184 // First w_corr is set to the ratio between
185 G4double w_corr = fOutsideWeightFactor;
186 if(fInModelWeightCorr)
187 {
188 w_corr =
189 G4AdjointCSManager::GetAdjointCSManager(
190 }
191 // Then another correction is needed because
192 // been used rather than the one consistent
193
194 G4double diffCS =
195 DiffCrossSectionPerAtomPrimToScatPrim(gamm
196 if(diffCS > 0.)
197 diffCS /= fDirectCS; // here we have the
198 // And we remultiply by the lambda of the fo
199 diffCS *= fDirectProcess->GetCrossSection(ga
200
201 w_corr *= diffCS / diffCSUsed;
202
203 G4double new_weight = aTrack.GetWeight() * w
204 fParticleChange->SetParentWeightByProcess(fa
205 fParticleChange->SetSecondaryWeightByProcess
206 fParticleChange->ProposeParentWeight(new_wei
207
208 G4double cos_th = 1. + electron_mass_c2 * (1
209 if(!isScatProjToProj)
210 {
211 G4double p_elec = theAdjointPrimary->GetTo
212 cos_th = (gammaE1 - gammaE2 * cos
213 }
214 G4double sin_th = 0.;
215 if(std::abs(cos_th) > 1.)
216 {
217 if(cos_th > 0.)
218 {
219 cos_th = 1.;
220 }
221 else
222 cos_th = -1.;
223 }
224 else
225 sin_th = std::sqrt(1. - cos_th * cos_th);
226
227 // gamma0 momentum
228 G4ThreeVector dir_parallel = theAdjointPrima
229 G4double phi = G4UniformRand()
230 G4ThreeVector gammaMomentum1 =
231 gammaE1 *
232 G4ThreeVector(std::cos(phi) * sin_th, std:
233 gammaMomentum1.rotateUz(dir_parallel);
234
235 if(!isScatProjToProj)
236 { // kill the primary and add a secondary
237 fParticleChange->ProposeTrackStatus(fStopA
238 fParticleChange->AddSecondary(
239 new G4DynamicParticle(fAdjEquivDirectPri
240 }
241 else
242 {
243 fParticleChange->ProposeEnergy(gammaE1);
244 fParticleChange->ProposeMomentumDirection(
245 }
246 }
247
248 //////////////////////////////////////////////
249 // The implementation here is correct for ener
250 // photoelectric and compton scattering the me
251 G4double G4AdjointComptonModel::DiffCrossSecti
252 G4double gamEnergy0, G4double kinEnergyElec,
253 {
254 G4double gamEnergy1 = gamEnergy0 - kinEner
255 G4double dSigmadEprod = 0.;
256 if(gamEnergy1 > 0.)
257 dSigmadEprod =
258 DiffCrossSectionPerAtomPrimToScatPrim(ga
259 return dSigmadEprod;
260 }
261
262 //////////////////////////////////////////////
263 G4double G4AdjointComptonModel::DiffCrossSecti
264 G4double gamEnergy0, G4double gamEnergy1, G4
265 {
266 // Based on Klein Nishina formula
267 // In the forward case (see G4KleinNishinaCo
268 // parametrised but the secondaries are samp
269 // differential cross section. The different
270 // is therefore the cross section multiplied
271 // differential Klein Nishina cross section
272
273 // Klein Nishina Cross Section
274 G4double epsilon = gamEnergy0 / el
275 G4double one_plus_two_epsi = 1. + 2. * epsil
276 if(gamEnergy1 > gamEnergy0 || gamEnergy1 < g
277 {
278 return 0.;
279 }
280
281 G4double CS = std::log(one_plus_two_epsi) *
282 (1. - 2. * (1. + epsilon) / (e
283 CS +=
284 4. / epsilon + 0.5 * (1. - 1. / (one_plus_
285 CS /= epsilon;
286 // Note that the pi*re2*Z factor is neglecte
287 // computing dCS_dE1/CS in the differential
288
289 // Klein Nishina Differential Cross Section
290 G4double epsilon1 = gamEnergy1 / electron_ma
291 G4double v = epsilon1 / epsilon;
292 G4double term1 = 1. + 1. / epsilon - 1. /
293 G4double dCS_dE1 = 1. / v + v + term1 * ter
294 dCS_dE1 *= 1. / epsilon / gamEnergy0;
295
296 // Normalised to the CS used in G4
297 fDirectCS = fDirectModel->ComputeCrossSectio
298 G4Gamma::Gamma(), gamEnergy0, Z, 0., 0., 0
299
300 dCS_dE1 *= fDirectCS / CS;
301
302 return dCS_dE1;
303 }
304
305 //////////////////////////////////////////////
306 G4double G4AdjointComptonModel::GetSecondAdjEn
307 G4double primAdjEnergy)
308 {
309 G4double inv_e_max = 1. / primAdjEnergy - 2.
310 G4double e_max = GetHighEnergyLimit();
311 if(inv_e_max > 0.)
312 e_max = std::min(1. / inv_e_max, e_max);
313 return e_max;
314 }
315
316 //////////////////////////////////////////////
317 G4double G4AdjointComptonModel::GetSecondAdjEn
318 G4double primAdjEnergy)
319 {
320 G4double half_e = primAdjEnergy / 2.;
321 return half_e + std::sqrt(half_e * (electron
322 }
323
324 //////////////////////////////////////////////
325 G4double G4AdjointComptonModel::AdjointCrossSe
326 const G4MaterialCutsCouple* aCouple, G4doubl
327 G4bool isScatProjToProj)
328 {
329 if(fUseMatrix)
330 return G4VEmAdjointModel::AdjointCrossSect
331
332 DefineCurrentMaterial(aCouple);
333
334 G4float Cross = 0.;
335 G4float Emax_proj = 0.;
336 G4float Emin_proj = 0.;
337 if(!isScatProjToProj)
338 {
339 Emax_proj = GetSecondAdjEnergyMaxForProdTo
340 Emin_proj = GetSecondAdjEnergyMinForProdTo
341 if(Emax_proj > Emin_proj)
342 {
343 Cross = 0.1 *
344 std::log((Emax_proj - G4float(pr
345 Emax_proj / (Emin_proj
346 (1. + 2. * std::log(G4float(1. +
347 }
348 }
349 else
350 {
351 Emax_proj = GetSecondAdjEnergyMaxForScatPr
352 Emin_proj = GetSecondAdjEnergyMinForScatPr
353 if(Emax_proj > Emin_proj)
354 {
355 Cross = 0.1 * std::log(Emax_proj / Emin_
356 }
357 }
358
359 Cross *= fCurrentMaterial->GetElectronDensit
360 fLastCS = Cross;
361 return double(Cross);
362 }
363