87 theIonEffChargeModel(0),
88 theNuclearStoppingModel(0),
89 theIonChuFluctuationModel(0),
90 theIonYangFluctuationModel(0),
91 protonTable(
"ICRU_R49p"),
92 antiprotonTable(
"ICRU_R49p"),
93 theNuclearTable(
"ICRU_R49"),
96 theMeanFreePathTable(0),
97 paramStepLimit (0.005),
98 pixeCrossSectionHandler(0)
124 G4String defaultPixeModel(
"ecpssr");
125 modelK = defaultPixeModel;
126 modelL = defaultPixeModel;
127 modelM = defaultPixeModel;
198 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable for "
209 G4cout <<
" 0: " << (*pv)[0]->GetProcessName() <<
" " << (*pv)[0]
210 <<
" 1: " << (*pv)[1]->GetProcessName() <<
" " << (*pv)[1]
255 for (
size_t j=0; j<numOfCouples; j++) {
307 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: "
308 <<
"Loss table is built "
321 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: "
322 <<
"DEDX table will be built "
337 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: end for "
350 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
355 if (particleDef == *proton)
382 for (
size_t j=0; j<numOfCouples; j++) {
403 paramB = ionloss/ionlossBB - 1.0 ;
410 if ( lowEdgeEnergy < highEnergy ) {
426 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
431 G4cout <<
"E(MeV)= " << lowEdgeEnergy/
MeV
432 <<
" dE/dx(MeV/mm)= " << ionloss*
mm/
MeV
452 G4cout <<
"G4hImpactIonisation::BuildLambdaTable for "
476 for (
size_t j=0 ; j < numOfCouples; j++) {
503 for (
G4int iel=0; iel<numberOfElements; iel++ )
505 Z = (
G4int) (*theElementVector)[iel]->GetZ();
513 sigma += theAtomicNumDensityVector[iel] * microCross ;
518 value = sigma<=0 ?
DBL_MAX : 1./sigma ;
550 G4double gamma = energy / particleMass;
551 G4double beta2 = 1. - 1. / (gamma * gamma);
557 if ( tMax > deltaCutInEnergy )
559 var = deltaCutInEnergy / tMax;
560 totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
567 totalCrossSection += 0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
570 else if (spin > 0.9 )
572 totalCrossSection += -std::log(var) /
573 (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
574 beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
582 return totalCrossSection ;
598 G4bool isOutRange =
false;
613 meanFreePath = (((*theMeanFreePathTable)(couple->
GetIndex()))->
617 return meanFreePath ;
641 G4double tScaled = kineticEnergy*massRatio ;
698 if(tScaled > highEnergy )
710 if (stepLimit > x) stepLimit =
x;
737 G4double tScaled = kineticEnergy * massRatio ;
745 eLoss = kineticEnergy ;
751 eLoss = stepLength *
fdEdx ;
756 eLoss = kineticEnergy ;
789 eLoss = stepLength *
fdEdx ;
797 if (eLoss < 0.0) eLoss = 0.0;
799 finalT = kineticEnergy - eLoss - nLoss;
806 if (eLoss < 0.0) eLoss = 0.0;
807 finalT = kineticEnergy - eLoss - nLoss;
822 eLoss = kineticEnergy-finalT;
852 G4cout <<
"p E(MeV)= " << kineticEnergy/
MeV
853 <<
" dE/dx(MeV/mm)= " << eLoss*
mm/
MeV
854 <<
" for " << material->
GetName()
858 if(eLoss < 0.0) eLoss = 0.0 ;
903 G4cout <<
"pbar E(MeV)= " << kineticEnergy/
MeV
904 <<
" dE/dx(MeV/mm)= " << eLoss*
mm/
MeV
905 <<
" for " << material->
GetName()
909 if(eLoss < 0.0) eLoss = 0.0 ;
925 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
927 G4double tau = kineticEnergy / particleMass ;
934 G4double beta2 = bg2/(gamma*gamma) ;
940 if ( deltaCut < tMax)
943 dLoss = ( beta2 * (x-1.) - std::log(x) ) *
twopi_mc2_rcl2 * electronDensity / beta2 ;
968 G4double pSquare = kineticEnergy *( totalEnergy +
mass) ;
969 G4double eSquare = totalEnergy * totalEnergy;
970 G4double betaSquare = pSquare / eSquare;
973 G4double gamma = kineticEnergy / mass + 1.;
981 if (deltaCut >= tMax)
998 gRej = 1.0 - betaSquare *
x ;
1000 else if (0.5 == spin)
1002 gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
1006 gRej = (1. - betaSquare *
x ) * (1. + x/(3.*xc)) +
1007 x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1008 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
1013 G4double deltaKineticEnergy = x * tMax;
1014 G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy *
1016 G4double totalMomentum = std::sqrt(pSquare) ;
1020 if ( cosTheta < -1. ) cosTheta = -1.;
1021 if ( cosTheta > 1. ) cosTheta = 1.;
1025 G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
1026 G4double dirX = sinTheta * std::cos(phi);
1027 G4double dirY = sinTheta * std::sin(phi);
1031 deltaDirection.
rotateUz(particleDirection);
1038 deltaDirection.
z());
1042 G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1043 size_t totalNumber = 1;
1049 size_t nSecondaries = 0;
1050 std::vector<G4DynamicParticle*>* secondaryVector = 0;
1111 if (finalKineticEnergy >= bindingEnergy)
1127 if (secondaryVector != 0)
1129 nSecondaries = secondaryVector->size();
1130 for (
size_t i = 0; i<nSecondaries; i++)
1148 if (e < finalKineticEnergy &&
1153 finalKineticEnergy -=
e;
1180 (*secondaryVector)[i] = 0;
1195 G4double finalPx = totalMomentum*particleDirection.
x() - deltaTotalMomentum*deltaDirection.
x();
1196 G4double finalPy = totalMomentum*particleDirection.
y() - deltaTotalMomentum*deltaDirection.
y();
1197 G4double finalPz = totalMomentum*particleDirection.
z() - deltaTotalMomentum*deltaDirection.
z();
1198 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1199 finalPx /= finalMomentum;
1200 finalPy /= finalMomentum;
1201 finalPz /= finalMomentum;
1207 eDeposit = finalKineticEnergy;
1208 finalKineticEnergy = 0.;
1210 particleDirection.
y(),
1211 particleDirection.
z());
1213 GetAtRestProcessVector()->size())
1236 if (secondaryVector != 0)
1238 for (
size_t l = 0; l < nSecondaries; l++)
1257 delete secondaryVector;
1365 if(0.5*
MeV > kinE) kinE = 0.5*
MeV ;
1367 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1368 if(0.0 >= beta2)
return 0.0;
1376 for (
G4int i = 0; i<numberOfElements; i++) {
1379 ZMaterial = (*theElementVector)[i]->GetZ();
1381 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1385 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1386 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1387 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1389 for(
G4int j=0; j<47; j++) {
1391 if( W < FTable[j][0] ) {
1394 FunctionOfW = FTable[0][1] ;
1397 FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1398 / (FTable[j][0] - FTable[j-1][0])
1407 BTerm += FunctionOfW /( std::sqrt(ZMaterial * X) *
X);
1429 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1430 G4double y = cSquare / (137.0*137.0*beta2) ;
1436 eLoss = 1.0 / (1.0 +
y) ;
1439 for(
G4int i=2; de>eLoss*0.01; i++) {
1440 de = 1.0/( i * (i*i +
y)) ;
1445 (
material->GetElectronDensity()) / beta2 ;
1467 static const G4double kappa = 10. ;
1471 G4int imaterial = couple->GetIndex() ;
1479 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1482 if(meanLoss < minLoss)
return meanLoss ;
1495 if(tMax > threshold) tMax = threshold;
1499 if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
1502 * electronDensity / beta2 ;
1505 if( beta2 > 3.0*theBohrBeta2*zeff ||
charge < 0.0) {
1506 siga = std::sqrt( siga * chargeSquare ) ;
1513 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1518 }
while (loss < 0. || loss > 2.0*meanLoss);
1523 static const G4double probLim = 0.01 ;
1524 static const G4double sumaLim = -std::log(probLim) ;
1531 G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1543 w1 = tMax/ipotFluct;
1546 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1548 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1549 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1550 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
1551 if(a1 < 0.0) a1 = 0.0;
1552 if(a2 < 0.0) a2 = 0.0;
1553 if(a3 < 0.0) a3 = 0.0;
1564 if(tMax == ipotFluct)
1570 siga=std::sqrt(a3) ;
1584 tMax = tMax-ipotFluct+e0 ;
1585 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
1589 siga=std::sqrt(a3) ;
1597 w = (tMax-e0)/tMax ;
1601 corrfac = dp3/
G4float(nmaxCont2) ;
1608 loss *= e0*corrfac ;
1618 siga=std::sqrt(a1) ;
1627 siga=std::sqrt(a2) ;
1633 loss = p1*e1Fluct+p2*e2Fluct;
1646 siga=std::sqrt(a3) ;
1660 rfac = dp3/(
G4float(nmaxCont2)+dp3);
1666 alfa = w1*
G4float(nmaxCont2+p3)/
1668 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1669 ea = na*ipotFluct*alfa1;
1670 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1678 w2 = alfa*ipotFluct;
1715 G4String comments =
" Knock-on electron cross sections . ";
1716 comments +=
"\n Good description above the mean excitation energy.\n";
1717 comments +=
" Delta ray energy sampled from differential Xsection.";
1722 <<
" in " <<
TotBin <<
" bins."
1723 <<
"\n Electronic stopping power model is "
1727 G4cout <<
"\n Parametrisation model for antiprotons is "
1732 G4cout <<
" Parametrization of the Barkas effect is switched on."
1748 for (
size_t j=0 ; j < numOfCouples; j++) {
1755 if(excitationEnergy > deltaCutNow) {
1759 G4cout <<
" material min.delta energy(keV) " <<
G4endl;
1764 << std::setw(15) << excitationEnergy/
keV <<
G4endl;