| 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 // GEANT4 Class header file 28 // GEANT4 Class header file
29 // 29 //
30 // 30 //
31 // File name: G4LindhardSorensenIonModel 31 // File name: G4LindhardSorensenIonModel
32 // 32 //
33 // Author: Alexander Bagulya & Vladimir 33 // Author: Alexander Bagulya & Vladimir Ivanchenko
34 // 34 //
35 // Creation date: 16.04.2018 35 // Creation date: 16.04.2018
36 // 36 //
37 // 37 //
38 // ------------------------------------------- 38 // -------------------------------------------------------------------
39 // 39 //
40 40
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 43
44 #include "G4LindhardSorensenIonModel.hh" 44 #include "G4LindhardSorensenIonModel.hh"
45 #include "Randomize.hh" 45 #include "Randomize.hh"
46 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
48 #include "G4Electron.hh" 48 #include "G4Electron.hh"
49 #include "G4LossTableManager.hh" 49 #include "G4LossTableManager.hh"
50 #include "G4EmCorrections.hh" 50 #include "G4EmCorrections.hh"
51 #include "G4ParticleChangeForLoss.hh" 51 #include "G4ParticleChangeForLoss.hh"
52 #include "G4Log.hh" 52 #include "G4Log.hh"
53 #include "G4DeltaAngle.hh" 53 #include "G4DeltaAngle.hh"
54 #include "G4LindhardSorensenData.hh" 54 #include "G4LindhardSorensenData.hh"
55 #include "G4BraggModel.hh" << 55 #include "G4BraggIonModel.hh"
56 #include "G4BetheBlochModel.hh" 56 #include "G4BetheBlochModel.hh"
57 #include "G4IonICRU73Data.hh" 57 #include "G4IonICRU73Data.hh"
58 #include "G4AutoLock.hh" <<
59 58
60 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 60
62 using namespace std; 61 using namespace std;
63 62
64 G4LindhardSorensenData* G4LindhardSorensenIonM 63 G4LindhardSorensenData* G4LindhardSorensenIonModel::lsdata = nullptr;
65 G4IonICRU73Data* G4LindhardSorensenIonModel::f 64 G4IonICRU73Data* G4LindhardSorensenIonModel::fIonData = nullptr;
66 << 65 std::vector<G4float>* G4LindhardSorensenIonModel::fact[] = {nullptr};
67 namespace <<
68 { <<
69 G4Mutex ionXSMutex = G4MUTEX_INITIALIZER; <<
70 } <<
71 66
72 G4LindhardSorensenIonModel::G4LindhardSorensen 67 G4LindhardSorensenIonModel::G4LindhardSorensenIonModel(const G4ParticleDefinition*,
73 68 const G4String& nam)
74 : G4VEmModel(nam), 69 : G4VEmModel(nam),
>> 70 particle(nullptr),
75 twoln10(2.0*G4Log(10.0)) 71 twoln10(2.0*G4Log(10.0))
76 { 72 {
>> 73 fParticleChange = nullptr;
77 theElectron = G4Electron::Electron(); 74 theElectron = G4Electron::Electron();
78 corr = G4LossTableManager::Instance()->EmCor 75 corr = G4LossTableManager::Instance()->EmCorrections();
79 nist = G4NistManager::Instance(); 76 nist = G4NistManager::Instance();
80 fBraggModel = new G4BraggModel(); << 77 fBraggModel = new G4BraggIonModel();
81 fBBModel = new G4BetheBlochModel(); 78 fBBModel = new G4BetheBlochModel();
82 fElimit = 2.0*CLHEP::MeV; 79 fElimit = 2.0*CLHEP::MeV;
83 } 80 }
84 81
85 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 83
87 G4LindhardSorensenIonModel::~G4LindhardSorense << 84 G4LindhardSorensenIonModel::~G4LindhardSorensenIonModel()
88 if(isFirst) { << 85 {}
89 delete lsdata; <<
90 delete fIonData; <<
91 lsdata = nullptr; <<
92 fIonData = nullptr; <<
93 } <<
94 } <<
95 86
96 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 88
98 void G4LindhardSorensenIonModel::Initialise(co 89 void G4LindhardSorensenIonModel::Initialise(const G4ParticleDefinition* p,
99 co 90 const G4DataVector& ptr)
100 { 91 {
101 fBraggModel->Initialise(p, ptr); 92 fBraggModel->Initialise(p, ptr);
102 fBBModel->Initialise(p, ptr); 93 fBBModel->Initialise(p, ptr);
103 SetParticle(p); 94 SetParticle(p);
>> 95 //G4cout << "G4LindhardSorensenIonModel::Initialise for "
>> 96 // << p->GetParticleName() << G4endl;
104 97
105 // always false before the run 98 // always false before the run
106 SetDeexcitationFlag(false); 99 SetDeexcitationFlag(false);
107 100
108 if(nullptr == fParticleChange) { 101 if(nullptr == fParticleChange) {
109 fParticleChange = GetParticleChangeForLoss 102 fParticleChange = GetParticleChangeForLoss();
110 if(UseAngularGeneratorFlag() && nullptr == 103 if(UseAngularGeneratorFlag() && nullptr == GetAngularDistribution()) {
111 SetAngularDistribution(new G4DeltaAngle( 104 SetAngularDistribution(new G4DeltaAngle());
112 } 105 }
113 } 106 }
114 if(nullptr == lsdata) { << 107 if(IsMaster()) {
115 G4AutoLock l(&ionXSMutex); <<
116 if(nullptr == lsdata) { 108 if(nullptr == lsdata) {
117 isFirst = true; <<
118 lsdata = new G4LindhardSorensenData(); 109 lsdata = new G4LindhardSorensenData();
>> 110 }
>> 111 if(nullptr == fIonData) {
119 fIonData = new G4IonICRU73Data(); 112 fIonData = new G4IonICRU73Data();
120 } 113 }
121 l.unlock(); <<
122 } <<
123 if(isFirst) { <<
124 fIonData->Initialise(); 114 fIonData->Initialise();
125 } 115 }
126 } 116 }
127 117
128 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 119
130 G4double 120 G4double
131 G4LindhardSorensenIonModel::GetChargeSquareRat 121 G4LindhardSorensenIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
132 122 const G4Material* mat,
133 123 G4double kinEnergy)
134 { 124 {
135 chargeSquare = corr->EffectiveChargeSquareRa << 125 return corr->EffectiveChargeSquareRatio(p,mat,kinEnergy);
136 return chargeSquare; <<
137 } 126 }
138 127
139 //....oooOO0OOooo........oooOO0OOooo........oo 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 129
141 G4double 130 G4double
142 G4LindhardSorensenIonModel::GetParticleCharge( 131 G4LindhardSorensenIonModel::GetParticleCharge(const G4ParticleDefinition* p,
143 132 const G4Material* mat,
144 133 G4double kinEnergy)
145 { 134 {
146 return corr->GetParticleCharge(p,mat,kinEner 135 return corr->GetParticleCharge(p,mat,kinEnergy);
147 } 136 }
148 137
149 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 139
151 void G4LindhardSorensenIonModel::SetupParamete 140 void G4LindhardSorensenIonModel::SetupParameters()
152 { 141 {
153 mass = particle->GetPDGMass(); 142 mass = particle->GetPDGMass();
154 spin = particle->GetPDGSpin(); 143 spin = particle->GetPDGSpin();
155 charge = particle->GetPDGCharge()*inveplus; 144 charge = particle->GetPDGCharge()*inveplus;
156 Zin = G4lrint(std::abs(charge)); 145 Zin = G4lrint(std::abs(charge));
157 chargeSquare = charge*charge; 146 chargeSquare = charge*charge;
158 eRatio = CLHEP::electron_mass_c2/mass; 147 eRatio = CLHEP::electron_mass_c2/mass;
159 pRatio = CLHEP::proton_mass_c2/mass; 148 pRatio = CLHEP::proton_mass_c2/mass;
160 const G4double aMag = 149 const G4double aMag =
161 1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CL 150 1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CLHEP::c_squared);
162 G4double magmom = particle->GetPDGMagneticMo 151 G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
163 magMoment2 = magmom*magmom - 1.0; 152 magMoment2 = magmom*magmom - 1.0;
164 G4double x = 0.8426*CLHEP::GeV; 153 G4double x = 0.8426*CLHEP::GeV;
165 if(spin == 0.0 && mass < GeV) { x = 0.736*CL 154 if(spin == 0.0 && mass < GeV) { x = 0.736*CLHEP::GeV; }
166 else if (Zin > 1) { x /= nist->GetA27(Zin); 155 else if (Zin > 1) { x /= nist->GetA27(Zin); }
167 156
168 formfact = 2.0*CLHEP::electron_mass_c2/(x*x) 157 formfact = 2.0*CLHEP::electron_mass_c2/(x*x);
169 tlimit = 2.0/formfact; 158 tlimit = 2.0/formfact;
170 } 159 }
171 160
172 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173 162
174 G4double G4LindhardSorensenIonModel::MinEnergy 163 G4double G4LindhardSorensenIonModel::MinEnergyCut(const G4ParticleDefinition*,
175 const 164 const G4MaterialCutsCouple* couple)
176 { 165 {
177 return couple->GetMaterial()->GetIonisation( 166 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
178 } 167 }
179 168
180 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181 170
182 G4double 171 G4double
183 G4LindhardSorensenIonModel::ComputeCrossSectio 172 G4LindhardSorensenIonModel::ComputeCrossSectionPerElectron(
184 173 const G4ParticleDefinition* p,
185 174 G4double kinEnergy,
186 175 G4double cutEnergy,
187 176 G4double maxKinEnergy)
188 { 177 {
189 // take into account formfactor 178 // take into account formfactor
190 G4double tmax = MaxSecondaryEnergy(p, kinEne 179 G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
191 G4double emax = std::min(tmax, maxKinEnergy) 180 G4double emax = std::min(tmax, maxKinEnergy);
192 G4double escaled = kinEnergy*pRatio; 181 G4double escaled = kinEnergy*pRatio;
193 G4double cross = (escaled <= fElimit) 182 G4double cross = (escaled <= fElimit)
194 ? fBraggModel->ComputeCrossSectionPerElect 183 ? fBraggModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax)
195 : fBBModel->ComputeCrossSectionPerElectron 184 : fBBModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax);
196 // G4cout << "LS: e= " << kinEnergy << " tmi 185 // G4cout << "LS: e= " << kinEnergy << " tmin= " << cutEnergy
197 // << " tmax= " << maxEnergy << " cro 186 // << " tmax= " << maxEnergy << " cross= " << cross << G4endl;
198 return cross; 187 return cross;
199 } 188 }
200 189
201 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202 191
203 G4double G4LindhardSorensenIonModel::ComputeCr 192 G4double G4LindhardSorensenIonModel::ComputeCrossSectionPerAtom(
204 con 193 const G4ParticleDefinition* p,
205 194 G4double kineticEnergy,
206 195 G4double Z, G4double,
207 196 G4double cutEnergy,
208 197 G4double maxEnergy)
209 { 198 {
210 return Z*ComputeCrossSectionPerElectron(p,ki 199 return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
211 } 200 }
212 201
213 //....oooOO0OOooo........oooOO0OOooo........oo 202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
214 203
215 G4double G4LindhardSorensenIonModel::CrossSect 204 G4double G4LindhardSorensenIonModel::CrossSectionPerVolume(
216 con 205 const G4Material* material,
217 con 206 const G4ParticleDefinition* p,
218 207 G4double kineticEnergy,
219 208 G4double cutEnergy,
220 209 G4double maxEnergy)
221 { 210 {
222 return material->GetElectronDensity() 211 return material->GetElectronDensity()
223 *ComputeCrossSectionPerElectron(p,kineticE 212 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
224 } 213 }
225 214
226 //....oooOO0OOooo........oooOO0OOooo........oo 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227 216
228 G4double 217 G4double
229 G4LindhardSorensenIonModel::ComputeDEDXPerVolu 218 G4LindhardSorensenIonModel::ComputeDEDXPerVolume(const G4Material* mat,
230 219 const G4ParticleDefinition* p,
231 220 G4double kinEnergy,
232 221 G4double cut)
233 { 222 {
234 // formfactor is taken into account in Corre 223 // formfactor is taken into account in CorrectionsAlongStep(..)
235 G4double tmax = MaxSecondaryEnergy(p, k 224 G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
236 G4double cutEnergy = std::min(std::min(cut,t 225 G4double cutEnergy = std::min(std::min(cut,tmax), tlimit);
237 226
238 G4double escaled = kinEnergy*pRatio; 227 G4double escaled = kinEnergy*pRatio;
239 G4double dedx = (escaled <= fElimit) 228 G4double dedx = (escaled <= fElimit)
240 ? fBraggModel->ComputeDEDXPerVolume(mat, p 229 ? fBraggModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy)
241 : fBBModel->ComputeDEDXPerVolume(mat, p, k 230 : fBBModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy);
242 231
243 //G4cout << "E(MeV)=" << kinEnergy/MeV << " 232 //G4cout << "E(MeV)=" << kinEnergy/MeV << " dedx=" << dedx
244 // << " " << mat->GetName() << " Ecut(MeV) << 233 // << " " << material->GetName() << " Ecut(MeV)=" << cutEnergy << G4endl;
245 return dedx; 234 return dedx;
246 } 235 }
247 236
248 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 238
250 void G4LindhardSorensenIonModel::CorrectionsAl 239 void G4LindhardSorensenIonModel::CorrectionsAlongStep(
251 const G4Mater 240 const G4MaterialCutsCouple* couple,
252 const G4Dynam 241 const G4DynamicParticle* dp,
253 const G4doubl 242 const G4double& length,
254 G4doubl 243 G4double& eloss)
255 { 244 {
256 // no correction at the last step 245 // no correction at the last step
257 const G4double preKinEnergy = dp->GetKinetic 246 const G4double preKinEnergy = dp->GetKineticEnergy();
258 if(eloss >= preKinEnergy) { return; } 247 if(eloss >= preKinEnergy) { return; }
259 248
260 const G4ParticleDefinition* p = dp->GetDefin 249 const G4ParticleDefinition* p = dp->GetDefinition();
261 SetParticle(p); <<
262 const G4Material* mat = couple->GetMaterial( 250 const G4Material* mat = couple->GetMaterial();
263 const G4double eDensity = mat->GetElectronDe 251 const G4double eDensity = mat->GetElectronDensity();
264 const G4double e = std::max(preKinEnergy - e << 252 const G4double e = std::max(preKinEnergy - eloss*0.5, preKinEnergy*0.75);
265 const G4double tmax = MaxSecondaryEnergy(p, 253 const G4double tmax = MaxSecondaryEnergy(p, e);
266 const G4double escaled = e*pRatio; 254 const G4double escaled = e*pRatio;
267 const G4double tau = e/mass; 255 const G4double tau = e/mass;
268 const G4double q2 = corr->EffectiveChargeSqu <<
269 const G4int Z = p->GetAtomicNumber(); <<
270 256
271 G4double res = 0.0; << 257 G4double q2 = corr->EffectiveChargeSquareRatio(p, mat, e);
>> 258 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
>> 259 const G4int Z = std::max(80, p->GetAtomicNumber());
>> 260
>> 261 G4double res;
272 if(escaled <= fElimit) { 262 if(escaled <= fElimit) {
273 // data from ICRU73 or ICRU90 << 263 res = fIonData->GetDEDX(mat, Z, escaled, G4Log(escaled));
274 if(Z > 2 && Z <= 80) { <<
275 res = fIonData->GetDEDX(mat, Z, escaled, <<
276 /* <<
277 G4cout << " GetDEDX for Z=" << Z << " in " <<
278 << " Escaled=" << escaled << " E=" <<
279 << e << " dEdx=" << res << G4endl; <<
280 */ <<
281 } <<
282 if(res > 0.0) { 264 if(res > 0.0) {
283 auto pcuts = couple->GetProductionCuts() << 265 G4double cut = couple->GetProductionCuts()->GetProductionCut(1);
284 G4double cut = (nullptr == pcuts) ? tmax <<
285 if(cut < tmax) { 266 if(cut < tmax) {
286 const G4double x = cut/tmax; << 267 const G4double x = cut/tmax;
287 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau 268 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x)
288 *q2*CLHEP::twopi_mc2_rcl2*eDensity; 269 *q2*CLHEP::twopi_mc2_rcl2*eDensity;
289 } 270 }
290 res *= length; <<
291 } else { 271 } else {
292 // simplified correction << 272 res = eloss*q2*corr->EffectiveChargeCorrection(p,mat,e)/chargeSquare;
293 res = eloss*q2/chargeSquare; <<
294 } 273 }
295 } else { 274 } else {
296 // Lindhard-Sorensen model <<
297 const G4double gam = tau + 1.0; 275 const G4double gam = tau + 1.0;
298 const G4double beta2 = tau * (tau+2.0)/(ga 276 const G4double beta2 = tau * (tau+2.0)/(gam*gam);
299 G4double deltaL0 = 2.0*corr->BarkasCorrect 277 G4double deltaL0 = 2.0*corr->BarkasCorrection(p, mat, e)*(charge-1.)/charge;
300 G4double deltaL = lsdata->GetDeltaL(Zin, g 278 G4double deltaL = lsdata->GetDeltaL(Zin, gam);
301 279
302 res = eloss + 280 res = eloss +
303 CLHEP::twopi_mc2_rcl2*q2*eDensity*(delta 281 CLHEP::twopi_mc2_rcl2*q2*eDensity*(deltaL+deltaL0)*length/beta2;
304 /* << 282 /*
305 G4cout << " E(GeV)=" << preKinEnergy/GeV << 283 G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)= "
306 << " L= " << eloss*beta2/(twopi_mc2_rcl2* << 284 << preKinEnergy/GeV << " eloss(MeV)= " << eloss
307 << " dL0= " << deltaL0 << 285 << " L= " << eloss*beta2/(twopi_mc2_rcl2*chargeSquare*eDensity*length)
308 << " dL= " << deltaL << " dE(MeV)=" << re << 286 << " dL0= " << deltaL0
309 */ << 287 << " dL= " << deltaL << G4endl;
310 } <<
311 if(res > preKinEnergy || 2*res < eloss) { re <<
312 /* <<
313 G4cout << " G4LindhardSorensenIonModel::Cor <<
314 << preKinEnergy/GeV << " eloss(MeV)=" <<
315 << " res(MeV)=" << res << G4endl; <<
316 */ 288 */
>> 289 }
>> 290 if(res > preKinEnergy) { res = preKinEnergy; }
>> 291 else if(res < 0.0) { res = eloss; }
>> 292 //G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)="
>> 293 // << preKinEnergy/GeV << " eloss(MeV)=" << eloss
>> 294 //<< " res(MeV)=" << res << G4endl;
317 eloss = res; 295 eloss = res;
318 } 296 }
319 297
320 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321 299
322 void G4LindhardSorensenIonModel::SampleSeconda 300 void G4LindhardSorensenIonModel::SampleSecondaries(
323 std::vector<G4DynamicParticle*>* v << 301 vector<G4DynamicParticle*>* vdp,
324 cons 302 const G4MaterialCutsCouple* couple,
325 cons 303 const G4DynamicParticle* dp,
326 G4do << 304 G4double minKinEnergy,
327 G4do 305 G4double maxEnergy)
328 { 306 {
329 G4double kineticEnergy = dp->GetKineticEnerg 307 G4double kineticEnergy = dp->GetKineticEnergy();
330 // take into account formfactor 308 // take into account formfactor
331 const G4double tmax = MaxSecondaryEnergy(dp- << 309 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
332 const G4double minKinEnergy = std::min(cut, << 310
333 const G4double maxKinEnergy = std::min(maxEn << 311 G4double maxKinEnergy = std::min(maxEnergy,tmax);
334 if(minKinEnergy >= maxKinEnergy) { return; } 312 if(minKinEnergy >= maxKinEnergy) { return; }
335 313
336 //G4cout << "G4LindhardSorensenIonModel::Sam 314 //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= "
337 // << minKinEnergy << " Emax= " << maxKinEne 315 // << minKinEnergy << " Emax= " << maxKinEnergy << G4endl;
338 316
339 G4double totEnergy = kineticEnergy + mass; << 317 G4double totEnergy = kineticEnergy + mass;
340 G4double etot2 = totEnergy*totEnergy; << 318 G4double etot2 = totEnergy*totEnergy;
341 G4double beta2 = kineticEnergy*(kineticEnerg << 319 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
342 320
343 G4double deltaKinEnergy, f; 321 G4double deltaKinEnergy, f;
344 G4double f1 = 0.0; 322 G4double f1 = 0.0;
345 G4double fmax = 1.0; 323 G4double fmax = 1.0;
346 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* 324 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
347 325
348 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 326 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
349 G4double rndm[2]; 327 G4double rndm[2];
350 328
351 // sampling without nuclear size effect 329 // sampling without nuclear size effect
352 do { 330 do {
353 rndmEngineMod->flatArray(2, rndm); 331 rndmEngineMod->flatArray(2, rndm);
354 deltaKinEnergy = minKinEnergy*maxKinEnergy 332 deltaKinEnergy = minKinEnergy*maxKinEnergy
355 /(minKinEnergy*(1.0 - rndm 333 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
356 334
357 f = 1.0 - beta2*deltaKinEnergy/tmax; 335 f = 1.0 - beta2*deltaKinEnergy/tmax;
358 if( 0.0 < spin ) { 336 if( 0.0 < spin ) {
359 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e 337 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
360 f += f1; 338 f += f1;
361 } 339 }
362 340
363 // Loop checking, 03-Aug-2015, Vladimir Iv 341 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
364 } while( fmax*rndm[1] > f); 342 } while( fmax*rndm[1] > f);
365 343
366 // projectile formfactor - suppresion of hig 344 // projectile formfactor - suppresion of high energy
367 // delta-electron production at high energy 345 // delta-electron production at high energy
368 346
369 G4double x = formfact*deltaKinEnergy; 347 G4double x = formfact*deltaKinEnergy;
370 if(x > 1.e-6) { 348 if(x > 1.e-6) {
371 349
372 G4double x1 = 1.0 + x; 350 G4double x1 = 1.0 + x;
373 G4double grej = 1.0/(x1*x1); 351 G4double grej = 1.0/(x1*x1);
374 if( 0.0 < spin ) { 352 if( 0.0 < spin ) {
375 G4double x2 = 0.5*electron_mass_c2*delta 353 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
376 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 354 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
377 } 355 }
378 if(grej > 1.1) { 356 if(grej > 1.1) {
379 G4cout << "### G4LindhardSorensenIonMode 357 G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej
380 << " " << dp->GetDefinition()->G 358 << " " << dp->GetDefinition()->GetParticleName()
381 << " Ekin(MeV)= " << kineticEner 359 << " Ekin(MeV)= " << kineticEnergy
382 << " delEkin(MeV)= " << deltaKinE 360 << " delEkin(MeV)= " << deltaKinEnergy
383 << G4endl; 361 << G4endl;
384 } 362 }
385 if(rndmEngineMod->flat() > grej) { return; 363 if(rndmEngineMod->flat() > grej) { return; }
386 } 364 }
387 365
388 G4ThreeVector deltaDirection; 366 G4ThreeVector deltaDirection;
389 367
390 if(UseAngularGeneratorFlag()) { 368 if(UseAngularGeneratorFlag()) {
391 369
392 const G4Material* mat = couple->GetMateri 370 const G4Material* mat = couple->GetMaterial();
393 G4int Z = SelectRandomAtomNumber(mat); 371 G4int Z = SelectRandomAtomNumber(mat);
394 372
395 deltaDirection = 373 deltaDirection =
396 GetAngularDistribution()->SampleDirectio 374 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
397 375
398 } else { 376 } else {
399 377
400 G4double deltaMomentum = 378 G4double deltaMomentum =
401 std::sqrt(deltaKinEnergy * (deltaKinEner 379 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
402 G4double cost = deltaKinEnergy * (totEnerg 380 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
403 (deltaMomentum * dp->GetTotalMomentum()) 381 (deltaMomentum * dp->GetTotalMomentum());
404 cost = std::min(cost, 1.0); 382 cost = std::min(cost, 1.0);
405 G4double sint = std::sqrt((1.0 - cost)*(1. 383 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
406 384
407 G4double phi = CLHEP::twopi*rndmEngineMod- 385 G4double phi = CLHEP::twopi*rndmEngineMod->flat();
408 386
409 deltaDirection.set(sint*std::cos(phi),sint 387 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
410 deltaDirection.rotateUz(dp->GetMomentumDir 388 deltaDirection.rotateUz(dp->GetMomentumDirection());
411 } 389 }
412 /* 390 /*
413 G4cout << "### G4LindhardSorensenIonModel 391 G4cout << "### G4LindhardSorensenIonModel "
414 << dp->GetDefinition()->GetParticle 392 << dp->GetDefinition()->GetParticleName()
415 << " Ekin(MeV)= " << kineticEnergy 393 << " Ekin(MeV)= " << kineticEnergy
416 << " delEkin(MeV)= " << deltaKinEne 394 << " delEkin(MeV)= " << deltaKinEnergy
417 << " tmin(MeV)= " << minKinEnergy 395 << " tmin(MeV)= " << minKinEnergy
418 << " tmax(MeV)= " << maxKinEnergy 396 << " tmax(MeV)= " << maxKinEnergy
419 << " dir= " << dp->GetMomentumDirec 397 << " dir= " << dp->GetMomentumDirection()
420 << " dirDelta= " << deltaDirection 398 << " dirDelta= " << deltaDirection
421 << G4endl; 399 << G4endl;
422 */ 400 */
423 // create G4DynamicParticle object for delta 401 // create G4DynamicParticle object for delta ray
424 auto delta = new G4DynamicParticle(theElectr << 402 G4DynamicParticle* delta =
>> 403 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
425 404
426 vdp->push_back(delta); 405 vdp->push_back(delta);
427 406
428 // Change kinematics of primary particle 407 // Change kinematics of primary particle
429 kineticEnergy -= deltaKinEnergy; 408 kineticEnergy -= deltaKinEnergy;
430 G4ThreeVector finalP = dp->GetMomentum() - d 409 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
431 finalP = finalP.unit(); << 410 finalP = finalP.unit();
432 411
433 fParticleChange->SetProposedKineticEnergy(ki 412 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
434 fParticleChange->SetProposedMomentumDirectio 413 fParticleChange->SetProposedMomentumDirection(finalP);
435 } 414 }
436 415
437 //....oooOO0OOooo........oooOO0OOooo........oo 416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
438 417
439 G4double 418 G4double
440 G4LindhardSorensenIonModel::MaxSecondaryEnergy 419 G4LindhardSorensenIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
441 420 G4double kinEnergy)
442 { 421 {
443 // here particle type is checked for any met 422 // here particle type is checked for any method
444 SetParticle(pd); 423 SetParticle(pd);
445 G4double tau = kinEnergy/mass; 424 G4double tau = kinEnergy/mass;
446 return 2.0*CLHEP::electron_mass_c2*tau*(tau 425 return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) /
447 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRati 426 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRatio);
448 } 427 }
449 428
450 //....oooOO0OOooo........oooOO0OOooo........oo 429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
451 430