| 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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4DynamicParticleIonisation
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 17.08.2024
37 //
38 // -------------------------------------------
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oo
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 #include "G4DynamicParticleIonisation.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4DynamicParticleFluctuation.hh"
47 #include "G4EmSecondaryParticleType.hh"
48 #include "G4Electron.hh"
49 #include "G4EmParameters.hh"
50 #include "G4EmProcessSubType.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4Material.hh"
55 #include "G4Step.hh"
56 #include "G4Track.hh"
57 #include "G4Log.hh"
58
59 namespace
60 {
61 constexpr G4double ekinLimit = 0.2*CLHEP::Me
62 const G4double twoln10 = 2*G4Log(10.0);
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oo
66
67 G4DynamicParticleIonisation::G4DynamicParticle
68 : G4VContinuousDiscreteProcess("dynPartIoni"
69 {
70 SetVerboseLevel(1);
71 SetProcessSubType(fDynamicIonisation);
72 theElectron = G4Electron::Electron();
73
74 lManager = G4LossTableManager::Instance();
75 lManager->Register(this);
76
77 fUrban = new G4DynamicParticleFluctuation();
78
79 // define these flags only once
80 auto param = G4EmParameters::Instance();
81 fFluct = param->LossFluctuation();
82 fLinLimit = 5*param->LinearLossLimit();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oo
86
87 G4DynamicParticleIonisation::~G4DynamicParticl
88 {
89 lManager->DeRegister(this);
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 void
95 G4DynamicParticleIonisation::BuildPhysicsTable
96 {
97 auto theCoupleTable = G4ProductionCutsTable:
98 fCuts = theCoupleTable->GetEnergyCutsVector(
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void G4DynamicParticleIonisation::PreStepIniti
104 {
105 fCouple = track.GetMaterialCutsCouple();
106 fMaterial = fCouple->GetMaterial();
107 auto dpart = track.GetDynamicParticle();
108 fEkinPreStep = dpart->GetKineticEnergy();
109 fMass = std::max(dpart->GetMass(), CLHEP::el
110 fCharge = dpart->GetCharge()/CLHEP::eplus;
111 fRatio = fMass/CLHEP::proton_mass_c2;
112 fLowestEkin = ekinLimit*fRatio;
113 G4double tau = fEkinPreStep/fMass;
114 G4double ratio = CLHEP::electron_mass_c2/fMa
115 fTmax = 2.0*CLHEP::electron_mass_c2*tau*(tau
116 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
117 fCut = (*fCuts)[fCouple->GetIndex()];
118 fCut = std::max(fCut, fMaterial->GetIonisati
119 fCut = std::min(fCut, fTmax);
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 G4double G4DynamicParticleIonisation::AlongSte
125 const G4Trac
126 G4GPILSelect
127 {
128 *selection = CandidateForSelection;
129
130 // no step limit for the time being
131 return DBL_MAX;
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oo
135
136 G4double G4DynamicParticleIonisation::PostStep
137 const G4Trac
138 G4ForceCondi
139 {
140 *condition = NotForced;
141 G4double x = DBL_MAX;
142 G4double xsec = 0.0;
143
144 PreStepInitialisation(track);
145
146 if (fCharge != 0.0) {
147 xsec = ComputeCrossSection(fEkinPreStep);
148 }
149
150 if (xsec <= 0.0) {
151 theNumberOfInteractionLengthLeft = -1.0;
152 currentInteractionLength = DBL_MAX;
153
154 } else {
155 if (theNumberOfInteractionLengthLeft < 0.0
156
157 theNumberOfInteractionLengthLeft = -G4Lo
158 theInitialNumberOfInteractionLength = th
159
160 } else if(currentInteractionLength < DBL_M
161
162 // subtract NumberOfInteractionLengthLef
163 theNumberOfInteractionLengthLeft -=
164 previousStepSize/currentInteractionLen
165
166 theNumberOfInteractionLengthLeft =
167 std::max(theNumberOfInteractionLengthL
168 }
169 currentInteractionLength = 1.0/xsec;
170 x = theNumberOfInteractionLengthLeft * cur
171 }
172 #ifdef G4VERBOSE
173 if (verboseLevel>2) {
174 G4cout << "G4DynamicParticleIonisation::Po
175 G4cout << " Process: " << GetProcessName(
176 << " for unknown particle Mass(GeV)
177 << " charge=" << fCharge
178 << " Material " << fMaterial->GetN
179 << " Ekin(MeV)=" << fEkinPreStep/C
180 << " MFP(cm)=" << currentInteracti
181 << " ProposedLength(cm)=" << x/CLH
182 }
183 #endif
184 return x;
185 }
186
187 //....oooOO0OOooo........oooOO0OOooo........oo
188
189 G4VParticleChange*
190 G4DynamicParticleIonisation::AlongStepDoIt(con
191 const G4Step& step)
192 {
193 fParticleChange.InitializeForAlongStep(track
194
195 // no energy loss
196 if (fCharge == 0.0) { return &fParticleChang
197
198 // stop low-energy object
199 if (fEkinPreStep <= fLowestEkin) {
200 fParticleChange.SetProposedKineticEnergy(0
201 fParticleChange.ProposeLocalEnergyDeposit(
202 return &fParticleChange;
203 }
204
205 G4double length = step.GetStepLength();
206 G4double dedxPre = ComputeDEDX(fEkinPreStep)
207 G4double eloss = dedxPre*length;
208 G4double ekinPostStep = fEkinPreStep - eloss
209
210 // correction for large step if it is not th
211 if (fEkinPreStep*fLinLimit < eloss && ekinPo
212 G4double dedxPost = ComputeDEDX(ekinPostSt
213 eloss = (eloss + dedxPost*length)*0.5;
214 }
215
216 // do not sample fluctuations at the last st
217 if (fFluct && fEkinPreStep > eloss) {
218 eloss = fUrban->SampleFluctuations(fCouple
219 fCut, fTmax, length, eloss);
220 }
221
222 ekinPostStep = fEkinPreStep - eloss;
223
224 // stop low-energy object
225 if (ekinPostStep <= fLowestEkin) {
226 fParticleChange.SetProposedKineticEnergy(0
227 fParticleChange.ProposeLocalEnergyDeposit(
228 } else {
229 fParticleChange.SetProposedKineticEnergy(e
230 fParticleChange.ProposeLocalEnergyDeposit(
231 }
232 return &fParticleChange;
233 }
234
235 //....oooOO0OOooo........oooOO0OOooo........oo
236
237 G4VParticleChange*
238 G4DynamicParticleIonisation::PostStepDoIt(cons
239 {
240 theNumberOfInteractionLengthLeft = -1.0;
241 fParticleChange.InitializeForPostStep(track)
242
243 auto dp = track.GetDynamicParticle();
244 G4double kinEnergy = dp->GetKineticEnergy();
245 const G4double totEnergy = kinEnergy + fMass
246 const G4double beta2 = kinEnergy*(kinEnergy
247
248 G4double deltaKinEnergy, f;
249
250 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
251 G4double rndm[2];
252
253 // sampling without nuclear size effect
254 do {
255 rndmEngineMod->flatArray(2, rndm);
256 deltaKinEnergy = fCut*fTmax/(fCut*(1.0 - r
257 f = 1.0 - beta2*deltaKinEnergy/fTmax;
258 // Loop checking, 14-Aug-2024, Vladimir Iv
259 } while( rndm[1] > f);
260
261 G4double deltaMomentum =
262 std::sqrt(deltaKinEnergy * (deltaKinEnergy
263 G4double cost = deltaKinEnergy * (totEnergy
264 (deltaMomentum * dp->GetTotalMomentum())
265 cost = std::min(cost, 1.0);
266 const G4double sint = std::sqrt((1.0 - cost)
267 const G4double phi = CLHEP::twopi*rndmEngine
268
269 G4ThreeVector deltaDirection(sint*std::cos(p
270 deltaDirection.rotateUz(dp->GetMomentumDirec
271
272 // create G4DynamicParticle object for delta
273 auto delta = new G4DynamicParticle(theElectr
274 auto t = new G4Track(delta, track.GetGlobalT
275 t->SetTouchableHandle(track.GetTouchableHand
276 t->SetCreatorModelID(fSecID);
277 fParticleChange.AddSecondary(t);
278
279 // Change kinematics of primary particle
280 kinEnergy -= deltaKinEnergy;
281 G4ThreeVector finalP = dp->GetMomentum() - d
282 finalP = finalP.unit();
283
284 fParticleChange.SetProposedKineticEnergy(kin
285 fParticleChange.SetProposedMomentumDirection
286 return &fParticleChange;
287 }
288
289 //....oooOO0OOooo........oooOO0OOooo........oo
290
291 G4double G4DynamicParticleIonisation::ComputeD
292 {
293 G4double tau = ekin/fMass;
294 G4double gam = tau + 1.0;
295 G4double bg2 = tau * (tau + 2.0);
296 G4double beta2 = bg2/(gam*gam);
297 G4double xc = fCut/fTmax;
298
299 G4double exc = fMaterial->GetIonisation()->
300 G4double exc2 = exc*exc;
301
302 // general Bethe-Bloch formula
303 G4double dedx = G4Log(2.0*CLHEP::electron_ma
304
305 // density correction
306 G4double x = G4Log(bg2)/twoln10;
307 dedx -= fMaterial->GetIonisation()->DensityC
308
309 // now compute the total ionization loss per
310 dedx *= CLHEP::twopi_mc2_rcl2*fCharge*fCharg
311 dedx = std::max(dedx, 0.0);
312 return dedx;
313 }
314
315 //....oooOO0OOooo........oooOO0OOooo........oo
316
317 G4double G4DynamicParticleIonisation::ComputeC
318 {
319 G4double cross = 0.0;
320 if (fCut < fTmax) {
321
322 G4double totEnergy = ekin + fMass;
323 G4double energy2 = totEnergy*totEnergy;
324 G4double beta2 = ekin*(ekin + 2.0*fMass)/e
325
326 cross = (fTmax - fCut)/(fCut*fTmax*beta2)
327 cross *= CLHEP::twopi_mc2_rcl2*fCharge*fCh
328 }
329 return cross;
330 }
331
332 //....oooOO0OOooo........oooOO0OOooo........oo
333
334 G4double G4DynamicParticleIonisation::GetMeanF
335
336 {
337 // Note: this method is not used at run-time
338 // It might be eventually refined late
339 *condition = NotForced;
340 return DBL_MAX;
341 }
342
343 //....oooOO0OOooo........oooOO0OOooo........oo
344
345 G4double G4DynamicParticleIonisation::GetConti
346 G4double, G4double&)
347 {
348 return DBL_MAX;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352
353 void G4DynamicParticleIonisation::ProcessDescr
354 {
355 out << "G4DynamicParticleIonisation: dynamic
356 }
357
358 //....oooOO0OOooo........oooOO0OOooo........oo
359