Geant4-11
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4LindhardSorensenIonModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// -------------------------------------------------------------------
27//
28// GEANT4 Class header file
29//
30//
31// File name: G4LindhardSorensenIonModel
32//
33// Author: Alexander Bagulya & Vladimir Ivanchenko
34//
35// Creation date: 16.04.2018
36//
37//
38// -------------------------------------------------------------------
39//
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43
45#include "Randomize.hh"
47#include "G4SystemOfUnits.hh"
48#include "G4Electron.hh"
49#include "G4LossTableManager.hh"
50#include "G4EmCorrections.hh"
52#include "G4Log.hh"
53#include "G4DeltaAngle.hh"
55#include "G4BraggIonModel.hh"
56#include "G4BetheBlochModel.hh"
57#include "G4IonICRU73Data.hh"
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61using namespace std;
62
65std::vector<G4float>* G4LindhardSorensenIonModel::fact[] = {nullptr};
66
68 const G4String& nam)
69 : G4VEmModel(nam),
70 particle(nullptr),
71 twoln10(2.0*G4Log(10.0))
72{
73 fParticleChange = nullptr;
79 fElimit = 2.0*CLHEP::MeV;
80}
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
85{}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
90 const G4DataVector& ptr)
91{
92 fBraggModel->Initialise(p, ptr);
93 fBBModel->Initialise(p, ptr);
94 SetParticle(p);
95 //G4cout << "G4LindhardSorensenIonModel::Initialise for "
96 // << p->GetParticleName() << G4endl;
97
98 // always false before the run
100
101 if(nullptr == fParticleChange) {
103 if(UseAngularGeneratorFlag() && nullptr == GetAngularDistribution()) {
105 }
106 }
107 if(IsMaster()) {
108 if(nullptr == lsdata) {
110 }
111 if(nullptr == fIonData) {
113 }
115 }
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119
122 const G4Material* mat,
123 G4double kinEnergy)
124{
125 return corr->EffectiveChargeSquareRatio(p,mat,kinEnergy);
126}
127
128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129
132 const G4Material* mat,
133 G4double kinEnergy)
134{
135 return corr->GetParticleCharge(p,mat,kinEnergy);
136}
137
138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139
141{
145 Zin = G4lrint(std::abs(charge));
149 const G4double aMag =
152 magMoment2 = magmom*magmom - 1.0;
153 G4double x = 0.8426*CLHEP::GeV;
154 if(spin == 0.0 && mass < GeV) { x = 0.736*CLHEP::GeV; }
155 else if (Zin > 1) { x /= nist->GetA27(Zin); }
156
158 tlimit = 2.0/formfact;
159}
160
161//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162
164 const G4MaterialCutsCouple* couple)
165{
167}
168
169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170
173 const G4ParticleDefinition* p,
174 G4double kinEnergy,
175 G4double cutEnergy,
176 G4double maxKinEnergy)
177{
178 // take into account formfactor
179 G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
180 G4double emax = std::min(tmax, maxKinEnergy);
181 G4double escaled = kinEnergy*pRatio;
182 G4double cross = (escaled <= fElimit)
183 ? fBraggModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax)
184 : fBBModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax);
185 // G4cout << "LS: e= " << kinEnergy << " tmin= " << cutEnergy
186 // << " tmax= " << maxEnergy << " cross= " << cross << G4endl;
187 return cross;
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
193 const G4ParticleDefinition* p,
194 G4double kineticEnergy,
196 G4double cutEnergy,
197 G4double maxEnergy)
198{
199 return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
200}
201
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203
205 const G4Material* material,
206 const G4ParticleDefinition* p,
207 G4double kineticEnergy,
208 G4double cutEnergy,
209 G4double maxEnergy)
210{
211 return material->GetElectronDensity()
212 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216
219 const G4ParticleDefinition* p,
220 G4double kinEnergy,
221 G4double cut)
222{
223 // formfactor is taken into account in CorrectionsAlongStep(..)
224 G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
225 G4double cutEnergy = std::min(std::min(cut,tmax), tlimit);
226
227 G4double escaled = kinEnergy*pRatio;
228 G4double dedx = (escaled <= fElimit)
229 ? fBraggModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy)
230 : fBBModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy);
231
232 //G4cout << "E(MeV)=" << kinEnergy/MeV << " dedx=" << dedx
233 // << " " << material->GetName() << " Ecut(MeV)=" << cutEnergy << G4endl;
234 return dedx;
235}
236
237//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
238
240 const G4MaterialCutsCouple* couple,
241 const G4DynamicParticle* dp,
242 const G4double& length,
243 G4double& eloss)
244{
245 // no correction at the last step
246 const G4double preKinEnergy = dp->GetKineticEnergy();
247 if(eloss >= preKinEnergy) { return; }
248
249 const G4ParticleDefinition* p = dp->GetDefinition();
250 const G4Material* mat = couple->GetMaterial();
251 const G4double eDensity = mat->GetElectronDensity();
252 const G4double e = std::max(preKinEnergy - eloss*0.5, preKinEnergy*0.75);
253 const G4double tmax = MaxSecondaryEnergy(p, e);
254 const G4double escaled = e*pRatio;
255 const G4double tau = e/mass;
256
259 const G4int Z = std::max(80, p->GetAtomicNumber());
260
261 G4double res;
262 if(escaled <= fElimit) {
263 res = fIonData->GetDEDX(mat, Z, escaled, G4Log(escaled));
264 if(res > 0.0) {
265 G4double cut = couple->GetProductionCuts()->GetProductionCut(1);
266 if(cut < tmax) {
267 const G4double x = cut/tmax;
268 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x)
269 *q2*CLHEP::twopi_mc2_rcl2*eDensity;
270 }
271 } else {
272 res = eloss*q2*corr->EffectiveChargeCorrection(p,mat,e)/chargeSquare;
273 }
274 } else {
275 const G4double gam = tau + 1.0;
276 const G4double beta2 = tau * (tau+2.0)/(gam*gam);
277 G4double deltaL0 = 2.0*corr->BarkasCorrection(p, mat, e)*(charge-1.)/charge;
278 G4double deltaL = lsdata->GetDeltaL(Zin, gam);
279
280 res = eloss +
281 CLHEP::twopi_mc2_rcl2*q2*eDensity*(deltaL+deltaL0)*length/beta2;
282 /*
283 G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)= "
284 << preKinEnergy/GeV << " eloss(MeV)= " << eloss
285 << " L= " << eloss*beta2/(twopi_mc2_rcl2*chargeSquare*eDensity*length)
286 << " dL0= " << deltaL0
287 << " dL= " << deltaL << G4endl;
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;
295 eloss = res;
296}
297
298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299
301 vector<G4DynamicParticle*>* vdp,
302 const G4MaterialCutsCouple* couple,
303 const G4DynamicParticle* dp,
304 G4double minKinEnergy,
305 G4double maxEnergy)
306{
307 G4double kineticEnergy = dp->GetKineticEnergy();
308 // take into account formfactor
309 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
310
311 G4double maxKinEnergy = std::min(maxEnergy,tmax);
312 if(minKinEnergy >= maxKinEnergy) { return; }
313
314 //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= "
315 // << minKinEnergy << " Emax= " << maxKinEnergy << G4endl;
316
317 G4double totEnergy = kineticEnergy + mass;
318 G4double etot2 = totEnergy*totEnergy;
319 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
320
321 G4double deltaKinEnergy, f;
322 G4double f1 = 0.0;
323 G4double fmax = 1.0;
324 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
325
326 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
327 G4double rndm[2];
328
329 // sampling without nuclear size effect
330 do {
331 rndmEngineMod->flatArray(2, rndm);
332 deltaKinEnergy = minKinEnergy*maxKinEnergy
333 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
334
335 f = 1.0 - beta2*deltaKinEnergy/tmax;
336 if( 0.0 < spin ) {
337 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
338 f += f1;
339 }
340
341 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
342 } while( fmax*rndm[1] > f);
343
344 // projectile formfactor - suppresion of high energy
345 // delta-electron production at high energy
346
347 G4double x = formfact*deltaKinEnergy;
348 if(x > 1.e-6) {
349
350 G4double x1 = 1.0 + x;
351 G4double grej = 1.0/(x1*x1);
352 if( 0.0 < spin ) {
353 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
354 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
355 }
356 if(grej > 1.1) {
357 G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej
358 << " " << dp->GetDefinition()->GetParticleName()
359 << " Ekin(MeV)= " << kineticEnergy
360 << " delEkin(MeV)= " << deltaKinEnergy
361 << G4endl;
362 }
363 if(rndmEngineMod->flat() > grej) { return; }
364 }
365
366 G4ThreeVector deltaDirection;
367
369
370 const G4Material* mat = couple->GetMaterial();
372
373 deltaDirection =
374 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
375
376 } else {
377
378 G4double deltaMomentum =
379 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
380 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
381 (deltaMomentum * dp->GetTotalMomentum());
382 cost = std::min(cost, 1.0);
383 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
384
385 G4double phi = CLHEP::twopi*rndmEngineMod->flat();
386
387 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
388 deltaDirection.rotateUz(dp->GetMomentumDirection());
389 }
390 /*
391 G4cout << "### G4LindhardSorensenIonModel "
392 << dp->GetDefinition()->GetParticleName()
393 << " Ekin(MeV)= " << kineticEnergy
394 << " delEkin(MeV)= " << deltaKinEnergy
395 << " tmin(MeV)= " << minKinEnergy
396 << " tmax(MeV)= " << maxKinEnergy
397 << " dir= " << dp->GetMomentumDirection()
398 << " dirDelta= " << deltaDirection
399 << G4endl;
400 */
401 // create G4DynamicParticle object for delta ray
402 G4DynamicParticle* delta =
403 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
404
405 vdp->push_back(delta);
406
407 // Change kinematics of primary particle
408 kineticEnergy -= deltaKinEnergy;
409 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
410 finalP = finalP.unit();
411
414}
415
416//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417
420 G4double kinEnergy)
421{
422 // here particle type is checked for any method
423 SetParticle(pd);
424 G4double tau = kinEnergy/mass;
425 return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) /
426 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRatio);
427}
428
429//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
static const G4double emax
G4double G4Log(G4double x)
Definition: G4Log.hh:226
static constexpr double GeV
Definition: G4SIunits.hh:203
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
Hep3Vector unit() const
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
virtual double flat()=0
virtual void flatArray(const int size, double *vect)=0
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double BarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetDEDX(const G4Material *, const G4int Z, const G4double e, const G4double loge) const
G4double GetMeanExcitationEnergy() const
G4double GetDeltaL(G4int Z, G4double gamma) const
static std::vector< G4float > * fact[MAXZION]
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) override
G4ParticleChangeForLoss * fParticleChange
static G4LindhardSorensenData * lsdata
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) override
G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
void SetParticle(const G4ParticleDefinition *p)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4LindhardSorensenIonModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="LindhardSorensen")
void CorrectionsAlongStep(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, const G4double &length, G4double &eloss) override
const G4ParticleDefinition * particle
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
static G4LossTableManager * Instance()
G4EmCorrections * EmCorrections()
const G4Material * GetMaterial() const
G4ProductionCuts * GetProductionCuts() const
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:222
G4double GetElectronDensity() const
Definition: G4Material.hh:213
G4double GetA27(G4int Z) const
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double GetPDGMagneticMoment() const
G4int GetAtomicNumber() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
G4double GetProductionCut(G4int index) const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:614
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:621
G4double inveplus
Definition: G4VEmModel.hh:431
G4bool IsMaster() const
Definition: G4VEmModel.hh:746
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.cc:297
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:823
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:628
G4bool UseAngularGeneratorFlag() const
Definition: G4VEmModel.hh:718
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:108
static constexpr double eplus
static constexpr double electron_mass_c2
static constexpr double proton_mass_c2
static constexpr double GeV
static constexpr double hbar_Planck
static constexpr double twopi
Definition: SystemOfUnits.h:56
static constexpr double MeV
static constexpr double c_squared
static constexpr double twopi_mc2_rcl2
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
string material
Definition: eplot.py:19
float electron_mass_c2
Definition: hepunit.py:273
int G4lrint(double ad)
Definition: templates.hh:134