| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 27
28 #include "G4ContinuousGainOfEnergy.hh" 28 #include "G4ContinuousGainOfEnergy.hh"
29 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" 30 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh" <<
38 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh" << 32 #include "G4Step.hh"
>> 33 #include "G4ParticleDefinition.hh"
40 #include "G4VEmModel.hh" 34 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh" << 35 #include "G4VEmFluctuationModel.hh"
42 #include "G4VParticleChange.hh" 36 #include "G4VParticleChange.hh"
>> 37 #include "G4AdjointCSManager.hh"
>> 38 #include "G4LossTableManager.hh"
>> 39 #include "G4SystemOfUnits.hh"
>> 40 #include "G4PhysicalConstants.hh"
43 41
44 ////////////////////////////////////////////// 42 ///////////////////////////////////////////////////////
45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEn << 43 //
46 << 44 G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name,
47 : G4VContinuousProcess(name, type) << 45 G4ProcessType type): G4VContinuousProcess(name, type)
48 {} << 46 {
49 47
50 ////////////////////////////////////////////// << 48
51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfE << 49 linLossLimit=0.05;
>> 50 lossFluctuationArePossible =true;
>> 51 lossFluctuationFlag=true;
>> 52 is_integral = false;
>> 53
>> 54 //Will be properly set in SetDirectParticle()
>> 55 IsIon=false;
>> 56 massRatio =1.;
>> 57 chargeSqRatio=1.;
>> 58 preStepChargeSqRatio=1.;
>> 59
>> 60 //Some initialization
>> 61 currentCoupleIndex=9999999;
>> 62 currentCutInRange=0.;
>> 63 currentMaterialIndex=9999999;
>> 64 currentTcut=0.;
>> 65 preStepKinEnergy=0.;
>> 66 preStepRange=0.;
>> 67 preStepScaledKinEnergy=0.;
>> 68
>> 69 currentCouple=0;
>> 70 }
52 71
53 ////////////////////////////////////////////// 72 ///////////////////////////////////////////////////////
54 void G4ContinuousGainOfEnergy::ProcessDescript << 73 //
>> 74 G4ContinuousGainOfEnergy::~G4ContinuousGainOfEnergy()
55 { 75 {
56 out << "Continuous process acting on adjoint << 76
57 "continuous gain of energy of charged << 77 }
58 "tracked back.\n"; << 78 ///////////////////////////////////////////////////////
>> 79 //
>> 80
>> 81 void G4ContinuousGainOfEnergy::PreparePhysicsTable(
>> 82 const G4ParticleDefinition& )
>> 83 {//theDirectEnergyLossProcess->PreparePhysicsTable(part);
>> 84
>> 85 ;
59 } 86 }
60 87
61 ////////////////////////////////////////////// 88 ///////////////////////////////////////////////////////
62 void G4ContinuousGainOfEnergy::SetDirectPartic << 89 //
63 { << 90
64 fDirectPartDef = p; << 91 void G4ContinuousGainOfEnergy::BuildPhysicsTable(const G4ParticleDefinition&)
65 if(fDirectPartDef->GetParticleType() == "nuc << 92 {//theDirectEnergyLossProcess->BuildPhysicsTable(part);
66 { << 93 ;
67 fIsIon = true; << 94 }
68 fMassRatio = proton_mass_c2 / fDirectPartD << 95
69 } << 96 ///////////////////////////////////////////////////////
>> 97 //
>> 98 void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p)
>> 99 {theDirectPartDef=p;
>> 100 if (theDirectPartDef->GetParticleType()== "nucleus") {
>> 101 IsIon=true;
>> 102 massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass();
>> 103 G4double q=theDirectPartDef->GetPDGCharge();
>> 104 chargeSqRatio=q*q;
>> 105
>> 106
>> 107 }
>> 108
70 } 109 }
71 110
72 ////////////////////////////////////////////// 111 ///////////////////////////////////////////////////////
>> 112 //
>> 113 //
73 G4VParticleChange* G4ContinuousGainOfEnergy::A 114 G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track,
74 << 115 const G4Step& step)
75 { 116 {
76 // Caution in this method the step length sh << 117
77 // A problem is that this is computed by the << 118 //Caution in this method the step length should be the true step length
78 // not know the energy at the end of the adj << 119 // A problem is that this is compute by the multiple scattering that does not know the energy at the end of the adjoint step. This energy is used during the
79 // during the forward sim. Nothing we can re << 120 //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method
80 // time. This is inherent to the MS method << 121 //
81 << 122
>> 123
>> 124
82 aParticleChange.Initialize(track); 125 aParticleChange.Initialize(track);
83 << 126
84 // Get the actual (true) Step length 127 // Get the actual (true) Step length
>> 128 //----------------------------------
85 G4double length = step.GetStepLength(); 129 G4double length = step.GetStepLength();
86 G4double degain = 0.0; << 130 G4double degain = 0.0;
87 << 131
>> 132
>> 133
88 // Compute this for weight change after cont 134 // Compute this for weight change after continuous energy loss
89 G4double DEDX_before = << 135 //-------------------------------------------------------------
90 fDirectEnergyLossProcess->GetDEDX(fPreStep << 136 G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple);
91 << 137
92 // For the fluctuation we generate a new dyn << 138
93 // = preEnergy+egain and then compute the fl << 139
94 // case. << 140 // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain
>> 141 // and then compute the fluctuation given in the direct case.
>> 142 //-----------------------------------------------------------------------
95 G4DynamicParticle* dynParticle = new G4Dynam 143 G4DynamicParticle* dynParticle = new G4DynamicParticle();
96 *dynParticle = *(track.Get << 144 *dynParticle = *(track.GetDynamicParticle());
97 dynParticle->SetDefinition(fDirectPartDef); << 145 dynParticle->SetDefinition(theDirectPartDef);
98 G4double Tkin = dynParticle->GetKineticEnerg << 146 G4double Tkin = dynParticle->GetKineticEnergy();
99 << 147
100 G4double dlength = length; << 148
101 if(Tkin != fPreStepKinEnergy && fIsIon) << 149 size_t n=1;
102 { << 150 if (is_integral ) n=10;
103 G4double chargeSqRatio = fCurrentModel->Ge << 151 n=1;
104 fDirectPartDef, fCurrentMaterial, Tkin); << 152 G4double dlength= length/n;
105 fDirectEnergyLossProcess->SetDynamicMassCh << 153 for (size_t i=0;i<n;i++) {
106 } << 154 if (Tkin != preStepKinEnergy && IsIon) {
107 << 155 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
108 G4double r = fDirectEnergyLossProcess->GetRa << 156 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
109 if(dlength <= fLinLossLimit * r) << 157
110 { << 158 }
111 degain = DEDX_before * dlength; << 159
112 } << 160 G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple);
113 else << 161 if( dlength <= linLossLimit * r ) {
114 { << 162 degain = DEDX_before*dlength;
115 G4double x = r + dlength; << 163 }
116 G4double E = fDirectEnergyLossProcess->Get << 164 else {
117 if(fIsIon) << 165 G4double x = r + dlength;
118 { << 166 //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple);
119 G4double chargeSqRatio = fCurrentModel-> << 167 G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
120 fDirectPartDef, fCurrentMaterial, E); << 168 if (IsIon){
121 fDirectEnergyLossProcess->SetDynamicMass << 169 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
122 G4double x1 = fDirectEnergyLossProcess-> << 170 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
123 << 171 G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
124 G4int ii = 0; << 172
125 constexpr G4int iimax = 100; << 173 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
126 while(std::abs(x - x1) > 0.01 * x) << 174 G4int ii=0;
127 { << 175 const G4int iimax = 100;
128 E = fDirectEnergyLossProcess->GetKinet << 176 while (std::abs(x-x1)>0.01*x) {
129 chargeSqRatio = fCurrentModel->GetChar << 177 E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
130 fDirectPartDef, fCurrentMaterial, E) << 178 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
131 fDirectEnergyLossProcess->SetDynamicMa << 179 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
132 << 180 x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
133 x1 = fDirectEnergyLossProcess->GetRang << 181 ++ii;
134 ++ii; << 182 if(ii >= iimax) { break; }
135 if(ii >= iimax) << 183 }
136 { << 184 }
137 break; << 185
138 } << 186 degain=E-Tkin;
139 } << 187
140 } << 188
141 << 189
142 degain = E - Tkin; << 190 }
143 } << 191 //G4cout<<degain<<G4endl;
144 G4double tmax = fCurrentModel->MaxSecondaryK << 192 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
145 fCurrentTcut = std::min(fCurrentTcut, tmax); << 193 tmax = std::min(tmax,currentTcut);
146 << 194
147 dynParticle->SetKineticEnergy(Tkin + degain) << 195
148 << 196 dynParticle->SetKineticEnergy(Tkin+degain);
149 // Corrections, which cannot be tabulated fo << 197
150 fCurrentModel->CorrectionsAlongStep(fCurrent << 198 // Corrections, which cannot be tabulated for ions
151 << 199 //----------------------------------------
152 // Sample fluctuations << 200 G4double esecdep=0;//not used in most models
153 G4double deltaE = 0.; << 201 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
154 if(fLossFluctuationFlag) << 202
155 { << 203 // Sample fluctuations
156 deltaE = fCurrentModel->GetModelOfFluctuat << 204 //-------------------
157 fCurrentCouple, dynParticle, fCurrentTcu << 205
158 - degain; << 206
159 } << 207 G4double deltaE =0.;
160 << 208 if (lossFluctuationFlag ) {
161 G4double egain = degain + deltaE; << 209 deltaE = currentModel->GetModelOfFluctuations()->
162 if(egain <= 0.) << 210 SampleFluctuations(currentCouple,dynParticle,tmax,dlength,degain)-degain;
163 egain = degain; << 211 }
164 Tkin += egain; << 212
165 dynParticle->SetKineticEnergy(Tkin); << 213 G4double egain=degain+deltaE;
>> 214 if (egain <=0) egain=degain;
>> 215 Tkin+=egain;
>> 216 dynParticle->SetKineticEnergy(Tkin);
>> 217 }
>> 218
>> 219
>> 220
166 221
>> 222
167 delete dynParticle; 223 delete dynParticle;
168 << 224
169 if(fIsIon) << 225 if (IsIon){
170 { << 226 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
171 G4double chargeSqRatio = fCurrentModel->Ge << 227 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
172 fDirectPartDef, fCurrentMaterial, Tkin); << 228
173 fDirectEnergyLossProcess->SetDynamicMassCh <<
174 } 229 }
175 << 230
176 G4double DEDX_after = fDirectEnergyLossProce << 231 G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
177 G4double weight_correction = DEDX_after / DE << 232
178 << 233
>> 234 G4double weight_correction=DEDX_after/DEDX_before;
>> 235
>> 236
179 aParticleChange.ProposeEnergy(Tkin); 237 aParticleChange.ProposeEnergy(Tkin);
180 238
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 239
185 G4double new_weight = << 240 //Caution!!!
186 weight_correction * step.GetPostStepPoint( << 241 // It is important to select the weight of the post_step_point
>> 242 // as the current weight and not the weight of the track, as t
>> 243 // the weight of the track is changed after having applied all
>> 244 // the along_step_do_it.
>> 245
>> 246 // G4double new_weight=weight_correction*track.GetWeight(); //old
>> 247 G4double new_weight=weight_correction*step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(fal 248 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weig 249 aParticleChange.ProposeParentWeight(new_weight);
189 250
>> 251
190 return &aParticleChange; 252 return &aParticleChange;
191 } <<
192 253
>> 254 }
193 ////////////////////////////////////////////// 255 ///////////////////////////////////////////////////////
>> 256 //
194 void G4ContinuousGainOfEnergy::SetLossFluctuat 257 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 { 258 {
196 if(val && !fLossFluctuationArePossible) << 259 if(val && !lossFluctuationArePossible) return;
197 return; << 260 lossFluctuationFlag = val;
198 fLossFluctuationFlag = val; <<
199 } 261 }
200 <<
201 ////////////////////////////////////////////// 262 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo << 263 //
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 264
240 G4double r1 = fDirectEnergyLossProcess->GetR <<
241 265
242 if(fIsIon) <<
243 fDirectEnergyLossProcess->SetDynamicMassCh <<
244 <<
245 266
246 return std::max(r1 - preStepRange, 0.001 * m << 267 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
>> 268 G4double , G4double , G4double& )
>> 269 {
>> 270 G4double x = DBL_MAX;
>> 271 x=.1*mm;
>> 272
>> 273
>> 274 DefineMaterial(track.GetMaterialCutsCouple());
>> 275
>> 276 preStepKinEnergy = track.GetKineticEnergy();
>> 277 preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
>> 278 currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
>> 279 G4double emax_model=currentModel->HighEnergyLimit();
>> 280 if (IsIon) {
>> 281 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
>> 282 preStepChargeSqRatio = chargeSqRatio;
>> 283 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 284 }
>> 285
>> 286
>> 287 G4double maxE =1.1*preStepKinEnergy;
>> 288 /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
>> 289 else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
>> 290 else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
>> 291
>> 292 if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
>> 293
>> 294 maxE=std::min(emax_model*1.001,maxE);
>> 295
>> 296 preStepRange = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
>> 297
>> 298 if (IsIon) {
>> 299 G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
>> 300 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
>> 301 }
>> 302
>> 303 G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
>> 304
>> 305 if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 306
>> 307
>> 308
>> 309 x=r1-preStepRange;
>> 310 x=std::max(r1-preStepRange,0.001*mm);
>> 311
>> 312 return x;
>> 313
>> 314
247 } 315 }
248 << 316 #include "G4EmCorrections.hh"
249 ////////////////////////////////////////////// 317 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC << 318 //
251 << 319
252 { << 320 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
253 G4double ChargeSqRatio = << 321 {
254 G4LossTableManager::Instance()->EmCorrecti << 322
255 fDirectPartDef, fCurrentMaterial, energy << 323 G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
256 if(fDirectEnergyLossProcess) << 324 if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
257 fDirectEnergyLossProcess->SetDynamicMassCh <<
258 } 325 }
259 326