69 mass = charge = chargeSquare = massRate = ratio = 0.0;
70 if(p) { SetParticle(p); }
73 lowestKinEnergy = 5.0*keV;
81 for (
G4int i = 0; i < 100; ++i)
85 for(
G4int i = 0; i < NQOELEM; ++i)
87 if(ZElementAvailable[i] > 0) {
88 indexZ[ZElementAvailable[i]] = i;
91 fParticleChange =
nullptr;
100 if(p != particle) SetParticle(p);
106 isInitialised =
true;
112 G4String pname = particle->GetParticleName();
115 denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
129 const G4double maxEnergy = std::min(tmax, maxKinEnergy);
130 const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate);
131 if(cutEnergy < maxEnergy) {
133 const G4double energy = kineticEnergy + mass;
134 const G4double energy2 = energy*energy;
135 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
136 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
137 - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
139 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
155 (p,kineticEnergy,cutEnergy,maxEnergy);
170 (p,kineticEnergy,cutEnergy,maxEnergy);
183 const G4double tkin = kineticEnergy/massRate;
184 const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate);
186 if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
187 else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
189 if (cutEnergy < tmax) {
191 const G4double tau = kineticEnergy/mass;
193 dedx += (
G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x) *
196 dedx = std::max(dedx, 0.0);
207 const G4double* theAtomicNumDensityVector =
215 for (std::size_t i=0; i<numberOfElements; ++i)
217 const G4Element* element = (*theElementVector)[i] ;
218 eloss += DEDXPerElement(element->
GetZasInt(), kineticEnergy)
219 * theAtomicNumDensityVector[i] * element->
GetZ();
229 G4int Z = std::min(AtomicNumber, 97);
230 G4int nbOfShells = std::max(GetNumberOfShells(Z), 1);
232 G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 );
234 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v;
236 G4double tau = kineticEnergy/proton_mass_c2;
241 G4double l0Term = 0, l1Term = 0, l2Term = 0;
243 for (
G4int nos = 0; nos < nbOfShells; ++nos){
245 G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) /
246 GetShellEnergy(Z,nos);
248 G4double shStrength = GetShellStrength(Z,nos);
250 G4double l0 = GetL0(NormalizedEnergy);
251 l0Term += shStrength * l0;
253 G4double l1 = GetL1(NormalizedEnergy);
254 l1Term += shStrength * l1;
256 G4double l2 = GetL2(NormalizedEnergy);
257 l2Term += shStrength * l2;
260 G4double dedx = 2*CLHEP::twopi_mc2_rcl2*chargeSquare*factorBethe[Z]*
261 (l0Term + charge*fBetheVelocity*l1Term
262 + chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2;
270 G4int nbOfTheShell)
const
273 if(idx == -1) { idx = denEffData->GetElementIndex(Z-1,
kStateUndefined); }
274 G4double PlasmaEnergy = denEffData->GetPlasmaEnergy(idx);
276 G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy;
280 * PlasmaEnergy2 / (Z*Z) ;
285 G4double ionTerm2 = ionTerm*ionTerm ;
287 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
289 return oscShellEnergy;
294G4int G4ICRU73QOModel::GetNumberOfShells(
G4int Z)
const
299 nShell = nbofShellsForElement[indexZ[Z]];
312 G4int idx = indexZ[Z];
315 shellEnergy = ShellEnergy[startElemIndex[idx] + nbOfTheShell]*CLHEP::eV;
317 shellEnergy = GetOscillatorEnergy(Z, nbOfTheShell);
329 G4int idx = indexZ[Z];
332 shellStrength = SubShellOccupation[startElemIndex[idx] + nbOfTheShell] / Z;
337 return shellStrength;
346 for(n = 0;
n < sizeL0;
n++) {
347 if( normEnergy < L0[n][0] )
break;
349 if(0 == n) {
n = 1; }
350 if(n >= sizeL0) {
n = sizeL0 - 1; }
354 G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) /
355 (L0[
n][0] - L0[
n-1][0]);
366 for(n = 0;
n < sizeL1;
n++) {
367 if( normEnergy < L1[n][0] )
break;
370 if(n >= sizeL1)
n = sizeL1 - 1 ;
374 G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) /
375 (L1[
n][0] - L1[
n-1][0]);
385 for(n = 0;
n < sizeL2;
n++) {
386 if( normEnergy < L2[n][0] )
break;
389 if(n >= sizeL2)
n = sizeL2 - 1 ;
393 G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) /
394 (L2[
n][0] - L2[
n-1][0]);
408 const G4double xmax = std::min(tmax, maxEnergy);
409 const G4double xmin = std::max(minEnergy, lowestKinEnergy*massRate);
410 if(xmin >= xmax) {
return; }
413 const G4double energy = kineticEnergy + mass;
414 const G4double energy2 = energy*energy;
415 const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
424 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
426 f = 1.0 - beta2*deltaKinEnergy/tmax;
429 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
430 <<
"Majorant " << grej <<
" < "
431 << f <<
" for e= " << deltaKinEnergy
450 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
451 G4double totMomentum = energy*sqrt(beta2);
452 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
453 (deltaMomentum * totMomentum);
454 if(cost > 1.0) { cost = 1.0; }
455 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
459 deltaDirection.
set(sint*cos(phi),sint*sin(phi), cost) ;
466 kineticEnergy -= deltaKinEnergy;
468 finalP = finalP.
unit();
470 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
471 fParticleChange->SetProposedMomentumDirection(finalP);
473 vdp->push_back(delta);
481 if(pd != particle) { SetParticle(pd); }
483 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
484 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
490const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] =
491 {1,2,4,6,7,8,10,13,14,-18,
492 22,26,28,29,32,36,42,47,
493 50,54,73,74,78,79,82,92};
495const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] =
496 {1,1,2,3,3,3,3,4,5,4,
500const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] =
501 {0,1,2,4,7,10,13,16,20,25,
502 29,34,39,44,49,55,59,65,
503 71,78,84,90,98,105,112,121};
508const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
517 1.623, 2.147, 6.259, 2.971,
518 1.631, 2.094, 6.588, 2.041, 1.646,
519 1.535, 8.655, 1.706, 6.104,
520 1.581, 8.358, 8.183, 2.000, 1.878,
521 1.516, 8.325, 8.461, 6.579, 1.119,
522 1.422, 7.81, 8.385, 8.216, 2.167,
523 1.458, 8.049, 8.79, 9.695, 1.008,
524 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
525 1.645, 7.765, 19.192, 7.398,
526 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
527 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
528 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
529 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
530 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
531 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
532 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
533 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
534 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
535 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
541const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
547 732.61, 100.646, 23.550,
548 965.1, 129.85, 31.60,
549 1525.9, 234.9, 56.18,
550 2701, 476.5, 150.42, 16.89,
551 3206.1, 586.4, 186.8, 23.52, 14.91,
552 5551.6, 472.43, 124.85, 22.332,
553 8554.6, 850.58, 93.47, 39.19, 19.46,
554 12254.7, 1279.29, 200.35, 49.19, 17.66,
555 14346.9, 1532.28, 262.71, 74.37, 23.03,
556 15438.5, 1667.96, 294.1, 70.69, 16.447,
557 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
558 24643, 2906.4, 366.85, 22.24,
559 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
560 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
561 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
562 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
563 88926, 18012, 3210, 575, 108.7, 30.8,
564 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
565 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
566 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
567 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
568 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
574const G4double G4ICRU73QOModel::L0[67][2] =
647const G4double G4ICRU73QOModel::L1[22][2] =
676const G4double G4ICRU73QOModel::L2[14][2] =
697const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
6980.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
6990.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
7000.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
7010.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
7020.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
7030.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
7040.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
7050.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
7060.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
7070.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};
std::vector< const G4Element * > G4ElementVector
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
std::vector< G4Material * > G4MaterialTable
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
static G4int GetNumberOfShells(G4int Z)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
static G4Electron * Electron()
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
G4ICRU73QOModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="ICRU73QO")
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double) override
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
const G4double * GetAtomicNumDensityVector() const
G4double GetElectronDensity() const
static G4MaterialTable * GetMaterialTable()
std::size_t GetNumberOfElements() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void SetHighEnergyLimit(G4double)
G4VEmAngularDistribution * GetAngularDistribution()
G4int SelectRandomAtomNumber(const G4Material *) const
void SetDeexcitationFlag(G4bool val)
G4VEmModel(const G4String &nam)
void SetAngularDistribution(G4VEmAngularDistribution *)
G4bool UseAngularGeneratorFlag() const
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4ParticleChangeForLoss * GetParticleChangeForLoss()