89 protonMassAMU(1.007276),
133 pname !=
"deuteron" && pname !=
"triton" &&
134 pname !=
"alpha+" && pname !=
"helium" &&
135 pname !=
"hydrogen") {
isIon =
true; }
176 if(cutEnergy < maxEnergy) {
180 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*
mass)/energy2;
181 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
182 - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
184 if( 0.0 <
spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
204 (p,kineticEnergy,cutEnergy,maxEnergy);
219 (p,kineticEnergy,cutEnergy,maxEnergy);
237 dedx =
DEDX(material, tkin);
240 if (cutEnergy < tmax) {
248 dedx += (
G4Log(x) + (1.0 -
x)*beta2) * twopi_mc2_rcl2
254 if (dedx < 0.0) { dedx = 0.0; }
274 if(xmin >= xmax) {
return; }
279 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*
mass)/energy2;
283 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
288 rndmEngineMod->flatArray(2, rndm);
289 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
291 f = 1.0 - beta2*deltaKinEnergy/tmax;
294 G4cout <<
"G4BraggModel::SampleSecondary Warning! "
295 <<
"Majorant " << grej <<
" < "
296 << f <<
" for e= " << deltaKinEnergy
301 }
while( grej*rndm[1] >= f );
315 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
316 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
318 if(cost > 1.0) { cost = 1.0; }
319 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
323 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
332 kineticEnergy -= deltaKinEnergy;
334 finalP = finalP.unit();
339 vdp->push_back(delta);
349 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
398 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
399 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
400 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
401 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
402 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
403 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
404 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
405 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
406 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
407 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
408 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
410 static const G4double atomicWeight[11] = {
411 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
412 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
417 }
else if ( T < 10000.0 ) {
421 ionloss = slow*shigh / (slow + shigh) ;
424 if ( ionloss < 0.0) ionloss = 0.0 ;
429 else if (T < 700.0) {
432 else if (T < 10000.0) {
464 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
465 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
466 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
467 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
468 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
469 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
470 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
471 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
472 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
473 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
475 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
476 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
477 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
478 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
479 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
480 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
481 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
482 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
483 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
484 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
486 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
487 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
488 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
489 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
490 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
491 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
492 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
493 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
494 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
495 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
497 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
498 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
499 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
500 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
501 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
502 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
503 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
504 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
505 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
506 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
508 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
509 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
510 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
511 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
512 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
513 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
515 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2},
516 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
517 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
518 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
520 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
521 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
522 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
523 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
524 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
525 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
526 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
527 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
528 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
529 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
531 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
532 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
533 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
534 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
535 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
536 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
537 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
538 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
539 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
540 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
542 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
543 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
544 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
545 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
546 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
547 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
548 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
549 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
551 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2},
552 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
554 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
555 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
556 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
557 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
558 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
559 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
560 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
561 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
562 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
563 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
565 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
566 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
572 if ( T < 40.0 && 5 == i) {
577 }
else if ( T < 10.0 ) {
584 G4double shigh =
G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
585 ionloss = slow*shigh*fac / (slow + shigh) ;
587 if ( ionloss < 0.0) { ionloss = 0.0; }
611 const G4double* theAtomicNumDensityVector =
624 }
else if(1 == numberOfElements) {
640 for (
G4int i=0; i<numberOfElements; i++) {
641 const G4Element* element = (*theElementVector)[i] ;
644 * theAtomicNumDensityVector[i] ;
646 * theAtomicNumDensityVector[i] ;
658 for (
G4int i=0; i<numberOfElements; i++)
660 const G4Element* element = (*theElementVector)[i] ;
662 * theAtomicNumDensityVector[i];
678 if (myFormula == chFormula ) {
return false; }
688 if( theState ==
kStateGas && myFormula == chFormula)
return false ;
692 static const G4double HeEff = 2.8735 ;
695 static const G4String nameOfMol[53] = {
696 "H_2O",
"C_2H_4O",
"C_3H_6O",
"C_2H_2",
"C_H_3OH",
697 "C_2H_5OH",
"C_3H_7OH",
"C_3H_4",
"NH_3",
"C_14H_10",
698 "C_6H_6",
"C_4H_10",
"C_4H_6",
"C_4H_8O",
"CCl_4",
699 "CF_4",
"C_6H_8",
"C_6H_12",
"C_6H_10O",
"C_6H_10",
700 "C_8H_16",
"C_5H_10",
"C_5H_8",
"C_3H_6-Cyclopropane",
"C_2H_4F_2",
701 "C_2H_2F_2",
"C_4H_8O_2",
"C_2H_6",
"C_2F_6",
"C_2H_6O",
702 "C_3H_6O",
"C_4H_10O",
"C_2H_4",
"C_2H_4O",
"C_2H_4S",
703 "SH_2",
"CH_4",
"CCLF_3",
"CCl_2F_2",
"CHCl_2F",
704 "(CH_3)_2S",
"N_2O",
"C_5H_10O",
"C_8H_6",
"(CH_2)_N",
705 "(C_3H_6)_N",
"(C_8H_8)_N",
"C_3H_8",
"C_3H_6-Propylene",
"C_3H_6O",
706 "C_3H_6S",
"C_4H_4S",
"C_7H_8"
710 66.1, 190.4, 258.7, 42.2, 141.5,
711 210.9, 279.6, 198.8, 31.0, 267.5,
712 122.8, 311.4, 260.3, 328.9, 391.3,
713 206.6, 374.0, 422.0, 432.0, 398.0,
714 554.0, 353.0, 326.0, 74.6, 220.5,
715 197.4, 362.0, 170.0, 330.5, 211.3,
716 262.3, 349.6, 51.3, 187.0, 236.9,
717 121.9, 35.8, 247.0, 292.6, 268.0,
718 262.3, 49.0, 398.9, 444.0, 22.91,
719 68.0, 155.0, 84.0, 74.2, 254.7,
723 static const G4double expCharge[53] = {
724 HeEff, HeEff, HeEff, 1.0, HeEff,
725 HeEff, HeEff, HeEff, 1.0, 1.0,
726 1.0, HeEff, HeEff, HeEff, HeEff,
727 HeEff, HeEff, HeEff, HeEff, HeEff,
728 HeEff, HeEff, HeEff, 1.0, HeEff,
729 HeEff, HeEff, HeEff, HeEff, HeEff,
730 HeEff, HeEff, 1.0, HeEff, HeEff,
731 HeEff, 1.0, HeEff, HeEff, HeEff,
732 HeEff, 1.0, HeEff, HeEff, 1.0,
733 1.0, 1.0, 1.0, 1.0, HeEff,
737 static const G4double numberOfAtomsPerMolecula[53] = {
738 3.0, 7.0, 10.0, 4.0, 6.0,
739 9.0, 12.0, 7.0, 4.0, 24.0,
740 12.0, 14.0, 10.0, 13.0, 5.0,
741 5.0, 14.0, 18.0, 17.0, 17.0,
742 24.0, 15.0, 13.0, 9.0, 8.0,
743 6.0, 14.0, 8.0, 8.0, 9.0,
744 10.0, 15.0, 6.0, 7.0, 7.0,
745 3.0, 5.0, 5.0, 5.0, 5.0,
746 9.0, 3.0, 16.0, 14.0, 3.0,
747 9.0, 16.0, 11.0, 9.0, 10.0,
754 if(chFormula == nameOfMol[i]) {
757 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
775 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
776 G4double gamma25 = 1.0 + 25.0*
keV /proton_mass_c2 ;
777 G4double gamma125 = 1.0 + 125.0*
keV/proton_mass_c2 ;
778 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ;
779 G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
780 G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
783 (1.0 +
G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
784 (1.0 +
G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
std::vector< G4Element * > G4ElementVector
G4bool HasMaterial(const G4Material *material)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
void SetExpStopPower125(G4double value)
static const G4int numberOfMolecula
const G4String & GetChemicalFormula() const
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4VEmAngularDistribution * GetAngularDistribution()
G4double GetDensity() const
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetElectronicDEDX(G4int idx, G4double energy) const
G4VEmFluctuationModel * GetModelOfFluctuations()
G4double GetChargeSquareRatio() const
const G4ElementVector * GetElementVector() const
const G4String & GetParticleName() const
static G4PSTARStopping * fPSTAR
G4ParticleChangeForLoss * fParticleChange
void SetHighEnergyLimit(G4double)
G4double GetTotalMomentum() const
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition * theElectron
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4double GetElectronDensity() const
G4bool UseAngularGeneratorFlag() const
G4EmCorrections * EmCorrections()
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
const G4String & GetParticleType() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4double * GetAtomicNumDensityVector() const
const G4Material * currentMaterial
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double GetTotNbOfAtomsPerVolume() const
G4double StoppingPower(const G4Material *material, G4double kineticEnergy)
const G4ParticleDefinition * particle
static const G4double factor
G4double GetPDGMass() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double energy(const ThreeVector &p, const G4double m)
static const G4double invLog10
void SetAngularDistribution(G4VEmAngularDistribution *)
const G4double x[NPOINTSGL]
void SetParticle(const G4ParticleDefinition *p)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double ChemicalFactor(G4double kineticEnergy, G4double eloss125) const
static G4Electron * Electron()
static const G4double fac
G4int GetIndex(const G4Material *) const
size_t GetNumberOfElements() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double theZieglerFactor
G4BraggModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="Bragg")
void SetDeexcitationFlag(G4bool val)
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4bool MolecIsInZiegler1988(const G4Material *material)
G4double DEDX(const G4Material *material, G4double kineticEnergy)
G4ThreeVector GetMomentum() const
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)
G4double ElectronicStoppingPower(G4double z, G4double kineticEnergy) const