| 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 #include "G4AdjointhIonisationModel.hh"
28
29 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointElectron.hh"
31 #include "G4AdjointProton.hh"
32 #include "G4BetheBlochModel.hh"
33 #include "G4BraggModel.hh"
34 #include "G4NistManager.hh"
35 #include "G4ParticleChange.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4Proton.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4TrackStatus.hh"
40
41 ////////////////////////////////////////////////////////////////////////////////
42 G4AdjointhIonisationModel::G4AdjointhIonisationModel(G4ParticleDefinition* pDef)
43 : G4VEmAdjointModel("Adjoint_hIonisation")
44 {
45 fUseMatrix = true;
46 fUseMatrixPerElement = true;
47 fApplyCutInRange = true;
48 fOneMatrixForAllElements = true;
49 fSecondPartSameType = false;
50
51 // The direct EM Model is taken as BetheBloch. It is only used for the
52 // computation of the differential cross section.
53 // The Bragg model could be used as an alternative as it offers the same
54 // differential cross section
55
56 fDirectModel = new G4BetheBlochModel(pDef);
57 fBraggDirectEMModel = new G4BraggModel(pDef);
58 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
59 fDirectPrimaryPart = pDef;
60
61 if(pDef == G4Proton::Proton())
62 {
63 fAdjEquivDirectPrimPart = G4AdjointProton::AdjointProton();
64 }
65
66 DefineProjectileProperty();
67 }
68
69 ////////////////////////////////////////////////////////////////////////////////
70 G4AdjointhIonisationModel::~G4AdjointhIonisationModel() {}
71
72 ////////////////////////////////////////////////////////////////////////////////
73 void G4AdjointhIonisationModel::SampleSecondaries(
74 const G4Track& aTrack, G4bool isScatProjToProj,
75 G4ParticleChange* fParticleChange)
76 {
77 if(!fUseMatrix)
78 return RapidSampleSecondaries(aTrack, isScatProjToProj, fParticleChange);
79
80 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
81
82 // Elastic inverse scattering
83 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
84 G4double adjointPrimP = theAdjointPrimary->GetTotalMomentum();
85
86 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
87 {
88 return;
89 }
90
91 // Sample secondary energy
92 G4double projectileKinEnergy =
93 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
94 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
95 adjointPrimKinEnergy, projectileKinEnergy,
96 isScatProjToProj);
97 // Caution!!! this weight correction should be always applied
98
99 // Kinematic:
100 // we consider a two body elastic scattering for the forward processes where
101 // the projectile knock on an e- at rest and gives it part of its energy
102 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass();
103 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
104 G4double projectileP2 =
105 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
106
107 // Companion
108 G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
109 if(isScatProjToProj)
110 {
111 companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
112 }
113 G4double companionTotalEnergy =
114 companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
115 G4double companionP2 =
116 companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
117
118 // Projectile momentum
119 G4double P_parallel =
120 (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
121 (2. * adjointPrimP);
122 G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
123 G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
124 G4double phi = G4UniformRand() * twopi;
125 G4ThreeVector projectileMomentum =
126 G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
127 projectileMomentum.rotateUz(dir_parallel);
128
129 if(!isScatProjToProj)
130 { // kill the primary and add a secondary
131 fParticleChange->ProposeTrackStatus(fStopAndKill);
132 fParticleChange->AddSecondary(
133 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
134 }
135 else
136 {
137 fParticleChange->ProposeEnergy(projectileKinEnergy);
138 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
139 }
140 }
141
142 ////////////////////////////////////////////////////////////////////////////////
143 void G4AdjointhIonisationModel::RapidSampleSecondaries(
144 const G4Track& aTrack, G4bool isScatProjToProj,
145 G4ParticleChange* fParticleChange)
146 {
147 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
148 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
149
150 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
151 G4double adjointPrimP = theAdjointPrimary->GetTotalMomentum();
152
153 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
154 {
155 return;
156 }
157
158 G4double projectileKinEnergy = 0.;
159 G4double eEnergy = 0.;
160 G4double newCS =
161 fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2 * fMass;
162 if(!isScatProjToProj)
163 { // 1/E^2 distribution
164
165 eEnergy = adjointPrimKinEnergy;
166 G4double Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
167 G4double Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
168 if(Emin >= Emax)
169 return;
170 G4double a = 1. / Emax;
171 G4double b = 1. / Emin;
172 newCS = newCS * (b - a) / eEnergy;
173
174 projectileKinEnergy = 1. / (b - (b - a) * G4UniformRand());
175 }
176 else
177 {
178 G4double Emax =
179 GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
180 G4double Emin =
181 GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, fTcutSecond);
182 if(Emin >= Emax)
183 return;
184 G4double diff1 = Emin - adjointPrimKinEnergy;
185 G4double diff2 = Emax - adjointPrimKinEnergy;
186
187 G4double t1 = adjointPrimKinEnergy * (1. / diff1 - 1. / diff2);
188 G4double t2 = adjointPrimKinEnergy * (1. / Emin - 1. / Emax);
189 G4double t3 = 2. * std::log(Emax / Emin);
190 G4double sum_t = t1 + t2 + t3;
191 newCS = newCS * sum_t / adjointPrimKinEnergy / adjointPrimKinEnergy;
192 G4double t = G4UniformRand() * sum_t;
193 if(t <= t1)
194 {
195 G4double q = G4UniformRand() * t1 / adjointPrimKinEnergy;
196 projectileKinEnergy = adjointPrimKinEnergy + 1. / (1. / diff1 - q);
197 }
198 else if(t <= t2)
199 {
200 G4double q = G4UniformRand() * t2 / adjointPrimKinEnergy;
201 projectileKinEnergy = 1. / (1. / Emin - q);
202 }
203 else
204 {
205 projectileKinEnergy = Emin * std::pow(Emax / Emin, G4UniformRand());
206 }
207 eEnergy = projectileKinEnergy - adjointPrimKinEnergy;
208 }
209
210 G4double diffCS_perAtom_Used = twopi_mc2_rcl2 * fMass * adjointPrimKinEnergy /
211 projectileKinEnergy / projectileKinEnergy /
212 eEnergy / eEnergy;
213
214 // Weight correction
215 // First w_corr is set to the ratio between adjoint total CS and fwd total CS
216 G4double w_corr =
217 G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
218
219 w_corr *= newCS / fLastCS;
220 // Then another correction is needed due to the fact that a biaised
221 // differential CS has been used rather than the one consistent with the
222 // direct model. Here we consider the true diffCS as the one obtained by the
223 // numerical differentiation over Tcut of the direct CS
224
225 G4double diffCS =
226 DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy, 1, 1);
227 w_corr *= diffCS / diffCS_perAtom_Used;
228
229 if (isScatProjToProj && fTcutSecond>0.005) w_corr=1.;
230
231 G4double new_weight = aTrack.GetWeight() * w_corr;
232 fParticleChange->SetParentWeightByProcess(false);
233 fParticleChange->SetSecondaryWeightByProcess(false);
234 fParticleChange->ProposeParentWeight(new_weight);
235
236 // Kinematic:
237 // we consider a two body elastic scattering for the forward processes where
238 // the projectile knocks on an e- at rest and gives it part of its energy
239 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass();
240 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
241 G4double projectileP2 =
242 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
243
244 // Companion
245 G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
246 if(isScatProjToProj)
247 {
248 companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
249 }
250 G4double companionTotalEnergy =
251 companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
252 G4double companionP2 =
253 companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
254
255 // Projectile momentum
256 G4double P_parallel =
257 (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
258 (2. * adjointPrimP);
259 G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
260 G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
261 G4double phi = G4UniformRand() * twopi;
262 G4ThreeVector projectileMomentum =
263 G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
264 projectileMomentum.rotateUz(dir_parallel);
265
266 if(!isScatProjToProj)
267 { // kill the primary and add a secondary
268 fParticleChange->ProposeTrackStatus(fStopAndKill);
269 fParticleChange->AddSecondary(
270 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
271 }
272 else
273 {
274 fParticleChange->ProposeEnergy(projectileKinEnergy);
275 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
276 }
277 }
278
279 ////////////////////////////////////////////////////////////////////////////////
280 G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
281 G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A)
282 { // Probably here the Bragg Model should be also used for
283 // kinEnergyProj/nuc < 2 MeV
284
285 G4double dSigmadEprod = 0.;
286 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
287 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
288
289 // the produced particle should have a kinetic energy smaller than the
290 // projectile
291 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
292 {
293 G4double Tmax = kinEnergyProj;
294 G4double E1 = kinEnergyProd;
295 //1.0006 factor seems to give the best diff CS, important impact on proton correction factor
296 G4double E2 = kinEnergyProd *1.0006;
297 G4double sigma1, sigma2;
298 if(kinEnergyProj > 2. * MeV)
299 {
300 sigma1 = fDirectModel->ComputeCrossSectionPerAtom(
301 fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
302 sigma2 = fDirectModel->ComputeCrossSectionPerAtom(
303 fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
304 }
305 else
306 {
307 sigma1 = fBraggDirectEMModel->ComputeCrossSectionPerAtom(
308 fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
309 sigma2 = fBraggDirectEMModel->ComputeCrossSectionPerAtom(
310 fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
311 }
312
313 dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
314 if(dSigmadEprod > 1.)
315 {
316 G4cout << "sigma1 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
317 << '\t' << sigma1 << G4endl;
318 G4cout << "sigma2 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
319 << '\t' << sigma2 << G4endl;
320 G4cout << "dsigma " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
321 << '\t' << dSigmadEprod << G4endl;
322 }
323
324 // correction of differential cross section at high energy to correct for
325 // the suppression of particle at secondary at high energy used in the Bethe
326 // Bloch Model. This correction consists of multiplying by g the probability
327 // function used to test the rejection of a secondary. Source code taken
328 // from G4BetheBlochModel::SampleSecondaries
329 G4double deltaKinEnergy = kinEnergyProd;
330
331 // projectile formfactor - suppression of high energy
332 // delta-electron production at high energy
333 G4double x = fFormFact * deltaKinEnergy;
334 if(x > 1.e-6)
335 {
336 G4double totEnergy = kinEnergyProj + fMass;
337 G4double etot2 = totEnergy * totEnergy;
338 G4double beta2 = kinEnergyProj * (kinEnergyProj + 2.0 * fMass) / etot2;
339 G4double f = 1.0 - beta2 * deltaKinEnergy / Tmax;
340 G4double f1 = 0.0;
341 if(0.5 == fSpin)
342 {
343 f1 = 0.5 * deltaKinEnergy * deltaKinEnergy / etot2;
344 f += f1;
345 }
346 G4double x1 = 1.0 + x;
347 G4double gg = 1.0 / (x1 * x1);
348 if(0.5 == fSpin)
349 {
350 G4double x2 = 0.5 * electron_mass_c2 * deltaKinEnergy / (fMass * fMass);
351 gg *= (1.0 + fMagMoment2 * (x2 - f1 / f) / (1.0 + x2));
352 }
353 if(gg > 1.0)
354 {
355 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g
356 << G4endl;
357 gg = 1.;
358 }
359 dSigmadEprod *= gg;
360 }
361 }
362
363 return dSigmadEprod;
364 }
365
366 ////////////////////////////////////////////////////////////////////////////////
367 void G4AdjointhIonisationModel::DefineProjectileProperty()
368 {
369 // Slightly modified code taken from G4BetheBlochModel::SetParticle
370 G4String pname = fDirectPrimaryPart->GetParticleName();
371
372 fMass = fDirectPrimaryPart->GetPDGMass();
373 fSpin = fDirectPrimaryPart->GetPDGSpin();
374 fMassRatio = electron_mass_c2 / fMass;
375 fOnePlusRatio2 = (1. + fMassRatio) * (1. + fMassRatio);
376 fOneMinusRatio2 = (1. - fMassRatio) * (1. - fMassRatio);
377 G4double magmom = fDirectPrimaryPart->GetPDGMagneticMoment() * fMass /
378 (0.5 * eplus * hbar_Planck * c_squared);
379 fMagMoment2 = magmom * magmom - 1.0;
380 fFormFact = 0.0;
381 if(fDirectPrimaryPart->GetLeptonNumber() == 0)
382 {
383 G4double x = 0.8426 * GeV;
384 if(fSpin == 0.0 && fMass < GeV)
385 {
386 x = 0.736 * GeV;
387 }
388 else if(fMass > GeV)
389 {
390 x /= G4NistManager::Instance()->GetZ13(fMass / proton_mass_c2);
391 }
392 fFormFact = 2.0 * electron_mass_c2 / (x * x);
393 }
394 }
395
396 ////////////////////////////////////////////////////////////////////////////////
397 G4double G4AdjointhIonisationModel::AdjointCrossSection(
398 const G4MaterialCutsCouple* aCouple, G4double primEnergy,
399 G4bool isScatProjToProj)
400 {
401 if(fUseMatrix)
402 return G4VEmAdjointModel::AdjointCrossSection(aCouple, primEnergy,
403 isScatProjToProj);
404 DefineCurrentMaterial(aCouple);
405
406 G4double Cross =
407 fCurrentMaterial->GetElectronDensity() * twopi_mc2_rcl2 * fMass;
408
409 if(!isScatProjToProj)
410 {
411 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(primEnergy);
412 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(primEnergy);
413 if(Emax_proj > Emin_proj && primEnergy > fTcutSecond)
414 {
415 Cross *= (1. / Emin_proj - 1. / Emax_proj) / primEnergy;
416 }
417 else
418 Cross = 0.;
419 }
420 else
421 {
422 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProj(primEnergy);
423 G4double Emin_proj =
424 GetSecondAdjEnergyMinForScatProjToProj(primEnergy, fTcutSecond);
425 G4double diff1 = Emin_proj - primEnergy;
426 G4double diff2 = Emax_proj - primEnergy;
427 G4double t1 =
428 (1. / diff1 + 1. / Emin_proj - 1. / diff2 - 1. / Emax_proj) / primEnergy;
429 G4double t2 =
430 2. * std::log(Emax_proj / Emin_proj) / primEnergy / primEnergy;
431 Cross *= (t1 + t2);
432 }
433 fLastCS = Cross;
434 return Cross;
435 }
436
437 //////////////////////////////////////////////////////////////////////////////
438 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProj(
439 G4double primAdjEnergy)
440 {
441 G4double Tmax = primAdjEnergy * fOnePlusRatio2 /
442 (fOneMinusRatio2 - 2. * fMassRatio * primAdjEnergy / fMass);
443 return Tmax;
444 }
445
446 //////////////////////////////////////////////////////////////////////////////
447 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProj(
448 G4double primAdjEnergy, G4double tcut)
449 {
450 return primAdjEnergy + tcut;
451 }
452
453 //////////////////////////////////////////////////////////////////////////////
454 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProj(G4double)
455 {
456 return GetHighEnergyLimit();
457 }
458
459 //////////////////////////////////////////////////////////////////////////////
460 G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProj(
461 G4double primAdjEnergy)
462 {
463 G4double Tmin =
464 (2. * primAdjEnergy - 4. * fMass +
465 std::sqrt(4. * primAdjEnergy * primAdjEnergy + 16. * fMass * fMass +
466 8. * primAdjEnergy * fMass * (1. / fMassRatio + fMassRatio))) /
467 4.;
468 return Tmin;
469 }
470