88 protonMassAMU(1.007276),
97 lowestKinEnergy = 1.0*
keV;
98 theZieglerFactor =
eV*
cm2*1.0e-15;
100 expStopPower125 = 0.0;
103 if(p) { SetParticle(p); }
104 else { SetParticle(theElectron); }
117 if(p != particle) { SetParticle(p); }
123 isInitialised =
true;
130 pname !=
"deuteron" && pname !=
"triton" &&
131 pname !=
"alpha+" && pname !=
"helium" &&
132 pname !=
"hydrogen") { isIon =
true; }
171 if(cutEnergy < maxEnergy) {
175 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
176 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*
G4Log(maxEnergy/cutEnergy)/tmax;
196 (p,kineticEnergy,cutEnergy,maxEnergy);
211 (p,kineticEnergy,cutEnergy,maxEnergy);
223 G4double tkin = kineticEnergy/massRate;
226 if(tkin < lowestKinEnergy) {
227 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
229 dedx = DEDX(material, tkin);
232 if (cutEnergy < tmax) {
246 if (dedx < 0.0) { dedx = 0.0; }
248 dedx *= chargeSquare;
266 if(xmin >= xmax) {
return; }
271 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
278 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
280 f = 1.0 - beta2*deltaKinEnergy/tmax;
283 G4cout <<
"G4BraggModel::SampleSecondary Warning! "
284 <<
"Majorant " << grej <<
" < "
285 << f <<
" for e= " << deltaKinEnergy
306 if(cost > 1.0) { cost = 1.0; }
307 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
311 deltaDirection.
set(sint*cos(phi),sint*sin(phi), cost) ;
320 kineticEnergy -= deltaKinEnergy;
322 finalP = finalP.
unit();
327 vdp->push_back(delta);
335 if(pd != particle) { SetParticle(pd); }
338 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
347 if(
"" == chFormula) {
return false; }
350 static const size_t numberOfMolecula = 11;
351 static const G4String molName[numberOfMolecula] = {
352 "Al_2O_3",
"CO_2",
"CH_4",
353 "(C_2H_4)_N-Polyethylene",
"(C_2H_4)_N-Polypropylene",
"(C_8H_8)_N",
354 "C_3H_8",
"SiO_2",
"H_2O",
355 "H_2O-Gas",
"Graphite" } ;
356 static const G4int idxPSTAR[numberOfMolecula] = {
365 if( theState ==
kStateGas &&
"H_2O" == chFormula) {
370 for (
size_t i=0; i<numberOfMolecula; ++i) {
371 if (chFormula == molName[i]) {
373 iPSTAR = idxPSTAR[i];
387 if (iMolecula >= 0) {
396 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
397 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
398 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
399 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
400 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
401 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
402 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
403 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
404 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
405 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
406 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
408 static const G4double atomicWeight[11] = {
409 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
410 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
413 ionloss = a[iMolecula][0] * sqrt(T) ;
415 }
else if ( T < 10000.0 ) {
418 + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
419 ionloss = slow*shigh / (slow + shigh) ;
422 if ( ionloss < 0.0) ionloss = 0.0 ;
423 if ( 10 == iMolecula ) {
425 ionloss *= (1.0+0.023+0.0066*
G4Log(T)*invLog10);
427 else if (T < 700.0) {
428 ionloss *=(1.0+0.089-0.0248*
G4Log(T-99.)*invLog10);
430 else if (T < 10000.0) {
431 ionloss *=(1.0+0.089-0.0248*
G4Log(700.-99.)*invLog10);
434 ionloss /= atomicWeight[iMolecula];
439 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
462 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
463 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
464 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
465 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
466 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
467 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
468 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
469 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
470 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
471 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
473 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
474 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
475 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
476 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
477 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
478 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
479 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
480 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
481 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
482 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
484 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
485 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
486 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
487 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
488 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
489 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
490 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
491 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
492 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
493 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
495 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
496 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
497 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
498 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
499 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
500 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
501 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
502 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
503 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
504 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
506 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
507 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
508 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
509 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
510 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
511 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
513 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2},
514 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
515 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
516 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
518 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
519 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
520 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
521 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
522 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
523 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
524 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
525 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
526 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
527 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
529 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
530 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
531 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
532 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
533 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
534 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
535 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
536 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
537 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
538 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
540 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
541 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
542 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
543 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
544 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
545 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
546 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
547 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
549 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2},
550 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
552 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
553 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
554 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
555 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
556 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
557 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
558 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
559 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
560 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
561 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
563 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
564 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
570 if ( T < 40.0 && 5 == i) {
575 }
else if ( T < 10.0 ) {
582 G4double shigh =
G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
583 ionloss = slow*shigh*fac / (slow + shigh) ;
585 if ( ionloss < 0.0) { ionloss = 0.0; }
597 if(material != currentMaterial) {
601 if( !HasMaterial(material) ) { iPSTAR = pstar.
GetIndex(material); }
609 const G4double* theAtomicNumDensityVector =
615 }
else if(iMolecula >= 0) {
617 eloss = StoppingPower(material, kineticEnergy)*
621 }
else if(1 == numberOfElements) {
624 eloss = ElectronicStoppingPower(z, kineticEnergy)
629 }
else if( MolecIsInZiegler1988(material) ) {
637 for (
G4int i=0; i<numberOfElements; i++) {
638 const G4Element* element = (*theElementVector)[i] ;
640 eloss += ElectronicStoppingPower(z,kineticEnergy)
641 * theAtomicNumDensityVector[i] ;
642 eloss125 += ElectronicStoppingPower(z,125.0*
keV)
643 * theAtomicNumDensityVector[i] ;
647 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
655 for (
G4int i=0; i<numberOfElements; i++)
657 const G4Element* element = (*theElementVector)[i] ;
658 eloss += ElectronicStoppingPower(element->
GetZ(), kineticEnergy)
659 * theAtomicNumDensityVector[i];
662 return eloss*theZieglerFactor;
675 if (myFormula == chFormula ) {
return false; }
685 if( theState ==
kStateGas && myFormula == chFormula)
return false ;
689 static const G4double HeEff = 2.8735 ;
691 static const size_t numberOfMolecula = 53;
692 static const G4String nameOfMol[53] = {
693 "H_2O",
"C_2H_4O",
"C_3H_6O",
"C_2H_2",
"C_H_3OH",
694 "C_2H_5OH",
"C_3H_7OH",
"C_3H_4",
"NH_3",
"C_14H_10",
695 "C_6H_6",
"C_4H_10",
"C_4H_6",
"C_4H_8O",
"CCl_4",
696 "CF_4",
"C_6H_8",
"C_6H_12",
"C_6H_10O",
"C_6H_10",
697 "C_8H_16",
"C_5H_10",
"C_5H_8",
"C_3H_6-Cyclopropane",
"C_2H_4F_2",
698 "C_2H_2F_2",
"C_4H_8O_2",
"C_2H_6",
"C_2F_6",
"C_2H_6O",
699 "C_3H_6O",
"C_4H_10O",
"C_2H_4",
"C_2H_4O",
"C_2H_4S",
700 "SH_2",
"CH_4",
"CCLF_3",
"CCl_2F_2",
"CHCl_2F",
701 "(CH_3)_2S",
"N_2O",
"C_5H_10O",
"C_8H_6",
"(CH_2)_N",
702 "(C_3H_6)_N",
"(C_8H_8)_N",
"C_3H_8",
"C_3H_6-Propylene",
"C_3H_6O",
703 "C_3H_6S",
"C_4H_4S",
"C_7H_8"
706 static const G4double expStopping[numberOfMolecula] = {
707 66.1, 190.4, 258.7, 42.2, 141.5,
708 210.9, 279.6, 198.8, 31.0, 267.5,
709 122.8, 311.4, 260.3, 328.9, 391.3,
710 206.6, 374.0, 422.0, 432.0, 398.0,
711 554.0, 353.0, 326.0, 74.6, 220.5,
712 197.4, 362.0, 170.0, 330.5, 211.3,
713 262.3, 349.6, 51.3, 187.0, 236.9,
714 121.9, 35.8, 247.0, 292.6, 268.0,
715 262.3, 49.0, 398.9, 444.0, 22.91,
716 68.0, 155.0, 84.0, 74.2, 254.7,
720 static const G4double expCharge[53] = {
721 HeEff, HeEff, HeEff, 1.0, HeEff,
722 HeEff, HeEff, HeEff, 1.0, 1.0,
723 1.0, HeEff, HeEff, HeEff, HeEff,
724 HeEff, HeEff, HeEff, HeEff, HeEff,
725 HeEff, HeEff, HeEff, 1.0, HeEff,
726 HeEff, HeEff, HeEff, HeEff, HeEff,
727 HeEff, HeEff, 1.0, HeEff, HeEff,
728 HeEff, 1.0, HeEff, HeEff, HeEff,
729 HeEff, 1.0, HeEff, HeEff, 1.0,
730 1.0, 1.0, 1.0, 1.0, HeEff,
734 static const G4double numberOfAtomsPerMolecula[53] = {
735 3.0, 7.0, 10.0, 4.0, 6.0,
736 9.0, 12.0, 7.0, 4.0, 24.0,
737 12.0, 14.0, 10.0, 13.0, 5.0,
738 5.0, 14.0, 18.0, 17.0, 17.0,
739 24.0, 15.0, 13.0, 9.0, 8.0,
740 6.0, 14.0, 8.0, 8.0, 9.0,
741 10.0, 15.0, 6.0, 7.0, 7.0,
742 3.0, 5.0, 5.0, 5.0, 5.0,
743 9.0, 3.0, 16.0, 14.0, 3.0,
744 9.0, 16.0, 11.0, 9.0, 10.0,
749 for (
size_t i=0; i<numberOfMolecula; i++)
751 if(chFormula == nameOfMol[i]) {
754 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
755 SetExpStopPower125(exp125);
775 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ;
776 G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
777 G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
779 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
780 (1.0 +
G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
781 (1.0 +
G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
void set(double x, double y, double z)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
G4BraggModel(const G4ParticleDefinition *p=0, const G4String &nam="Bragg")
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
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)
G4VEmFluctuationModel * GetModelOfFluctuations()
G4double GetChargeSquareRatio() const
const G4ElementVector * GetElementVector() const
const G4String & GetParticleName() const
void SetHighEnergyLimit(G4double)
G4double GetTotalMomentum() const
G4int GetIndex(const G4Material *)
G4GLOB_DLL std::ostream G4cout
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
Hep3Vector & rotateUz(const Hep3Vector &)
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
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double GetTotNbOfAtomsPerVolume() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double GetElectronicDEDX(G4int idx, G4double energy)
void SetAngularDistribution(G4VEmAngularDistribution *)
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)
static G4Electron * Electron()
size_t GetNumberOfElements() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void SetDeexcitationFlag(G4bool val)
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
G4ThreeVector GetMomentum() const
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)