| Geant4 Cross Reference |
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 //
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::AdjointGamma();
50 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
51 fDirectPrimaryPart = G4Gamma::Gamma();
52 fSecondPartSameType = false;
53 fDirectModel =
54 new G4KleinNishinaCompton(G4Gamma::Gamma(), "ComptonDirectModel");
55 }
56
57 ////////////////////////////////////////////////////////////////////////////////
58 G4AdjointComptonModel::~G4AdjointComptonModel() {}
59
60 ////////////////////////////////////////////////////////////////////////////////
61 void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
62 G4bool isScatProjToProj,
63 G4ParticleChange* fParticleChange)
64 {
65 if(!fUseMatrix)
66 return RapidSampleSecondaries(aTrack, isScatProjToProj, fParticleChange);
67
68 // A recall of the compton scattering law:
69 // Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th))
70 // Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1
71 // and Egamma2_min= Egamma2(cos_th=-1) =
72 // Egamma1/(1+2.(Egamma1/E0_electron))
73 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
74 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
75 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
76 {
77 return;
78 }
79
80 // Sample secondary energy
81 G4double gammaE1;
82 gammaE1 =
83 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
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. / gammaE1 - 1. / gammaE2);
92 if(!isScatProjToProj)
93 {
94 cos_th =
95 (gammaE1 - gammaE2 * cos_th) / theAdjointPrimary->GetTotalMomentum();
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 = theAdjointPrimary->GetMomentumDirection();
113 G4double phi = G4UniformRand() * twopi;
114 G4ThreeVector gammaMomentum1 =
115 gammaE1 *
116 G4ThreeVector(std::cos(phi) * sin_th, std::sin(phi) * sin_th, cos_th);
117 gammaMomentum1.rotateUz(dir_parallel);
118
119 // correct the weight of particles before adding the secondary
120 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
121 adjointPrimKinEnergy, gammaE1, isScatProjToProj);
122
123 if(!isScatProjToProj)
124 { // kill the primary and add a secondary
125 fParticleChange->ProposeTrackStatus(fStopAndKill);
126 fParticleChange->AddSecondary(
127 new G4DynamicParticle(fAdjEquivDirectPrimPart, gammaMomentum1));
128 }
129 else
130 {
131 fParticleChange->ProposeEnergy(gammaE1);
132 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
133 }
134 }
135
136 ////////////////////////////////////////////////////////////////////////////////
137 void G4AdjointComptonModel::RapidSampleSecondaries(
138 const G4Track& aTrack, G4bool isScatProjToProj,
139 G4ParticleChange* fParticleChange)
140 {
141 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
142 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
143
144 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
145
146 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
147 {
148 return;
149 }
150
151 G4double diffCSUsed =
152 0.1 * fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2;
153 G4double gammaE1 = 0.;
154 G4double gammaE2 = 0.;
155 if(!isScatProjToProj)
156 {
157 G4double Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
158 G4double Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
159 if(Emin >= Emax)
160 return;
161 G4double f1 = (Emin - adjointPrimKinEnergy) / Emin;
162 G4double f2 = (Emax - adjointPrimKinEnergy) / Emax / f1;
163 gammaE1 = adjointPrimKinEnergy / (1. - f1 * std::pow(f2, G4UniformRand()));
164 gammaE2 = gammaE1 - adjointPrimKinEnergy;
165 diffCSUsed =
166 diffCSUsed *
167 (1. + 2. * std::log(1. + electron_mass_c2 / adjointPrimKinEnergy)) *
168 adjointPrimKinEnergy / gammaE1 / gammaE2;
169 }
170 else
171 {
172 G4double Emax =
173 GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
174 G4double Emin =
175 GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, fTcutSecond);
176 if(Emin >= Emax)
177 return;
178 gammaE2 = adjointPrimKinEnergy;
179 gammaE1 = Emin * std::pow(Emax / Emin, G4UniformRand());
180 diffCSUsed = diffCSUsed / gammaE1;
181 }
182
183 // Weight correction
184 // First w_corr is set to the ratio between adjoint total CS and fwd total CS
185 G4double w_corr = fOutsideWeightFactor;
186 if(fInModelWeightCorr)
187 {
188 w_corr =
189 G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
190 }
191 // Then another correction is needed because a biased differential CS has
192 // been used rather than the one consistent with the direct model
193
194 G4double diffCS =
195 DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2, 1, 0.);
196 if(diffCS > 0.)
197 diffCS /= fDirectCS; // here we have the normalised diffCS
198 // And we remultiply by the lambda of the forward process
199 diffCS *= fDirectProcess->GetCrossSection(gammaE1, fCurrentCouple);
200
201 w_corr *= diffCS / diffCSUsed;
202
203 G4double new_weight = aTrack.GetWeight() * w_corr;
204 fParticleChange->SetParentWeightByProcess(false);
205 fParticleChange->SetSecondaryWeightByProcess(false);
206 fParticleChange->ProposeParentWeight(new_weight);
207
208 G4double cos_th = 1. + electron_mass_c2 * (1. / gammaE1 - 1. / gammaE2);
209 if(!isScatProjToProj)
210 {
211 G4double p_elec = theAdjointPrimary->GetTotalMomentum();
212 cos_th = (gammaE1 - gammaE2 * cos_th) / p_elec;
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 = theAdjointPrimary->GetMomentumDirection();
229 G4double phi = G4UniformRand() * twopi;
230 G4ThreeVector gammaMomentum1 =
231 gammaE1 *
232 G4ThreeVector(std::cos(phi) * sin_th, std::sin(phi) * sin_th, cos_th);
233 gammaMomentum1.rotateUz(dir_parallel);
234
235 if(!isScatProjToProj)
236 { // kill the primary and add a secondary
237 fParticleChange->ProposeTrackStatus(fStopAndKill);
238 fParticleChange->AddSecondary(
239 new G4DynamicParticle(fAdjEquivDirectPrimPart, gammaMomentum1));
240 }
241 else
242 {
243 fParticleChange->ProposeEnergy(gammaE1);
244 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
245 }
246 }
247
248 ////////////////////////////////////////////////////////////////////////////////
249 // The implementation here is correct for energy loss process, for the
250 // photoelectric and compton scattering the method should be redefined
251 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
252 G4double gamEnergy0, G4double kinEnergyElec, G4double Z, G4double A)
253 {
254 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
255 G4double dSigmadEprod = 0.;
256 if(gamEnergy1 > 0.)
257 dSigmadEprod =
258 DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0, gamEnergy1, Z, A);
259 return dSigmadEprod;
260 }
261
262 ////////////////////////////////////////////////////////////////////////////////
263 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
264 G4double gamEnergy0, G4double gamEnergy1, G4double Z, G4double)
265 {
266 // Based on Klein Nishina formula
267 // In the forward case (see G4KleinNishinaCompton) the cross section is
268 // parametrised but the secondaries are sampled from the Klein Nishina
269 // differential cross section. The differential cross section used here
270 // is therefore the cross section multiplied by the normalised
271 // differential Klein Nishina cross section
272
273 // Klein Nishina Cross Section
274 G4double epsilon = gamEnergy0 / electron_mass_c2;
275 G4double one_plus_two_epsi = 1. + 2. * epsilon;
276 if(gamEnergy1 > gamEnergy0 || gamEnergy1 < gamEnergy0 / one_plus_two_epsi)
277 {
278 return 0.;
279 }
280
281 G4double CS = std::log(one_plus_two_epsi) *
282 (1. - 2. * (1. + epsilon) / (epsilon * epsilon));
283 CS +=
284 4. / epsilon + 0.5 * (1. - 1. / (one_plus_two_epsi * one_plus_two_epsi));
285 CS /= epsilon;
286 // Note that the pi*re2*Z factor is neglected because it is suppressed when
287 // computing dCS_dE1/CS in the differential cross section
288
289 // Klein Nishina Differential Cross Section
290 G4double epsilon1 = gamEnergy1 / electron_mass_c2;
291 G4double v = epsilon1 / epsilon;
292 G4double term1 = 1. + 1. / epsilon - 1. / epsilon1;
293 G4double dCS_dE1 = 1. / v + v + term1 * term1 - 1.;
294 dCS_dE1 *= 1. / epsilon / gamEnergy0;
295
296 // Normalised to the CS used in G4
297 fDirectCS = fDirectModel->ComputeCrossSectionPerAtom(
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::GetSecondAdjEnergyMaxForScatProjToProj(
307 G4double primAdjEnergy)
308 {
309 G4double inv_e_max = 1. / primAdjEnergy - 2. / electron_mass_c2;
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::GetSecondAdjEnergyMinForProdToProj(
318 G4double primAdjEnergy)
319 {
320 G4double half_e = primAdjEnergy / 2.;
321 return half_e + std::sqrt(half_e * (electron_mass_c2 + half_e));
322 }
323
324 ////////////////////////////////////////////////////////////////////////////////
325 G4double G4AdjointComptonModel::AdjointCrossSection(
326 const G4MaterialCutsCouple* aCouple, G4double primEnergy,
327 G4bool isScatProjToProj)
328 {
329 if(fUseMatrix)
330 return G4VEmAdjointModel::AdjointCrossSection(aCouple, primEnergy,
331 isScatProjToProj);
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 = GetSecondAdjEnergyMaxForProdToProj(primEnergy);
340 Emin_proj = GetSecondAdjEnergyMinForProdToProj(primEnergy);
341 if(Emax_proj > Emin_proj)
342 {
343 Cross = 0.1 *
344 std::log((Emax_proj - G4float(primEnergy)) * Emin_proj /
345 Emax_proj / (Emin_proj - primEnergy)) *
346 (1. + 2. * std::log(G4float(1. + electron_mass_c2 / primEnergy)));
347 }
348 }
349 else
350 {
351 Emax_proj = GetSecondAdjEnergyMaxForScatProjToProj(primEnergy);
352 Emin_proj = GetSecondAdjEnergyMinForScatProjToProj(primEnergy, 0.);
353 if(Emax_proj > Emin_proj)
354 {
355 Cross = 0.1 * std::log(Emax_proj / Emin_proj);
356 }
357 }
358
359 Cross *= fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2;
360 fLastCS = Cross;
361 return double(Cross);
362 }
363