88 theIonEffChargeModel(0),
89 theNuclearStoppingModel(0),
90 theIonChuFluctuationModel(0),
91 theIonYangFluctuationModel(0),
92 protonTable(
"ICRU_R49p"),
93 antiprotonTable(
"ICRU_R49p"),
94 theNuclearTable(
"ICRU_R49"),
97 theMeanFreePathTable(0),
98 paramStepLimit (0.005),
99 pixeCrossSectionHandler(0)
125 G4String defaultPixeModel(
"ecpssr");
126 modelK = defaultPixeModel;
127 modelL = defaultPixeModel;
128 modelM = defaultPixeModel;
199 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable for "
210 G4cout <<
" 0: " << (*pv)[0]->GetProcessName() <<
" " << (*pv)[0]
211 <<
" 1: " << (*pv)[1]->GetProcessName() <<
" " << (*pv)[1]
256 for (
size_t j=0; j<numOfCouples; j++) {
308 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: "
309 <<
"Loss table is built "
322 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: "
323 <<
"DEDX table will be built "
338 G4cout <<
"G4hImpactIonisation::BuildPhysicsTable: end for "
351 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
356 if (particleDef == *proton)
383 for (
size_t j=0; j<numOfCouples; j++) {
404 paramB = ionloss/ionlossBB - 1.0 ;
411 if ( lowEdgeEnergy < highEnergy ) {
427 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
432 G4cout <<
"E(MeV)= " << lowEdgeEnergy/
MeV
433 <<
" dE/dx(MeV/mm)= " << ionloss*
mm/
MeV
453 G4cout <<
"G4hImpactIonisation::BuildLambdaTable for "
477 for (
size_t j=0 ; j < numOfCouples; j++) {
504 for (
G4int iel=0; iel<numberOfElements; iel++ )
506 Z = (
G4int) (*theElementVector)[iel]->GetZ();
514 sigma += theAtomicNumDensityVector[iel] * microCross ;
519 value = sigma<=0 ?
DBL_MAX : 1./sigma ;
551 G4double gamma = energy / particleMass;
552 G4double beta2 = 1. - 1. / (gamma * gamma);
553 G4double var = electron_mass_c2 / particleMass;
554 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
558 if ( tMax > deltaCutInEnergy )
560 var = deltaCutInEnergy / tMax;
561 totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
568 totalCrossSection += 0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
571 else if (spin > 0.9 )
573 totalCrossSection += -std::log(var) /
574 (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
575 beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
577 totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
583 return totalCrossSection ;
599 G4bool isOutRange =
false;
614 meanFreePath = (((*theMeanFreePathTable)(couple->
GetIndex()))->
618 return meanFreePath ;
642 G4double tScaled = kineticEnergy*massRatio ;
699 if(tScaled > highEnergy )
711 if (stepLimit > x) stepLimit = x;
738 G4double tScaled = kineticEnergy * massRatio ;
746 eLoss = kineticEnergy ;
752 eLoss = stepLength *
fdEdx ;
757 eLoss = kineticEnergy ;
790 eLoss = stepLength *
fdEdx ;
798 if (eLoss < 0.0) eLoss = 0.0;
800 finalT = kineticEnergy - eLoss - nLoss;
807 if (eLoss < 0.0) eLoss = 0.0;
808 finalT = kineticEnergy - eLoss - nLoss;
823 eLoss = kineticEnergy-finalT;
853 G4cout <<
"p E(MeV)= " << kineticEnergy/
MeV
854 <<
" dE/dx(MeV/mm)= " << eLoss*
mm/
MeV
855 <<
" for " << material->
GetName()
859 if(eLoss < 0.0) eLoss = 0.0 ;
904 G4cout <<
"pbar E(MeV)= " << kineticEnergy/
MeV
905 <<
" dE/dx(MeV/mm)= " << eLoss*
mm/
MeV
906 <<
" for " << material->
GetName()
910 if(eLoss < 0.0) eLoss = 0.0 ;
926 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
928 G4double tau = kineticEnergy / particleMass ;
929 G4double rateMass = electron_mass_c2/particleMass ;
935 G4double beta2 = bg2/(gamma*gamma) ;
936 G4double tMax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
941 if ( deltaCut < tMax)
944 dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
968 G4double totalEnergy = kineticEnergy + mass ;
969 G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
970 G4double eSquare = totalEnergy * totalEnergy;
971 G4double betaSquare = pSquare / eSquare;
974 G4double gamma = kineticEnergy / mass + 1.;
975 G4double r = electron_mass_c2 / mass;
976 G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
982 if (deltaCut >= tMax)
999 gRej = 1.0 - betaSquare * x ;
1001 else if (0.5 == spin)
1003 gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
1007 gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
1008 x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1009 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
1014 G4double deltaKineticEnergy = x * tMax;
1015 G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy *
1016 (deltaKineticEnergy + 2. * electron_mass_c2 ));
1017 G4double totalMomentum = std::sqrt(pSquare) ;
1018 G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
1021 if ( cosTheta < -1. ) cosTheta = -1.;
1022 if ( cosTheta > 1. ) cosTheta = 1.;
1026 G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
1027 G4double dirX = sinTheta * std::cos(phi);
1028 G4double dirY = sinTheta * std::sin(phi);
1032 deltaDirection.rotateUz(particleDirection);
1039 deltaDirection.z());
1043 G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1044 size_t totalNumber = 1;
1050 size_t nSecondaries = 0;
1051 std::vector<G4DynamicParticle*>* secondaryVector = 0;
1112 if (finalKineticEnergy >= bindingEnergy)
1128 if (secondaryVector != 0)
1130 nSecondaries = secondaryVector->size();
1131 for (
size_t i = 0; i<nSecondaries; i++)
1149 if (e < finalKineticEnergy &&
1154 finalKineticEnergy -= e;
1181 (*secondaryVector)[i] = 0;
1196 G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
1197 G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
1198 G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
1199 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1200 finalPx /= finalMomentum;
1201 finalPy /= finalMomentum;
1202 finalPz /= finalMomentum;
1208 eDeposit = finalKineticEnergy;
1209 finalKineticEnergy = 0.;
1211 particleDirection.y(),
1212 particleDirection.z());
1214 GetAtRestProcessVector()->size())
1237 if (secondaryVector != 0)
1239 for (
size_t l = 0; l < nSecondaries; l++)
1258 delete secondaryVector;
1366 if(0.5*
MeV > kinE) kinE = 0.5*
MeV ;
1367 G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1368 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1369 if(0.0 >= beta2)
return 0.0;
1374 const G4ElementVector* theElementVector = material->GetElementVector();
1375 G4int numberOfElements = material->GetNumberOfElements();
1377 for (
G4int i = 0; i<numberOfElements; i++) {
1380 ZMaterial = (*theElementVector)[i]->GetZ();
1382 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1386 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1387 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1388 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1390 for(
G4int j=0; j<47; j++) {
1392 if( W < FTable[j][0] ) {
1395 FunctionOfW = FTable[0][1] ;
1398 FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1399 / (FTable[j][0] - FTable[j-1][0])
1408 BTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1411 BTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1429 G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1430 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1431 G4double y = cSquare / (137.0*137.0*beta2) ;
1437 eLoss = 1.0 / (1.0 + y) ;
1440 for(
G4int i=2; de>eLoss*0.01; i++) {
1441 de = 1.0/( i * (i*i + y)) ;
1445 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1446 (material->GetElectronDensity()) / beta2 ;
1468 static const G4double kappa = 10. ;
1469 static const G4double theBohrBeta2 = 50.0 *
keV/proton_mass_c2 ;
1478 G4double tkin = particle->GetKineticEnergy();
1479 G4double particleMass = particle->GetMass() ;
1480 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1483 if(meanLoss < minLoss)
return meanLoss ;
1489 G4double rmass = electron_mass_c2/particleMass;
1493 G4double tMax = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1496 if(tMax > threshold) tMax = threshold;
1500 if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
1502 siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1503 * electronDensity / beta2 ;
1506 if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1507 siga = std::sqrt( siga * chargeSquare ) ;
1512 G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1513 G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1514 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1519 }
while (loss < 0. || loss > 2.0*meanLoss);
1524 static const G4double probLim = 0.01 ;
1525 static const G4double sumaLim = -std::log(probLim) ;
1532 G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1544 w1 = tMax/ipotFluct;
1545 w2 = std::log(2.*electron_mass_c2*tau2);
1547 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1549 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1550 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1551 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
1552 if(a1 < 0.0) a1 = 0.0;
1553 if(a2 < 0.0) a2 = 0.0;
1554 if(a3 < 0.0) a3 = 0.0;
1565 if(tMax == ipotFluct)
1571 siga=std::sqrt(a3) ;
1585 tMax = tMax-ipotFluct+e0 ;
1586 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
1590 siga=std::sqrt(a3) ;
1598 w = (tMax-e0)/tMax ;
1602 corrfac = dp3/
G4float(nmaxCont2) ;
1609 loss *= e0*corrfac ;
1619 siga=std::sqrt(a1) ;
1628 siga=std::sqrt(a2) ;
1634 loss = p1*e1Fluct+p2*e2Fluct;
1647 siga=std::sqrt(a3) ;
1661 rfac = dp3/(
G4float(nmaxCont2)+dp3);
1667 alfa = w1*
G4float(nmaxCont2+p3)/
1669 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1670 ea = na*ipotFluct*alfa1;
1671 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1679 w2 = alfa*ipotFluct;
1716 G4String comments =
" Knock-on electron cross sections . ";
1717 comments +=
"\n Good description above the mean excitation energy.\n";
1718 comments +=
" Delta ray energy sampled from differential Xsection.";
1723 <<
" in " <<
TotBin <<
" bins."
1724 <<
"\n Electronic stopping power model is "
1728 G4cout <<
"\n Parametrisation model for antiprotons is "
1733 G4cout <<
" Parametrization of the Barkas effect is switched on."
1749 for (
size_t j=0 ; j < numOfCouples; j++) {
1756 if(excitationEnergy > deltaCutNow) {
1760 G4cout <<
" material min.delta energy(keV) " <<
G4endl;
1765 << std::setw(15) << excitationEnergy/
keV <<
G4endl;
G4double condition(const G4ErrorSymMatrix &m)
G4ParticleDefinition * GetDefinition() const
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
G4double antiprotonLowEnergy
void BuildLossTable(const G4ParticleDefinition &aParticleType)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
void SetProtonElectronicStoppingPowerModel(const G4String &dedxTable)
static G4ThreadLocal G4PhysicsTable * theProperTimepTable
static G4ThreadLocal G4double LowestKineticEnergy
G4long G4Poisson(G4double mean)
static G4ThreadLocal G4PhysicsTable * theLabTimepTable
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
std::vector< G4Element * > G4ElementVector
void insert(G4PhysicsVector *)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void ActivateAugerElectronProduction(G4bool val)
G4double GetStepLength() const
const G4DynamicParticle * GetDynamicParticle() const
G4double GetEnergy2fluct() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4double DeltaRaysEnergy(const G4MaterialCutsCouple *couple, G4double kineticEnergy, G4double particleMass) const
G4double GetConstraints(const G4DynamicParticle *particle, const G4MaterialCutsCouple *couple)
G4double GetProductionCut(G4int index) const
const G4String & GetName() const
G4VLowEnergyModel * theNuclearStoppingModel
static G4Proton * ProtonDefinition()
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double GetLogEnergy2fluct() const
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double protonHighEnergy
G4double AntiProtonParametrisedDEDX(const G4MaterialCutsCouple *couple, G4double kineticEnergy) const
G4ParticleDefinition * GetDefinition() const
static G4ThreadLocal G4PhysicsTable ** RecorderOfpbarProcess
G4double BindingEnergy() const
G4double GetLogMeanExcEnergy() const
G4hImpactIonisation(const G4String &processName="hImpactIoni")
const G4String & GetParticleSubType() const
static G4ThreadLocal G4double finalRange
G4double GetLowEdgeEnergy(size_t binNumber) const
G4VLowEnergyModel * theIonEffChargeModel
G4ProcessManager * GetProcessManager() const
const G4ElementVector * GetElementVector() const
const G4String & GetParticleName() const
const G4double paramStepLimit
G4VLowEnergyModel * theIonChuFluctuationModel
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double MicroscopicCrossSection(const G4ParticleDefinition &aParticleType, G4double kineticEnergy, G4double atomicNumber, G4double deltaCutInEnergy) const
virtual G4bool IsInCharge(const G4DynamicParticle *particle, const G4Material *material) const =0
G4VLowEnergyModel * theIonYangFluctuationModel
static G4ThreadLocal G4int TotBin
G4double GetEnergy0fluct() const
G4double MinKineticEnergy
static G4ThreadLocal G4PhysicsTable ** RecorderOfpProcess
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
void SetCutForSecondaryPhotons(G4double cut)
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
static G4ThreadLocal G4PhysicsTable * theRangepTable
void BuildLambdaTable(const G4ParticleDefinition &aParticleType)
G4double GetElectronDensity() const
void InitializeParametrisation()
static G4AntiProton * AntiProton()
G4PixeCrossSectionHandler * pixeCrossSectionHandler
G4AtomicDeexcitation atomicDeexcitation
void SetCutForAugerElectrons(G4double cut)
G4double BarkasTerm(const G4Material *material, G4double kineticEnergy) const
G4double ElectronicLossFluctuation(const G4DynamicParticle *particle, const G4MaterialCutsCouple *material, G4double meanLoss, G4double step) const
const G4ThreeVector & GetMomentumDirection() const
virtual G4double HighEnergyLimit(const G4ParticleDefinition *aParticle, const G4Material *material) const =0
G4double antiprotonHighEnergy
G4double ProtonParametrisedDEDX(const G4MaterialCutsCouple *couple, G4double kineticEnergy) const
G4double GetCharge() const
void PutValue(size_t index, G4double theValue)
void SetElectronicStoppingPowerModel(const G4ParticleDefinition *aParticle, const G4String &dedxTable)
static G4Proton * Proton()
G4int SelectRandomAtom(const G4Material *material, G4double e) const
const G4String & GetParticleType() const
G4double minElectronEnergy
static G4ThreadLocal G4double dRoverRange
void ActivateAugerElectronProduction(G4bool val)
void LoadShellData(const G4String &dataFile)
static G4ThreadLocal G4int CounterOfpProcess
G4PhysicsTable * theMeanFreePathTable
const G4String & GetProcessName() const
const G4double * GetAtomicNumDensityVector() const
void SetKineticEnergy(G4double aEnergy)
static void BuildDEDXTable(const G4ParticleDefinition &aParticleType)
static G4ThreadLocal G4double HighestKineticEnergy
G4double GetTotNbOfAtomsPerVolume() const
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static G4ThreadLocal G4int CounterOfpbarProcess
G4int SelectRandomShell(G4int Z, G4double e) const
G4double GetLogEnergy1fluct() const
G4VLowEnergyModel * protonModel
G4double GetRateionexcfluct() const
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
G4bool CutsWhereModified()
static G4ThreadLocal G4bool EnlossFlucFlag
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
virtual G4double TheValue(const G4DynamicParticle *particle, const G4Material *material)=0
G4VLowEnergyModel * antiprotonModel
G4double GetPDGSpin() const
G4double GetMeanExcitationEnergy() const
G4PhysicsTable * theLossTable
static G4Electron * Electron()
G4double GetF2fluct() const
void SetAntiProtonElectronicStoppingPowerModel(const G4String &dedxTable)
size_t GetNumberOfElements() const
std::vector< G4DynamicParticle * > * GenerateParticles(G4int Z, G4int shellId)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeTrackStatus(G4TrackStatus status)
static void Register(const G4ParticleDefinition *p, const G4PhysicsTable *tDEDX, const G4PhysicsTable *tRange, const G4PhysicsTable *tInverseRange, const G4PhysicsTable *tLabTime, const G4PhysicsTable *tProperTime, G4double lowestKineticEnergy, G4double highestKineticEnergy, G4double massRatio, G4int NumberOfBins)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, enum G4ForceCondition *condition)
static const double eplus
G4double BlochTerm(const G4Material *material, G4double kineticEnergy, G4double cSquare) const
G4double GetPDGCharge() const
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
G4double bindingEnergy(G4int A, G4int Z)
static G4ThreadLocal G4PhysicsTable * theInverseRangepTable
static G4AtomicTransitionManager * Instance()
G4ProductionCuts * GetProductionCuts() const
G4double ComputeDEDX(const G4ParticleDefinition *aParticle, const G4MaterialCutsCouple *couple, G4double kineticEnergy)
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &Step)
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
G4VParticleChange * AlongStepDoIt(const G4Track &trackData, const G4Step &stepData)
G4double GetF1fluct() const
static G4ThreadLocal G4PhysicsTable * theDEDXpTable
G4ProcessVector * GetProcessList() const
G4VLowEnergyModel * betheBlochModel
const G4Material * GetMaterial() const
G4double GetEnergy1fluct() const
void PrintInfoDefinition() const