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) {
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) {
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());
232 if(Z > 97) { Z = 97; }
233 G4int nbOfShells = GetNumberOfShells(Z);
234 if(nbOfShells < 1) { nbOfShells = 1; }
238 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/
v;
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; }
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
400 G4double totMomentum = energy*sqrt(beta2);
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); }
434 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
440 const 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};
444 const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] = {1,1,2,3,3,3,3,4,5,4,
448 const 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};
455 const 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
488 const 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
521 const G4double G4ICRU73QOModel::L0[67][2] =
594 const G4double G4ICRU73QOModel::L1[22][2] =
623 const G4double G4ICRU73QOModel::L2[14][2] =
644 const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
645 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979,
646 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96,
647 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821,
648 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481,
649 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459,
650 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828,
651 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221,
652 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226,
653 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676,
654 0.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};