71 mass = charge = chargeSquare = massRate = ratio = 0.0;
72 if(p) { SetParticle(p); }
75 lowestKinEnergy = 5.0*keV;
83 for (
G4int i = 0; i < 100; ++i)
87 for(
G4int i = 0; i < NQOELEM; ++i)
89 if(ZElementAvailable[i] > 0) {
90 indexZ[ZElementAvailable[i]] = i;
93 fParticleChange =
nullptr;
102 if(p != particle) SetParticle(p);
108 isInitialised =
true;
117 denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
131 G4double maxEnergy = std::min(tmax,maxKinEnergy);
132 if(cutEnergy < maxEnergy) {
134 G4double energy = kineticEnergy + mass;
136 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
137 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
138 - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
140 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
156 (p,kineticEnergy,cutEnergy,maxEnergy);
171 (p,kineticEnergy,cutEnergy,maxEnergy);
184 G4double tkin = kineticEnergy/massRate;
186 if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
187 else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
189 if (cutEnergy < tmax) {
197 dedx += chargeSquare*(
G4Log(x) + (1.0 - x)*beta2 ) * twopi_mc2_rcl2
200 if(dedx < 0.0) { dedx = 0.0; }
211 const G4double* theAtomicNumDensityVector =
219 for (
G4int i=0; i<numberOfElements; ++i)
221 const G4Element* element = (*theElementVector)[i] ;
222 eloss += DEDXPerElement(element->
GetZasInt(), kineticEnergy)
223 * theAtomicNumDensityVector[i] * element->
GetZ();
233 G4int Z = std::min(AtomicNumber, 97);
234 G4int nbOfShells = std::max(GetNumberOfShells(Z), 1);
236 G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 );
238 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v;
240 G4double tau = kineticEnergy/proton_mass_c2;
245 G4double l0Term = 0, l1Term = 0, l2Term = 0;
247 for (
G4int nos = 0; nos < nbOfShells; ++nos){
249 G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) /
250 GetShellEnergy(Z,nos);
252 G4double shStrength = GetShellStrength(Z,nos);
254 G4double l0 = GetL0(NormalizedEnergy);
255 l0Term += shStrength * l0;
257 G4double l1 = GetL1(NormalizedEnergy);
258 l1Term += shStrength * l1;
260 G4double l2 = GetL2(NormalizedEnergy);
261 l2Term += shStrength * l2;
264 G4double dedx = 2*CLHEP::twopi_mc2_rcl2*chargeSquare*factorBethe[Z]*
265 (l0Term + charge*fBetheVelocity*l1Term
266 + chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2;
274 G4int nbOfTheShell)
const
280 G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy;
284 * PlasmaEnergy2 / (Z*Z) ;
289 G4double ionTerm2 = ionTerm*ionTerm ;
291 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
293 return oscShellEnergy;
298G4int G4ICRU73QOModel::GetNumberOfShells(
G4int Z)
const
303 nShell = nbofShellsForElement[indexZ[Z]];
316 G4int idx = indexZ[Z];
319 shellEnergy = ShellEnergy[startElemIndex[idx] + nbOfTheShell]*CLHEP::eV;
321 shellEnergy = GetOscillatorEnergy(Z, nbOfTheShell);
333 G4int idx = indexZ[Z];
336 shellStrength = SubShellOccupation[startElemIndex[idx] + nbOfTheShell] / Z;
341 return shellStrength;
350 for(n = 0;
n < sizeL0;
n++) {
351 if( normEnergy < L0[n][0] )
break;
353 if(0 == n) {
n = 1; }
354 if(n >= sizeL0) {
n = sizeL0 - 1; }
358 G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) /
359 (L0[
n][0] - L0[
n-1][0]);
370 for(n = 0;
n < sizeL1;
n++) {
371 if( normEnergy < L1[n][0] )
break;
374 if(n >= sizeL1)
n = sizeL1 - 1 ;
378 G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) /
379 (L1[
n][0] - L1[
n-1][0]);
389 for(n = 0;
n < sizeL2;
n++) {
390 if( normEnergy < L2[n][0] )
break;
393 if(n >= sizeL2)
n = sizeL2 - 1 ;
397 G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) /
398 (L2[
n][0] - L2[
n-1][0]);
421 G4double xmax = std::min(tmax, maxEnergy);
422 if(xmin >= xmax) {
return; }
425 G4double energy = kineticEnergy + mass;
427 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
436 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
438 f = 1.0 - beta2*deltaKinEnergy/tmax;
441 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
442 <<
"Majorant " << grej <<
" < "
443 << f <<
" for e= " << deltaKinEnergy
462 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
463 G4double totMomentum = energy*sqrt(beta2);
464 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
465 (deltaMomentum * totMomentum);
466 if(cost > 1.0) { cost = 1.0; }
467 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
471 deltaDirection.
set(sint*cos(phi),sint*sin(phi), cost) ;
479 kineticEnergy -= deltaKinEnergy;
481 finalP = finalP.
unit();
486 vdp->push_back(delta);
494 if(pd != particle) { SetParticle(pd); }
496 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
497 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
503const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] =
504 {1,2,4,6,7,8,10,13,14,-18,
505 22,26,28,29,32,36,42,47,
506 50,54,73,74,78,79,82,92};
508const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] =
509 {1,1,2,3,3,3,3,4,5,4,
513const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] =
514 {0,1,2,4,7,10,13,16,20,25,
515 29,34,39,44,49,55,59,65,
516 71,78,84,90,98,105,112,121};
521const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
530 1.623, 2.147, 6.259, 2.971,
531 1.631, 2.094, 6.588, 2.041, 1.646,
532 1.535, 8.655, 1.706, 6.104,
533 1.581, 8.358, 8.183, 2.000, 1.878,
534 1.516, 8.325, 8.461, 6.579, 1.119,
535 1.422, 7.81, 8.385, 8.216, 2.167,
536 1.458, 8.049, 8.79, 9.695, 1.008,
537 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
538 1.645, 7.765, 19.192, 7.398,
539 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
540 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
541 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
542 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
543 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
544 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
545 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
546 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
547 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
548 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
554const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
560 732.61, 100.646, 23.550,
561 965.1, 129.85, 31.60,
562 1525.9, 234.9, 56.18,
563 2701, 476.5, 150.42, 16.89,
564 3206.1, 586.4, 186.8, 23.52, 14.91,
565 5551.6, 472.43, 124.85, 22.332,
566 8554.6, 850.58, 93.47, 39.19, 19.46,
567 12254.7, 1279.29, 200.35, 49.19, 17.66,
568 14346.9, 1532.28, 262.71, 74.37, 23.03,
569 15438.5, 1667.96, 294.1, 70.69, 16.447,
570 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
571 24643, 2906.4, 366.85, 22.24,
572 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
573 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
574 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
575 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
576 88926, 18012, 3210, 575, 108.7, 30.8,
577 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
578 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
579 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
580 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
581 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
587const G4double G4ICRU73QOModel::L0[67][2] =
660const G4double G4ICRU73QOModel::L1[22][2] =
689const G4double G4ICRU73QOModel::L2[14][2] =
710const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
7110.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
7120.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
7130.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
7140.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
7150.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
7160.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
7170.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
7180.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
7190.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
7200.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};
std::vector< G4Element * > G4ElementVector
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
std::vector< G4Material * > G4MaterialTable
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)
G4double GetPlasmaEnergy(G4int idx) const
G4int GetElementIndex(G4int Z, G4State mState) const
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
static G4Electron * Electron()
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length) override
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) final
G4ICRU73QOModel(const G4ParticleDefinition *p=0, const G4String &nam="ICRU73QO")
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double) override
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
size_t GetNumberOfElements() const
const G4double * GetAtomicNumDensityVector() const
G4double GetElectronDensity() const
static G4MaterialTable * GetMaterialTable()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4String & GetParticleName() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void SetHighEnergyLimit(G4double)
G4VEmAngularDistribution * GetAngularDistribution()
G4int SelectRandomAtomNumber(const G4Material *)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
G4bool UseAngularGeneratorFlag() const
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4ParticleChangeForLoss * GetParticleChangeForLoss()