| 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 "G4ContinuousGainOfEnergy.hh"
29
30 #include "G4EmCorrections.hh"
31 #include "G4LossTableManager.hh"
32 #include "G4Material.hh"
33 #include "G4MaterialCutsCouple.hh"
34 #include "G4ParticleChange.hh"
35 #include "G4ParticleDefinition.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh"
40 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh"
42 #include "G4VParticleChange.hh"
43
44 //////////////////////////////////////////////
45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEn
46
47 : G4VContinuousProcess(name, type)
48 {}
49
50 //////////////////////////////////////////////
51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfE
52
53 //////////////////////////////////////////////
54 void G4ContinuousGainOfEnergy::ProcessDescript
55 {
56 out << "Continuous process acting on adjoint
57 "continuous gain of energy of charged
58 "tracked back.\n";
59 }
60
61 //////////////////////////////////////////////
62 void G4ContinuousGainOfEnergy::SetDirectPartic
63 {
64 fDirectPartDef = p;
65 if(fDirectPartDef->GetParticleType() == "nuc
66 {
67 fIsIon = true;
68 fMassRatio = proton_mass_c2 / fDirectPartD
69 }
70 }
71
72 //////////////////////////////////////////////
73 G4VParticleChange* G4ContinuousGainOfEnergy::A
74
75 {
76 // Caution in this method the step length sh
77 // A problem is that this is computed by the
78 // not know the energy at the end of the adj
79 // during the forward sim. Nothing we can re
80 // time. This is inherent to the MS method
81
82 aParticleChange.Initialize(track);
83
84 // Get the actual (true) Step length
85 G4double length = step.GetStepLength();
86 G4double degain = 0.0;
87
88 // Compute this for weight change after cont
89 G4double DEDX_before =
90 fDirectEnergyLossProcess->GetDEDX(fPreStep
91
92 // For the fluctuation we generate a new dyn
93 // = preEnergy+egain and then compute the fl
94 // case.
95 G4DynamicParticle* dynParticle = new G4Dynam
96 *dynParticle = *(track.Get
97 dynParticle->SetDefinition(fDirectPartDef);
98 G4double Tkin = dynParticle->GetKineticEnerg
99
100 G4double dlength = length;
101 if(Tkin != fPreStepKinEnergy && fIsIon)
102 {
103 G4double chargeSqRatio = fCurrentModel->Ge
104 fDirectPartDef, fCurrentMaterial, Tkin);
105 fDirectEnergyLossProcess->SetDynamicMassCh
106 }
107
108 G4double r = fDirectEnergyLossProcess->GetRa
109 if(dlength <= fLinLossLimit * r)
110 {
111 degain = DEDX_before * dlength;
112 }
113 else
114 {
115 G4double x = r + dlength;
116 G4double E = fDirectEnergyLossProcess->Get
117 if(fIsIon)
118 {
119 G4double chargeSqRatio = fCurrentModel->
120 fDirectPartDef, fCurrentMaterial, E);
121 fDirectEnergyLossProcess->SetDynamicMass
122 G4double x1 = fDirectEnergyLossProcess->
123
124 G4int ii = 0;
125 constexpr G4int iimax = 100;
126 while(std::abs(x - x1) > 0.01 * x)
127 {
128 E = fDirectEnergyLossProcess->GetKinet
129 chargeSqRatio = fCurrentModel->GetChar
130 fDirectPartDef, fCurrentMaterial, E)
131 fDirectEnergyLossProcess->SetDynamicMa
132
133 x1 = fDirectEnergyLossProcess->GetRang
134 ++ii;
135 if(ii >= iimax)
136 {
137 break;
138 }
139 }
140 }
141
142 degain = E - Tkin;
143 }
144 G4double tmax = fCurrentModel->MaxSecondaryK
145 fCurrentTcut = std::min(fCurrentTcut, tmax);
146
147 dynParticle->SetKineticEnergy(Tkin + degain)
148
149 // Corrections, which cannot be tabulated fo
150 fCurrentModel->CorrectionsAlongStep(fCurrent
151
152 // Sample fluctuations
153 G4double deltaE = 0.;
154 if(fLossFluctuationFlag)
155 {
156 deltaE = fCurrentModel->GetModelOfFluctuat
157 fCurrentCouple, dynParticle, fCurrentTcu
158 - degain;
159 }
160
161 G4double egain = degain + deltaE;
162 if(egain <= 0.)
163 egain = degain;
164 Tkin += egain;
165 dynParticle->SetKineticEnergy(Tkin);
166
167 delete dynParticle;
168
169 if(fIsIon)
170 {
171 G4double chargeSqRatio = fCurrentModel->Ge
172 fDirectPartDef, fCurrentMaterial, Tkin);
173 fDirectEnergyLossProcess->SetDynamicMassCh
174 }
175
176 G4double DEDX_after = fDirectEnergyLossProce
177 G4double weight_correction = DEDX_after / DE
178
179 aParticleChange.ProposeEnergy(Tkin);
180
181 // Caution!!! It is important to select the
182 // as the current weight and not the weight
183 // the track is changed after having applied
184
185 G4double new_weight =
186 weight_correction * step.GetPostStepPoint(
187 aParticleChange.SetParentWeightByProcess(fal
188 aParticleChange.ProposeParentWeight(new_weig
189
190 return &aParticleChange;
191 }
192
193 //////////////////////////////////////////////
194 void G4ContinuousGainOfEnergy::SetLossFluctuat
195 {
196 if(val && !fLossFluctuationArePossible)
197 return;
198 fLossFluctuationFlag = val;
199 }
200
201 //////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo
203
204
205 {
206 DefineMaterial(track.GetMaterialCutsCouple()
207
208 fPreStepKinEnergy = track.GetKineticEnergy()
209 fCurrentModel = fDirectEnergyLossProcess
210 track.GetKineticEnergy() * fMassRatio, fCu
211 G4double emax_model = fCurrentMode
212 G4double preStepChargeSqRatio = 0.;
213 if(fIsIon)
214 {
215 G4double chargeSqRatio = fCurrentModel->Ge
216 fDirectPartDef, fCurrentMaterial, fPreSt
217 preStepChargeSqRatio = chargeSqRatio;
218 fDirectEnergyLossProcess->SetDynamicMassCh
219
220 }
221
222 G4double maxE = 1.1 * fPreStepKinEnergy;
223
224 if(fPreStepKinEnergy < fCurrentTcut)
225 maxE = std::min(fCurrentTcut, maxE);
226
227 maxE = std::min(emax_model * 1.001, maxE);
228
229 G4double preStepRange =
230 fDirectEnergyLossProcess->GetRange(fPreSte
231
232 if(fIsIon)
233 {
234 G4double chargeSqRatioAtEmax = fCurrentMod
235 fDirectPartDef, fCurrentMaterial, maxE);
236 fDirectEnergyLossProcess->SetDynamicMassCh
237
238 }
239
240 G4double r1 = fDirectEnergyLossProcess->GetR
241
242 if(fIsIon)
243 fDirectEnergyLossProcess->SetDynamicMassCh
244
245
246 return std::max(r1 - preStepRange, 0.001 * m
247 }
248
249 //////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC
251
252 {
253 G4double ChargeSqRatio =
254 G4LossTableManager::Instance()->EmCorrecti
255 fDirectPartDef, fCurrentMaterial, energy
256 if(fDirectEnergyLossProcess)
257 fDirectEnergyLossProcess->SetDynamicMassCh
258 }
259