| 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 #include "G4AdjointIonIonisationModel.hh"
28
29 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointElectron.hh"
31 #include "G4BetheBlochModel.hh"
32 #include "G4BraggIonModel.hh"
33 #include "G4GenericIon.hh"
34 #include "G4NistManager.hh"
35 #include "G4ParticleChange.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4TrackStatus.hh"
39
40 //////////////////////////////////////////////
41 G4AdjointIonIonisationModel::G4AdjointIonIonis
42 : G4VEmAdjointModel("Adjoint_IonIonisation")
43 {
44 fUseMatrix = true;
45 fUseMatrixPerElement = true;
46 fApplyCutInRange = true;
47 fOneMatrixForAllElements = true;
48 fSecondPartSameType = false;
49
50 // The direct EM Model is taken as BetheBloc
51 // computation of the differential cross sec
52 // The Bragg model could be used as an alter
53 // differential cross section
54
55 fBetheBlochDirectEMModel = new G4BetheBloch
56 fBraggIonDirectEMModel = new G4BraggIonMo
57 fAdjEquivDirectSecondPart = G4AdjointElectro
58 fDirectPrimaryPart = nullptr;
59 }
60
61 //////////////////////////////////////////////
62 G4AdjointIonIonisationModel::~G4AdjointIonIoni
63
64 //////////////////////////////////////////////
65 void G4AdjointIonIonisationModel::SampleSecond
66 const G4Track& aTrack, G4bool isScatProjToPr
67 G4ParticleChange* fParticleChange)
68 {
69 const G4DynamicParticle* theAdjointPrimary =
70
71 // Elastic inverse scattering
72 G4double adjointPrimKinEnergy = theAdjointPr
73 G4double adjointPrimP = theAdjointPr
74
75 if(adjointPrimKinEnergy > GetHighEnergyLimit
76 {
77 return;
78 }
79
80 // Sample secondary energy
81 G4double projectileKinEnergy =
82 SampleAdjSecEnergyFromCSMatrix(adjointPrim
83 // Caution !!!this weight correction should
84 CorrectPostStepWeight(fParticleChange, aTrac
85 adjointPrimKinEnergy,
86 isScatProjToProj);
87
88 // Kinematics:
89 // we consider a two body elastic scattering
90 // the projectile knock on an e- at rest and
91 G4double projectileM0 = fAdjEquivDi
92 G4double projectileTotalEnergy = projectileM
93 G4double projectileP2 =
94 projectileTotalEnergy * projectileTotalEne
95
96 // Companion
97 G4double companionM0 = fAdjEquivDirectPrimPa
98 if(isScatProjToProj)
99 {
100 companionM0 = fAdjEquivDirectSecondPart->G
101 }
102 G4double companionTotalEnergy =
103 companionM0 + projectileKinEnergy - adjoin
104 G4double companionP2 =
105 companionTotalEnergy * companionTotalEnerg
106
107 // Projectile momentum
108 G4double P_parallel =
109 (adjointPrimP * adjointPrimP + projectileP
110 (2. * adjointPrimP);
111 G4double P_perp = std::sqrt(projectileP2 - P
112 G4ThreeVector dir_parallel = theAdjointPrima
113 G4double phi = G4UniformRand()
114 G4ThreeVector projectileMomentum =
115 G4ThreeVector(P_perp * std::cos(phi), P_pe
116 projectileMomentum.rotateUz(dir_parallel);
117
118 if(!isScatProjToProj)
119 { // kill the primary and add a secondary
120 fParticleChange->ProposeTrackStatus(fStopA
121 fParticleChange->AddSecondary(
122 new G4DynamicParticle(fAdjEquivDirectPri
123 }
124 else
125 {
126 fParticleChange->ProposeEnergy(projectileK
127 fParticleChange->ProposeMomentumDirection(
128 }
129 }
130
131 //////////////////////////////////////////////
132 G4double G4AdjointIonIonisationModel::DiffCros
133 G4double kinEnergyProj, G4double kinEnergyPr
134 {
135 // Probably that here the Bragg Model should
136 // kinEnergyProj/nuc<2MeV
137 G4double dSigmadEprod = 0.;
138 G4double Emax_proj = GetSecondAdjEnergyMa
139 G4double Emin_proj = GetSecondAdjEnergyMi
140
141 G4double kinEnergyProjScaled = fMassRatio *
142
143 // the produced particle should have a kinet
144 // projectile
145 if(kinEnergyProj > Emin_proj && kinEnergyPro
146 {
147 G4double Tmax = kinEnergyProj;
148
149 G4double E1 = kinEnergyProd;
150 G4double E2 = kinEnergyProd * 1.000001;
151 G4double dE = (E2 - E1);
152 G4double sigma1, sigma2;
153 fDirectModel = fBraggIonDirectEMModel;
154 if(kinEnergyProjScaled > 2. * MeV && !fUse
155 fDirectModel = fBetheBlochDirectEMModel;
156 sigma1 = fDirectModel->ComputeCrossSection
157 fDirectPrimaryPart, kinEnergyProj, Z, A,
158 sigma2 = fDirectModel->ComputeCrossSection
159 fDirectPrimaryPart, kinEnergyProj, Z, A,
160
161 dSigmadEprod = (sigma1 - sigma2) / dE;
162
163 if(dSigmadEprod > 1.)
164 {
165 G4cout << "sigma1 " << kinEnergyProj / M
166 << '\t' << sigma1 << G4endl;
167 G4cout << "sigma2 " << kinEnergyProj / M
168 << '\t' << sigma2 << G4endl;
169 G4cout << "dsigma " << kinEnergyProj / M
170 << '\t' << dSigmadEprod << G4endl
171 }
172
173 if(fDirectModel == fBetheBlochDirectEMMode
174 {
175 // correction of differential cross sect
176 // the suppression of particle at second
177 // Bethe Bloch Model. This correction co
178 // probability function used to test the
179 // code taken from G4BetheBlochModel::Sa
180
181 G4double deltaKinEnergy = kinEnergyProd;
182
183 G4double x = fFormFact * deltaKinEnergy;
184 if(x > 1.e-6)
185 {
186 G4double totEnergy = kinEnergyProj + f
187 G4double etot2 = totEnergy * totEn
188 G4double beta2 = kinEnergyProj * (kinE
189 G4double f1 = 0.0;
190 G4double f = 1.0 - beta2 * deltaKi
191 if(0.5 == fSpin)
192 {
193 f1 = 0.5 * deltaKinEnergy * deltaKin
194 f += f1;
195 }
196 G4double x1 = 1.0 + x;
197 G4double gg = 1.0 / (x1 * x1);
198 if(0.5 == fSpin)
199 {
200 G4double x2 =
201 0.5 * electron_mass_c2 * deltaKinE
202 gg *= (1.0 + fMagMoment2 * (x2 - f1
203 }
204 if(gg > 1.0)
205 {
206 G4cout << "### G4BetheBlochModel in
207 << G4endl;
208 gg = 1.;
209 }
210 dSigmadEprod *= gg;
211 }
212 }
213 }
214 return dSigmadEprod;
215 }
216
217 //////////////////////////////////////////////
218 void G4AdjointIonIonisationModel::SetIon(G4Par
219 G4Par
220 {
221 fDirectPrimaryPart = fwd_ion;
222 fAdjEquivDirectPrimPart = adj_ion;
223
224 DefineProjectileProperty();
225 }
226
227 //////////////////////////////////////////////
228 void G4AdjointIonIonisationModel::CorrectPostS
229 G4ParticleChange* fParticleChange, G4double
230 G4double adjointPrimKinEnergy, G4double proj
231 {
232 // It is needed because the direct cross sec
233 // differential cross section is not the one
234 // the direct model where the GenericIon stu
235 // of effective charge. In the direct model
236 // not reflect the integral cross section. T
237 // we used to compute the differential CS ma
238 // the forward case despite the fact that it
239 // the FWD case. For this reason an extra we
240 // end.
241
242 G4double new_weight = old_weight;
243
244 // the correction of CS due to the problem e
245 G4double kinEnergyProjScaled = fMassRatio *
246 fDirectModel = fBraggIonDire
247 if(kinEnergyProjScaled > 2. * MeV && !fUseOn
248 fDirectModel = fBetheBlochDirectEMModel;
249 G4double UsedFwdCS = fDirectModel->ComputeCr
250 fDirectPrimaryPart, projectileKinEnergy, 1
251 G4double chargeSqRatio = 1.;
252 if(fChargeSquare > 1.)
253 chargeSqRatio = fDirectModel->GetChargeSqu
254 fDirectPrimaryPart, fCurrentMaterial, pr
255 G4double CorrectFwdCS =
256 chargeSqRatio * fDirectModel->ComputeCross
257 G4GenericIon::GenericIon
258 fTcutSecond, 1.e20);
259 // May be some check is needed if UsedFwdCS
260 // avoid a secondary to be produced,
261 if(UsedFwdCS > 0.)
262 new_weight *= CorrectFwdCS / UsedFwdCS;
263
264 // additional CS correction needed for cross
265 // May be wrong for ions. Most of the time n
266 new_weight *=
267 G4AdjointCSManager::GetAdjointCSManager()-
268 fCsBiasingFactor;
269
270 new_weight *= projectileKinEnergy / adjointP
271
272 fParticleChange->SetParentWeightByProcess(fa
273 fParticleChange->SetSecondaryWeightByProcess
274 fParticleChange->ProposeParentWeight(new_wei
275 }
276
277 //////////////////////////////////////////////
278 void G4AdjointIonIonisationModel::DefineProjec
279 {
280 // Slightly modified code taken from G4Bethe
281 G4String pname = fDirectPrimaryPart->GetPart
282
283 fMass = fDirectPrimaryPart->GetPDG
284 fMassRatio = G4GenericIon::GenericIon()
285 fSpin = fDirectPrimaryPart->GetPDG
286 G4double q = fDirectPrimaryPart->GetPDG
287 fChargeSquare = q * q;
288 fRatio = electron_mass_c2 / fMass;
289 fOnePlusRatio2 = (1. + fRatio) * (1. + fRat
290 fOneMinusRatio2 = (1. - fRatio) * (1. - fRat
291 G4double magmom = fDirectPrimaryPart->GetPDG
292 (0.5 * eplus * hbar_Planck
293 fMagMoment2 = magmom * magmom - 1.0;
294 if(fDirectPrimaryPart->GetLeptonNumber() ==
295 {
296 G4double x = 0.8426 * GeV;
297 if(fSpin == 0.0 && fMass < GeV)
298 {
299 x = 0.736 * GeV;
300 }
301 else if(fMass > GeV)
302 {
303 x /= G4NistManager::Instance()->GetZ13(f
304 }
305 fFormFact = 2.0 * electron_mass_c2 / (x *
306 }
307 }
308
309 //////////////////////////////////////////////
310 G4double G4AdjointIonIonisationModel::GetSecon
311 G4double primAdjEnergy)
312 {
313 return primAdjEnergy * fOnePlusRatio2 /
314 (fOneMinusRatio2 - 2. * fRatio * prim
315 }
316
317 //////////////////////////////////////////////
318 G4double G4AdjointIonIonisationModel::GetSecon
319 G4double primAdjEnergy, G4double tcut)
320 {
321 return primAdjEnergy + tcut;
322 }
323
324 //////////////////////////////////////////////
325 G4double G4AdjointIonIonisationModel::GetSecon
326 G4double)
327 {
328 return GetHighEnergyLimit();
329 }
330
331 //////////////////////////////////////////////
332 G4double G4AdjointIonIonisationModel::GetSecon
333 G4double primAdjEnergy)
334 {
335 return (2. * primAdjEnergy - 4. * fMass +
336 std::sqrt(4. * primAdjEnergy * primA
337 8. * primAdjEnergy * fMass
338 4.;
339 }
340