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;
105 if(p != particle) SetParticle(p);
111 isInitialised =
true;
116 denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
130 G4double maxEnergy = std::min(tmax,maxKinEnergy);
131 if(cutEnergy < maxEnergy) {
133 G4double energy = kineticEnergy + mass;
135 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
136 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
138 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
154 (p,kineticEnergy,cutEnergy,maxEnergy);
169 (p,kineticEnergy,cutEnergy,maxEnergy);
182 G4double tkin = kineticEnergy/massRate;
184 if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
185 else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
187 if (cutEnergy < tmax) {
195 dedx += chargeSquare*( log(x) + (1.0 - x)*beta2 ) * twopi_mc2_rcl2
198 if(dedx < 0.0) { dedx = 0.0; }
209 const G4double* theAtomicNumDensityVector =
217 for (
G4int i=0; i<numberOfElements; ++i)
219 const G4Element* element = (*theElementVector)[i] ;
220 eloss += DEDXPerElement(
G4int(element->
GetZ()), kineticEnergy)
221 * theAtomicNumDensityVector[i] *
G4int(element->
GetZ());
231 G4int Z = AtomicNumber;
232 if(Z > 97) { Z = 97; }
233 G4int nbOfShells = GetNumberOfShells(Z);
234 if(nbOfShells < 1) { nbOfShells = 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;
283 * PlasmaEnergy2 / (Z*Z) ;
286 G4double ionTerm2 = ionTerm*ionTerm ;
288 G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
290 return oscShellEnergy;
299 for(n = 0;
n < sizeL0;
n++) {
300 if( normEnergy < L0[n][0] )
break;
302 if(0 == n) {
n = 1; }
303 if(n >= sizeL0) {
n = sizeL0 - 1; }
307 G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) /
308 (L0[
n][0] - L0[
n-1][0]);
319 for(n = 0;
n < sizeL1;
n++) {
320 if( normEnergy < L1[n][0] )
break;
323 if(n >= sizeL1)
n = sizeL1 - 1 ;
327 G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) /
328 (L1[
n][0] - L1[
n-1][0]);
338 for(n = 0;
n < sizeL2;
n++) {
339 if( normEnergy < L2[n][0] )
break;
342 if(n >= sizeL2)
n = sizeL2 - 1 ;
346 G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) /
347 (L2[
n][0] - L2[
n-1][0]);
370 G4double xmax = std::min(tmax, maxEnergy);
371 if(xmin >= xmax) {
return; }
374 G4double energy = kineticEnergy + mass;
376 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
385 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
387 f = 1.0 - beta2*deltaKinEnergy/tmax;
390 G4cout <<
"G4ICRU73QOModel::SampleSecondary Warning! "
391 <<
"Majorant " << grej <<
" < "
392 << f <<
" for e= " << deltaKinEnergy
399 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
400 G4double totMomentum = energy*sqrt(beta2);
401 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
402 (deltaMomentum * totMomentum);
403 if(cost > 1.0) { cost = 1.0; }
404 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
408 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
412 kineticEnergy -= deltaKinEnergy;
413 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
414 finalP = finalP.
unit();
423 vdp->push_back(delta);
431 if(pd != particle) { SetParticle(pd); }
433 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
434 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
440const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] = {1,2,4,6,7,8,10,13,14,-18,
441 22,26,28,29,32,36,42,47,
442 50,54,73,74,78,79,82,92};
444const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] = {1,1,2,3,3,3,3,4,5,4,
448const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] = {0,1,2,4,7,10,13,16,20,25,
449 29,34,39,44,49,55,59,65,
450 71,78,84,90,98,105,112,121};
455const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
464 1.623, 2.147, 6.259, 2.971,
465 1.631, 2.094, 6.588, 2.041, 1.646,
466 1.535, 8.655, 1.706, 6.104,
467 1.581, 8.358, 8.183, 2.000, 1.878,
468 1.516, 8.325, 8.461, 6.579, 1.119,
469 1.422, 7.81, 8.385, 8.216, 2.167,
470 1.458, 8.049, 8.79, 9.695, 1.008,
471 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
472 1.645, 7.765, 19.192, 7.398,
473 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
474 1.295, 6.219, 18.751, 8.748, 10.184, 1.803,
475 1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948,
476 1.563, 6.312, 21.868, 5.762, 11.245, 7.250,
477 0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284,
478 1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108,
479 1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025,
480 1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322,
481 2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000,
482 2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000
488const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
494 732.61, 100.646, 23.550,
495 965.1, 129.85, 31.60,
496 1525.9, 234.9, 56.18,
497 2701, 476.5, 150.42, 16.89,
498 3206.1, 586.4, 186.8, 23.52, 14.91,
499 5551.6, 472.43, 124.85, 22.332,
500 8554.6, 850.58, 93.47, 39.19, 19.46,
501 12254.7, 1279.29, 200.35, 49.19, 17.66,
502 14346.9, 1532.28, 262.71, 74.37, 23.03,
503 15438.5, 1667.96, 294.1, 70.69, 16.447,
504 19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95,
505 24643, 2906.4, 366.85, 22.24,
506 34394, 4365.3, 589.36, 129.42, 35.59, 18.42,
507 43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63,
508 49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54,
509 58987, 8159, 1296.6, 356.75, 101.03, 16.52,
510 88926, 18012, 3210, 575, 108.7, 30.8,
511 115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25,
512 128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04,
513 131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575,
514 154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17,
515 167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43
521const G4double G4ICRU73QOModel::L0[67][2] =
594const G4double G4ICRU73QOModel::L1[22][2] =
623const G4double G4ICRU73QOModel::L2[14][2] =
644const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
6450.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
6460.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
6470.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
6480.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
6490.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
6500.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
6510.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
6520.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
6530.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
6540.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};
std::vector< G4Element * > G4ElementVector
std::vector< G4Material * > G4MaterialTable
G4DLLIMPORT std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
G4double GetPlasmaEnergy(G4int idx)
G4int GetElementIndex(G4int Z, G4State mState)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
virtual ~G4ICRU73QOModel()
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double)
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
G4ICRU73QOModel(const G4ParticleDefinition *p=0, const G4String &nam="ICRU73QO")
static const G4MaterialTable * GetMaterialTable()
const G4ElementVector * GetElementVector() const
size_t GetNumberOfElements() const
const G4double * GetAtomicNumDensityVector() const
G4double GetElectronDensity() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4String & GetParticleName() const
void SetHighEnergyLimit(G4double)
void SetDeexcitationFlag(G4bool val)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
G4ParticleChangeForLoss * GetParticleChangeForLoss()