84 for (
G4int i = 0; i < 100; ++i)
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);
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] ;
229 * theAtomicNumDensityVector[i] *
G4int(element->
GetZ());
239 G4int Z = AtomicNumber;
240 if(Z > 97) { Z = 97; }
242 if(nbOfShells < 1) { nbOfShells = 1; }
244 G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 );
246 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v;
248 G4double tau = kineticEnergy/proton_mass_c2;
253 G4double l0Term = 0, l1Term = 0, l2Term = 0;
255 for (
G4int nos = 0; nos < nbOfShells; ++nos){
257 G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) /
263 l0Term += shStrength * l0;
266 l1Term += shStrength * l1;
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) ;
296 G4double ionTerm2 = ionTerm*ionTerm ;
298 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
300 return oscShellEnergy;
309 for(n = 0; n <
sizeL0; n++) {
310 if( normEnergy <
L0[n][0] )
break;
312 if(0 == n) { n = 1; }
313 if(n >= sizeL0) { n = sizeL0 - 1; }
317 G4double bethe = l0p + (l0 - l0p) * ( normEnergy -
L0[n-1][0]) /
318 (
L0[
n][0] -
L0[n-1][0]);
329 for(n = 0; n <
sizeL1; n++) {
330 if( normEnergy <
L1[n][0] )
break;
333 if(n >= sizeL1) n = sizeL1 - 1 ;
337 G4double barkas= l1p + (l1 - l1p) * ( normEnergy -
L1[n-1][0]) /
338 (
L1[
n][0] -
L1[n-1][0]);
348 for(n = 0; n <
sizeL2; n++) {
349 if( normEnergy <
L2[n][0] )
break;
352 if(n >= sizeL2) n = sizeL2 - 1 ;
356 G4double bloch = l2p + (l2 - l2p) * ( normEnergy -
L2[n-1][0]) /
357 (
L2[
n][0] -
L2[n-1][0]);
381 if(xmin >= xmax) {
return; }
386 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*
mass)/energy2;
395 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
397 f = 1.0 - beta2*deltaKinEnergy/tmax;
400 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
401 <<
"Majorant " << grej <<
" < "
402 << f <<
" for e= " << deltaKinEnergy
420 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
421 G4double totMomentum = energy*sqrt(beta2);
422 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
423 (deltaMomentum * totMomentum);
424 if(cost > 1.0) { cost = 1.0; }
425 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
429 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
430 deltaDirection.rotateUz(direction);
437 kineticEnergy -= deltaKinEnergy;
439 finalP = finalP.unit();
444 vdp->push_back(delta);
454 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
461 const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] = {1,2,4,6,7,8,10,13,14,-18,
462 22,26,28,29,32,36,42,47,
463 50,54,73,74,78,79,82,92};
465 const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] = {1,1,2,3,3,3,3,4,5,4,
469 const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] = {0,1,2,4,7,10,13,16,20,25,
470 29,34,39,44,49,55,59,65,
471 71,78,84,90,98,105,112,121};
485 1.623, 2.147, 6.259, 2.971,
486 1.631, 2.094, 6.588, 2.041, 1.646,
487 1.535, 8.655, 1.706, 6.104,
488 1.581, 8.358, 8.183, 2.000, 1.878,
489 1.516, 8.325, 8.461, 6.579, 1.119,
490 1.422, 7.81, 8.385, 8.216, 2.167,
491 1.458, 8.049, 8.79, 9.695, 1.008,
492 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
493 1.645, 7.765, 19.192, 7.398,
494 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
495 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
496 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
497 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
498 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
499 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
500 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
501 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
502 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
503 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
515 732.61, 100.646, 23.550,
516 965.1, 129.85, 31.60,
517 1525.9, 234.9, 56.18,
518 2701, 476.5, 150.42, 16.89,
519 3206.1, 586.4, 186.8, 23.52, 14.91,
520 5551.6, 472.43, 124.85, 22.332,
521 8554.6, 850.58, 93.47, 39.19, 19.46,
522 12254.7, 1279.29, 200.35, 49.19, 17.66,
523 14346.9, 1532.28, 262.71, 74.37, 23.03,
524 15438.5, 1667.96, 294.1, 70.69, 16.447,
525 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
526 24643, 2906.4, 366.85, 22.24,
527 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
528 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
529 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
530 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
531 88926, 18012, 3210, 575, 108.7, 30.8,
532 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
533 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
534 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
535 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
536 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
666 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
667 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
668 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
669 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
670 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
671 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
672 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
673 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
674 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
675 0.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
static const G4int nbofShellsForElement[NQOELEM]
static const G4double factorBethe[99]
G4double DEDXPerElement(G4int Z, G4double kineticEnergy)
G4double GetL1(G4double normEnergy) const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
std::vector< G4Element * > G4ElementVector
static const G4double SubShellOccupation[NQODATA]
G4double DEDX(const G4Material *material, G4double kineticEnergy)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
static const G4double L0[67][2]
G4int GetElementIndex(G4int Z, G4State mState)
static const G4double L2[14][2]
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
G4VEmAngularDistribution * GetAngularDistribution()
G4double GetL2(G4double normEnergy) const
static const G4double L1[22][2]
G4ParticleChangeForLoss * fParticleChange
const G4ElementVector * GetElementVector() const
G4ICRU73QOModel(const G4ParticleDefinition *p=0, const G4String &nam="ICRU73QO")
const G4String & GetParticleName() const
G4double GetShellStrength(G4int Z, G4int nbOfTheShell) const
void SetHighEnergyLimit(G4double)
G4ParticleDefinition * theElectron
G4double GetShellEnergy(G4int Z, G4int nbOfTheShell) const
G4GLOB_DLL std::ostream G4cout
G4DensityEffectData * denEffData
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double GetElectronDensity() const
G4bool UseAngularGeneratorFlag() const
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double GetPlasmaEnergy(G4int idx)
void SetParticle(const G4ParticleDefinition *p)
virtual ~G4ICRU73QOModel()
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double GetOscillatorEnergy(G4int Z, G4int nbOfTheShell) const
const G4double * GetAtomicNumDensityVector() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double G4Log(G4double x)
static const G4int NQOELEM
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double energy(const ThreeVector &p, const G4double m)
void SetAngularDistribution(G4VEmAngularDistribution *)
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()
size_t GetNumberOfElements() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
void SetDeexcitationFlag(G4bool val)
static const G4double ShellEnergy[NQODATA]
G4int GetNumberOfShells(G4int Z) const
static const G4int ZElementAvailable[NQOELEM]
const G4ParticleDefinition * particle
static const G4int startElemIndex[NQOELEM]
G4double GetL0(G4double normEnergy) const
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)