| 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 // $Id$
26 // 27 //
27 28
28 #include "G4ContinuousGainOfEnergy.hh" 29 #include "G4ContinuousGainOfEnergy.hh"
29 30
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" 31 #include "G4PhysicalConstants.hh"
37 #include "G4Step.hh" <<
38 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
39 #include "G4VEmFluctuationModel.hh" << 33 #include "G4Step.hh"
>> 34 #include "G4ParticleDefinition.hh"
40 #include "G4VEmModel.hh" 35 #include "G4VEmModel.hh"
41 #include "G4VEnergyLossProcess.hh" << 36 #include "G4VEmFluctuationModel.hh"
42 #include "G4VParticleChange.hh" 37 #include "G4VParticleChange.hh"
>> 38 #include "G4UnitsTable.hh"
>> 39 #include "G4AdjointCSManager.hh"
>> 40 #include "G4LossTableManager.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 while (std::abs(x-x1)>0.01*x) {
125 constexpr G4int iimax = 100; << 173 E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple);
126 while(std::abs(x - x1) > 0.01 * x) << 174 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E);
127 { << 175 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
128 E = fDirectEnergyLossProcess->GetKinet << 176 x1= theDirectEnergyLossProcess->GetRange(E, currentCouple);
129 chargeSqRatio = fCurrentModel->GetChar << 177
130 fDirectPartDef, fCurrentMaterial, E) << 178 }
131 fDirectEnergyLossProcess->SetDynamicMa << 179 }
132 << 180
133 x1 = fDirectEnergyLossProcess->GetRang << 181 degain=E-Tkin;
134 ++ii; << 182
135 if(ii >= iimax) << 183
136 { << 184
137 break; << 185 }
138 } << 186 //G4cout<<degain<<G4endl;
139 } << 187 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
140 } << 188 tmax = std::min(tmax,currentTcut);
141 << 189
142 degain = E - Tkin; << 190
143 } << 191 dynParticle->SetKineticEnergy(Tkin+degain);
144 G4double tmax = fCurrentModel->MaxSecondaryK << 192
145 fCurrentTcut = std::min(fCurrentTcut, tmax); << 193 // Corrections, which cannot be tabulated for ions
146 << 194 //----------------------------------------
147 dynParticle->SetKineticEnergy(Tkin + degain) << 195 G4double esecdep=0;//not used in most models
148 << 196 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength);
149 // Corrections, which cannot be tabulated fo << 197
150 fCurrentModel->CorrectionsAlongStep(fCurrent << 198 // Sample fluctuations
151 << 199 //-------------------
152 // Sample fluctuations << 200
153 G4double deltaE = 0.; << 201
154 if(fLossFluctuationFlag) << 202 G4double deltaE =0.;
155 { << 203 if (lossFluctuationFlag ) {
156 deltaE = fCurrentModel->GetModelOfFluctuat << 204 deltaE = currentModel->GetModelOfFluctuations()->
157 fCurrentCouple, dynParticle, fCurrentTcu << 205 SampleFluctuations(currentMaterial,dynParticle,tmax,dlength,degain)-degain;
158 - degain; << 206 }
159 } << 207
160 << 208 G4double egain=degain+deltaE;
161 G4double egain = degain + deltaE; << 209 if (egain <=0) egain=degain;
162 if(egain <= 0.) << 210 Tkin+=egain;
163 egain = degain; << 211 dynParticle->SetKineticEnergy(Tkin);
164 Tkin += egain; << 212 }
165 dynParticle->SetKineticEnergy(Tkin); << 213
>> 214
>> 215
166 216
>> 217
167 delete dynParticle; 218 delete dynParticle;
168 << 219
169 if(fIsIon) << 220 if (IsIon){
170 { << 221 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin);
171 G4double chargeSqRatio = fCurrentModel->Ge << 222 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio);
172 fDirectPartDef, fCurrentMaterial, Tkin); << 223
173 fDirectEnergyLossProcess->SetDynamicMassCh <<
174 } 224 }
175 << 225
176 G4double DEDX_after = fDirectEnergyLossProce << 226 G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple);
177 G4double weight_correction = DEDX_after / DE << 227
178 << 228
>> 229 G4double weight_correction=DEDX_after/DEDX_before;
>> 230
>> 231
179 aParticleChange.ProposeEnergy(Tkin); 232 aParticleChange.ProposeEnergy(Tkin);
180 233
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 234
185 G4double new_weight = << 235 //Caution!!!
186 weight_correction * step.GetPostStepPoint( << 236 // It is important to select the weight of the post_step_point
>> 237 // as the current weight and not the weight of the track, as t
>> 238 // the weight of the track is changed after having applied all
>> 239 // the along_step_do_it.
>> 240
>> 241 // G4double new_weight=weight_correction*track.GetWeight(); //old
>> 242 G4double new_weight=weight_correction*step.GetPostStepPoint()->GetWeight();
187 aParticleChange.SetParentWeightByProcess(fal 243 aParticleChange.SetParentWeightByProcess(false);
188 aParticleChange.ProposeParentWeight(new_weig 244 aParticleChange.ProposeParentWeight(new_weight);
189 245
>> 246
190 return &aParticleChange; 247 return &aParticleChange;
191 } <<
192 248
>> 249 }
193 ////////////////////////////////////////////// 250 ///////////////////////////////////////////////////////
>> 251 //
194 void G4ContinuousGainOfEnergy::SetLossFluctuat 252 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val)
195 { 253 {
196 if(val && !fLossFluctuationArePossible) << 254 if(val && !lossFluctuationArePossible) return;
197 return; << 255 lossFluctuationFlag = val;
198 fLossFluctuationFlag = val; <<
199 } 256 }
200 <<
201 ////////////////////////////////////////////// 257 ///////////////////////////////////////////////////////
202 G4double G4ContinuousGainOfEnergy::GetContinuo << 258 //
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 259
240 G4double r1 = fDirectEnergyLossProcess->GetR <<
241 260
242 if(fIsIon) <<
243 fDirectEnergyLossProcess->SetDynamicMassCh <<
244 <<
245 261
246 return std::max(r1 - preStepRange, 0.001 * m << 262 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track,
>> 263 G4double , G4double , G4double& )
>> 264 {
>> 265 G4double x = DBL_MAX;
>> 266 x=.1*mm;
>> 267
>> 268
>> 269 DefineMaterial(track.GetMaterialCutsCouple());
>> 270
>> 271 preStepKinEnergy = track.GetKineticEnergy();
>> 272 preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio;
>> 273 currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex);
>> 274 G4double emax_model=currentModel->HighEnergyLimit();
>> 275 if (IsIon) {
>> 276 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy);
>> 277 preStepChargeSqRatio = chargeSqRatio;
>> 278 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 279 }
>> 280
>> 281
>> 282 G4double maxE =1.1*preStepKinEnergy;
>> 283 /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy;
>> 284 else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy;
>> 285 else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/
>> 286
>> 287 if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE);
>> 288
>> 289 maxE=std::min(emax_model*1.001,maxE);
>> 290
>> 291 preStepRange = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple);
>> 292
>> 293 if (IsIon) {
>> 294 G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE);
>> 295 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax);
>> 296 }
>> 297
>> 298 G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple);
>> 299
>> 300 if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio);
>> 301
>> 302
>> 303
>> 304 x=r1-preStepRange;
>> 305 x=std::max(r1-preStepRange,0.001*mm);
>> 306
>> 307 return x;
>> 308
>> 309
247 } 310 }
248 << 311 #include "G4EmCorrections.hh"
249 ////////////////////////////////////////////// 312 ///////////////////////////////////////////////////////
250 void G4ContinuousGainOfEnergy::SetDynamicMassC << 313 //
251 << 314
252 { << 315 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy)
253 G4double ChargeSqRatio = << 316 {
254 G4LossTableManager::Instance()->EmCorrecti << 317
255 fDirectPartDef, fCurrentMaterial, energy << 318 G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy);
256 if(fDirectEnergyLossProcess) << 319 if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio);
257 fDirectEnergyLossProcess->SetDynamicMassCh <<
258 } 320 }
259 321