72 mass = charge = chargeSquare = massRate = ratio = 0.0;
73 if(p) { SetParticle(p); }
76 lowestKinEnergy = 5.0*
keV;
84 for (
G4int i = 0; i < 100; ++i)
88 for(
G4int i = 0; i < NQOELEM; ++i)
90 if(ZElementAvailable[i] > 0) {
91 indexZ[ZElementAvailable[i]] = i;
94 fParticleChange =
nullptr;
108 if(p != particle) SetParticle(p);
114 isInitialised =
true;
123 denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
138 if(cutEnergy < maxEnergy) {
142 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
143 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
144 - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
162 (p,kineticEnergy,cutEnergy,maxEnergy);
177 (p,kineticEnergy,cutEnergy,maxEnergy);
190 G4double tkin = kineticEnergy/massRate;
192 if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
193 else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
195 if (cutEnergy < tmax) {
206 if(dedx < 0.0) { dedx = 0.0; }
217 const G4double* theAtomicNumDensityVector =
225 for (
G4int i=0; i<numberOfElements; ++i)
227 const G4Element* element = (*theElementVector)[i] ;
228 eloss += DEDXPerElement(element->
GetZasInt(), kineticEnergy)
229 * theAtomicNumDensityVector[i] * element->
GetZ();
240 if(Z > 97) { Z = 97; }
241 G4int nbOfShells = GetNumberOfShells(Z);
242 if(nbOfShells < 1) { nbOfShells = 1; }
253 G4double l0Term = 0, l1Term = 0, l2Term = 0;
255 for (
G4int nos = 0; nos < nbOfShells; ++nos){
258 GetShellEnergy(Z,nos);
260 G4double shStrength = GetShellStrength(Z,nos);
262 G4double l0 = GetL0(NormalizedEnergy);
263 l0Term += shStrength * l0;
265 G4double l1 = GetL1(NormalizedEnergy);
266 l1Term += shStrength * l1;
268 G4double l2 = GetL2(NormalizedEnergy);
269 l2Term += shStrength * l2;
273 (l0Term + charge*fBetheVelocity*l1Term
274 + chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2;
282 G4int nbOfTheShell)
const
288 G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy;
292 * PlasmaEnergy2 / (Z*
Z) ;
297 G4double ionTerm2 = ionTerm*ionTerm ;
299 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
301 return oscShellEnergy;
310 for(n = 0; n < sizeL0; n++) {
311 if( normEnergy < L0[n][0] )
break;
313 if(0 == n) { n = 1; }
314 if(n >= sizeL0) { n = sizeL0 - 1; }
318 G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) /
319 (L0[
n][0] - L0[n-1][0]);
330 for(n = 0; n < sizeL1; n++) {
331 if( normEnergy < L1[n][0] )
break;
334 if(n >= sizeL1) n = sizeL1 - 1 ;
338 G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) /
339 (L1[
n][0] - L1[n-1][0]);
349 for(n = 0; n < sizeL2; n++) {
350 if( normEnergy < L2[n][0] )
break;
353 if(n >= sizeL2) n = sizeL2 - 1 ;
357 G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) /
358 (L2[
n][0] - L2[n-1][0]);
382 if(xmin >= xmax) {
return; }
387 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
396 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
398 f = 1.0 - beta2*deltaKinEnergy/tmax;
401 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
402 <<
"Majorant " << grej <<
" < "
403 << f <<
" for e= " << deltaKinEnergy
423 G4double totMomentum = energy*sqrt(beta2);
425 (deltaMomentum * totMomentum);
426 if(cost > 1.0) { cost = 1.0; }
427 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
431 deltaDirection.
set(sint*cos(phi),sint*sin(phi), cost) ;
439 kineticEnergy -= deltaKinEnergy;
441 finalP = finalP.
unit();
446 vdp->push_back(delta);
454 if(pd != particle) { SetParticle(pd); }
457 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
463 const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] =
464 {1,2,4,6,7,8,10,13,14,-18,
465 22,26,28,29,32,36,42,47,
466 50,54,73,74,78,79,82,92};
468 const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] =
469 {1,1,2,3,3,3,3,4,5,4,
473 const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] =
474 {0,1,2,4,7,10,13,16,20,25,
475 29,34,39,44,49,55,59,65,
476 71,78,84,90,98,105,112,121};
481 const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
490 1.623, 2.147, 6.259, 2.971,
491 1.631, 2.094, 6.588, 2.041, 1.646,
492 1.535, 8.655, 1.706, 6.104,
493 1.581, 8.358, 8.183, 2.000, 1.878,
494 1.516, 8.325, 8.461, 6.579, 1.119,
495 1.422, 7.81, 8.385, 8.216, 2.167,
496 1.458, 8.049, 8.79, 9.695, 1.008,
497 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
498 1.645, 7.765, 19.192, 7.398,
499 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
500 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
501 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
502 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
503 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
504 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
505 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
506 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
507 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
508 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
514 const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
520 732.61, 100.646, 23.550,
521 965.1, 129.85, 31.60,
522 1525.9, 234.9, 56.18,
523 2701, 476.5, 150.42, 16.89,
524 3206.1, 586.4, 186.8, 23.52, 14.91,
525 5551.6, 472.43, 124.85, 22.332,
526 8554.6, 850.58, 93.47, 39.19, 19.46,
527 12254.7, 1279.29, 200.35, 49.19, 17.66,
528 14346.9, 1532.28, 262.71, 74.37, 23.03,
529 15438.5, 1667.96, 294.1, 70.69, 16.447,
530 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
531 24643, 2906.4, 366.85, 22.24,
532 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
533 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
534 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
535 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
536 88926, 18012, 3210, 575, 108.7, 30.8,
537 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
538 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
539 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
540 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
541 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
547 const G4double G4ICRU73QOModel::L0[67][2] =
620 const G4double G4ICRU73QOModel::L1[22][2] =
649 const G4double G4ICRU73QOModel::L2[14][2] =
670 const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
671 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
672 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
673 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
674 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
675 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
676 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
677 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
678 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
679 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
680 0.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};
void set(double x, double y, double z)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
G4ParticleChangeForLoss * GetParticleChangeForLoss()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double) override
std::vector< G4Element * > G4ElementVector
static constexpr double proton_mass_c2
G4double GetKineticEnergy() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
static constexpr double twopi_mc2_rcl2
G4double GetPlasmaEnergy(G4int idx) const
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
G4VEmAngularDistribution * GetAngularDistribution()
G4int GetElementIndex(G4int Z, G4State mState) const
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length) override
const G4ElementVector * GetElementVector() const
G4ICRU73QOModel(const G4ParticleDefinition *p=0, const G4String &nam="ICRU73QO")
const G4String & GetParticleName() const
static constexpr double twopi
static constexpr double electron_mass_c2
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) final
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4double GetElectronDensity() const
G4bool UseAngularGeneratorFlag() const
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual ~G4ICRU73QOModel()
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4double * GetAtomicNumDensityVector() const
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double energy(const ThreeVector &p, const G4double m)
void SetAngularDistribution(G4VEmAngularDistribution *)
static constexpr double c_light
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
static constexpr double MeV
size_t GetNumberOfElements() const
void SetDeexcitationFlag(G4bool val)
static constexpr double fine_structure_const
static constexpr double keV
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)