71 mass = charge = chargeSquare = massRate = ratio = 0.0;
72 if(p) { SetParticle(p); }
75 lowestKinEnergy = 5.0*
keV;
83 for (
G4int i = 0; i < 100; ++i)
87 for(
G4int i = 0; i < NQOELEM; ++i)
89 if(ZElementAvailable[i] > 0) {
90 indexZ[ZElementAvailable[i]] = i;
107 if(p != particle) SetParticle(p);
113 isInitialised =
true;
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;
141 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
157 (p,kineticEnergy,cutEnergy,maxEnergy);
172 (p,kineticEnergy,cutEnergy,maxEnergy);
185 G4double tkin = kineticEnergy/massRate;
187 if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
188 else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
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] ;
223 eloss += DEDXPerElement(
G4int(element->
GetZ()), kineticEnergy)
224 * theAtomicNumDensityVector[i] *
G4int(element->
GetZ());
234 G4int Z = AtomicNumber;
235 if(Z > 97) { Z = 97; }
236 G4int nbOfShells = GetNumberOfShells(Z);
237 if(nbOfShells < 1) { nbOfShells = 1; }
241 G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/
v;
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) /
253 GetShellEnergy(Z,nos);
255 G4double shStrength = GetShellStrength(Z,nos);
257 G4double l0 = GetL0(NormalizedEnergy);
258 l0Term += shStrength * l0;
260 G4double l1 = GetL1(NormalizedEnergy);
261 l1Term += shStrength * l1;
263 G4double l2 = GetL2(NormalizedEnergy);
264 l2Term += shStrength * l2;
267 G4double dedx = 2*CLHEP::twopi_mc2_rcl2*chargeSquare*factorBethe[Z]*
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
405 G4double totMomentum = energy*sqrt(beta2);
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) ;
417 kineticEnergy -= deltaKinEnergy;
418 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
419 finalP = finalP.
unit();
428 vdp->push_back(delta);
436 if(pd != particle) { SetParticle(pd); }
439 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
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};
460 const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
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
493 const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
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
526 const G4double G4ICRU73QOModel::L0[67][2] =
599 const G4double G4ICRU73QOModel::L1[22][2] =
628 const G4double G4ICRU73QOModel::L2[14][2] =
649 const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
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)
G4ParticleChangeForLoss * GetParticleChangeForLoss()
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
G4int GetElementIndex(G4int Z, G4State mState)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
const G4ElementVector * GetElementVector() const
G4ICRU73QOModel(const G4ParticleDefinition *p=0, const G4String &nam="ICRU73QO")
const G4String & GetParticleName() const
void SetHighEnergyLimit(G4double)
double precision function energy(A, Z)
G4GLOB_DLL std::ostream G4cout
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
Hep3Vector & rotateUz(const Hep3Vector &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double GetPlasmaEnergy(G4int idx)
virtual ~G4ICRU73QOModel()
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4double * GetAtomicNumDensityVector() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
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)
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)