83 for (
G4int i = 0; i < 100; ++i)
118 denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
133 if(cutEnergy < maxEnergy) {
137 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*
mass)/energy2;
138 cross = 1.0/cutEnergy - 1.0/maxEnergy
139 - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
157 (p,kineticEnergy,cutEnergy,maxEnergy);
172 (p,kineticEnergy,cutEnergy,maxEnergy);
190 if (cutEnergy < tmax) {
201 if(dedx < 0.0) { dedx = 0.0; }
212 const G4double* theAtomicNumDensityVector =
220 for (
G4int i=0; i<numberOfElements; ++i)
222 const G4Element* element = (*theElementVector)[i] ;
224 * theAtomicNumDensityVector[i] *
G4int(element->
GetZ());
234 G4int Z = AtomicNumber;
235 if(Z > 97) { Z = 97; }
237 if(nbOfShells < 1) { nbOfShells = 1; }
239 G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 );
241 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v;
243 G4double tau = kineticEnergy/proton_mass_c2;
248 G4double l0Term = 0, l1Term = 0, l2Term = 0;
250 for (
G4int nos = 0; nos < nbOfShells; ++nos){
252 G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) /
258 l0Term += shStrength * l0;
261 l1Term += shStrength * l1;
264 l2Term += shStrength * l2;
268 (l0Term +
charge*fBetheVelocity*l1Term
269 +
chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2;
277 G4int nbOfTheShell)
const
283 G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy;
287 * PlasmaEnergy2 / (Z*Z) ;
291 G4double ionTerm2 = ionTerm*ionTerm ;
293 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
295 return oscShellEnergy;
304 for(n = 0; n <
sizeL0; n++) {
305 if( normEnergy <
L0[n][0] )
break;
307 if(0 == n) { n = 1; }
308 if(n >= sizeL0) { n = sizeL0 - 1; }
312 G4double bethe = l0p + (l0 - l0p) * ( normEnergy -
L0[n-1][0]) /
313 (
L0[
n][0] -
L0[n-1][0]);
324 for(n = 0; n <
sizeL1; n++) {
325 if( normEnergy <
L1[n][0] )
break;
328 if(n >= sizeL1) n = sizeL1 - 1 ;
332 G4double barkas= l1p + (l1 - l1p) * ( normEnergy -
L1[n-1][0]) /
333 (
L1[
n][0] -
L1[n-1][0]);
343 for(n = 0; n <
sizeL2; n++) {
344 if( normEnergy <
L2[n][0] )
break;
347 if(n >= sizeL2) n = sizeL2 - 1 ;
351 G4double bloch = l2p + (l2 - l2p) * ( normEnergy -
L2[n-1][0]) /
352 (
L2[
n][0] -
L2[n-1][0]);
376 if(xmin >= xmax) {
return; }
381 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*
mass)/energy2;
390 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
392 f = 1.0 - beta2*deltaKinEnergy/tmax;
395 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
396 <<
"Majorant " << grej <<
" < "
397 << f <<
" for e= " << deltaKinEnergy
404 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
405 G4double totMomentum = energy*sqrt(beta2);
406 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
407 (deltaMomentum * totMomentum);
408 if(cost > 1.0) { cost = 1.0; }
409 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
413 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
414 deltaDirection.rotateUz(direction);
417 kineticEnergy -= deltaKinEnergy;
418 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
419 finalP = finalP.unit();
428 vdp->push_back(delta);
438 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
445 const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] = {1,2,4,6,7,8,10,13,14,-18,
446 22,26,28,29,32,36,42,47,
447 50,54,73,74,78,79,82,92};
449 const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] = {1,1,2,3,3,3,3,4,5,4,
453 const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] = {0,1,2,4,7,10,13,16,20,25,
454 29,34,39,44,49,55,59,65,
455 71,78,84,90,98,105,112,121};
469 1.623, 2.147, 6.259, 2.971,
470 1.631, 2.094, 6.588, 2.041, 1.646,
471 1.535, 8.655, 1.706, 6.104,
472 1.581, 8.358, 8.183, 2.000, 1.878,
473 1.516, 8.325, 8.461, 6.579, 1.119,
474 1.422, 7.81, 8.385, 8.216, 2.167,
475 1.458, 8.049, 8.79, 9.695, 1.008,
476 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
477 1.645, 7.765, 19.192, 7.398,
478 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
479 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
480 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
481 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
482 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
483 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
484 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
485 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
486 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
487 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
499 732.61, 100.646, 23.550,
500 965.1, 129.85, 31.60,
501 1525.9, 234.9, 56.18,
502 2701, 476.5, 150.42, 16.89,
503 3206.1, 586.4, 186.8, 23.52, 14.91,
504 5551.6, 472.43, 124.85, 22.332,
505 8554.6, 850.58, 93.47, 39.19, 19.46,
506 12254.7, 1279.29, 200.35, 49.19, 17.66,
507 14346.9, 1532.28, 262.71, 74.37, 23.03,
508 15438.5, 1667.96, 294.1, 70.69, 16.447,
509 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
510 24643, 2906.4, 366.85, 22.24,
511 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
512 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
513 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
514 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
515 88926, 18012, 3210, 575, 108.7, 30.8,
516 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
517 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
518 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
519 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
520 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
650 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
651 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
652 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
653 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
654 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
655 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
656 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
657 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
658 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
659 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
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 G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double GetElectronDensity() 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)
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)