77 G4Exception(
"G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()",
"InvalidCall",
138 delete theRangeTable; }
143 for (
G4int J=0; J<numOfCouples; J++)
147 highestKineticEnergy,TotBin);
150 theRangeTable->
insert(aVector);
152 return theRangeTable ;
168 G4double energy1 = lowestKineticEnergy;
175 for (
G4int j=1; j<TotBin; j++) {
178 G4double de = (energy2 - energy1) * del ;
181 for (
G4int i=1; i<
n; i++) {
184 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
199 G4double dtau,Value,taui,ti,lossi,ci;
204 for (
G4int i=0; i<=nbin; i++)
208 lossi = physicsVector->
GetValue(ti,isOut);
230 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
236 for (
G4int i=0; i<=nbin; i++)
238 ui = ltaulow+dltau*i;
241 lossi = physicsVector->
GetValue(ti,isOut);
251 Value += ci*taui/lossi;
271 delete theLabTimeTable; }
275 for (
G4int J=0; J<numOfCouples; J++)
280 highestKineticEnergy,TotBin);
283 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
284 theLabTimeTable->
insert(aVector);
288 return theLabTimeTable ;
302 if(theProperTimeTable)
304 delete theProperTimeTable; }
308 for (
G4int J=0; J<numOfCouples; J++)
313 highestKineticEnergy,TotBin);
316 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
317 theProperTimeTable->
insert(aVector);
321 return theProperTimeTable ;
336 G4double tlim=5.*
keV,parlowen=0.4,ppar=0.5-parlowen ;
337 G4double losslim,clim,taulim,timelim,
338 LowEdgeEnergy,tau,Value ;
343 losslim = physicsVector->
GetValue(tlim,isOut);
357 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
369 }
while (tau<=taulim) ;
371 for (
G4int j=i; j<TotBin; j++)
395 G4double tlim=5.*
keV,parlowen=0.4,ppar=0.5-parlowen ;
396 G4double losslim,clim,taulim,timelim,
397 LowEdgeEnergy,tau,Value ;
403 losslim = physicsVector->
GetValue(tlim,isOut);
417 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
429 }
while (tau<=taulim) ;
431 for (
G4int j=i; j<TotBin; j++)
450 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
456 for (
G4int i=0; i<=nbin; i++)
458 ui = ltaulow+dltau*i;
461 lossi = physicsVector->
GetValue(ti,isOut);
483 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
489 for (
G4int i=0; i<=nbin; i++)
491 ui = ltaulow+dltau*i;
494 lossi = physicsVector->
GetValue(ti,isOut);
526 if(theInverseRangeTable)
528 delete theInverseRangeTable; }
532 for (
G4int i=0; i<numOfCouples; i++)
542 rhigh *= std::exp(std::log(rhigh/rlow)/((
G4double)(nbins-1)));
554 for (
size_t j=1; j<nbins; j++) {
558 for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
561 if(range2 >= range || ihigh == nbins-1) {
569 G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
573 theInverseRangeTable->
insert(v);
576 return theInverseRangeTable ;
591 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
592 G4double Tbin = lowestKineticEnergy/RTable ;
594 G4int binnumber = -1 ;
598 for(
G4int i=0; i<TotBin; i++)
601 if( rangebin < LowEdgeRange )
607 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
609 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
613 KineticEnergy = lowestKineticEnergy ;
614 else if(binnumber == TotBin-1)
615 KineticEnergy = highestKineticEnergy ;
618 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
619 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
620 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
622 KineticEnergy = (LowEdgeRange -
C )/B ;
625 discr = B*B - 4.*A*(C-LowEdgeRange);
626 discr = discr>0. ? std::sqrt(discr) : 0.;
627 KineticEnergy = 0.5*(discr-
B)/A ;
631 aVector->
PutValue(i,KineticEnergy) ;
647 if(theRangeCoeffATable)
649 delete theRangeCoeffATable; }
652 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
656 G4double w1 = RTable/
w , w2 = -RTable*R1/
w , w3 = R2/
w ;
657 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
661 for (
G4int J=0; J<numOfCouples; J++)
663 G4int binmax=TotBin ;
666 Ti = lowestKineticEnergy ;
669 for (
G4int i=0; i<TotBin; i++)
671 Ri = rangeVector->
GetValue(Ti,isOut) ;
677 Rim = rangeVector->
GetValue(Tim,isOut);
684 Rip = rangeVector->
GetValue(Tip,isOut);
686 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
692 theRangeCoeffATable->
insert(aVector);
694 return theRangeCoeffATable ;
709 if(theRangeCoeffBTable)
711 delete theRangeCoeffBTable; }
714 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
718 G4double w1 = -R1/
w , w2 = R1*(R2+1.)/
w , w3 = -R2*R1/
w ;
719 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
723 for (
G4int J=0; J<numOfCouples; J++)
725 G4int binmax=TotBin ;
728 Ti = lowestKineticEnergy ;
731 for (
G4int i=0; i<TotBin; i++)
733 Ri = rangeVector->
GetValue(Ti,isOut) ;
739 Rim = rangeVector->
GetValue(Tim,isOut);
746 Rip = rangeVector->
GetValue(Tip,isOut);
748 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
753 theRangeCoeffBTable->
insert(aVector);
755 return theRangeCoeffBTable ;
770 if(theRangeCoeffCTable)
772 delete theRangeCoeffCTable; }
775 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
779 G4double w1 = 1./
w , w2 = -RTable*R1/
w , w3 = RTable*R2/
w ;
780 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
784 for (
G4int J=0; J<numOfCouples; J++)
786 G4int binmax=TotBin ;
789 Ti = lowestKineticEnergy ;
792 for (
G4int i=0; i<TotBin; i++)
794 Ri = rangeVector->
GetValue(Ti,isOut) ;
800 Rim = rangeVector->
GetValue(Tim,isOut);
807 Rip = rangeVector->
GetValue(Tip,isOut);
809 Value = w1*Rip + w2*Ri + w3*Rim ;
814 theRangeCoeffCTable->
insert(aVector);
816 return theRangeCoeffCTable ;
829 static const G4double probLim = 0.01 ;
830 static const G4double sumaLim = -std::log(probLim) ;
853 beta2,suma,e0,loss,lossc,
w;
857 G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
863 if(MeanLoss < minLoss)
return MeanLoss ;
871 ->GetEnergyCutsVector(1)))[
imat];
878 if(Tm > threshold) Tm = threshold;
879 beta2 = tau2/(tau1*tau1);
882 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*
ipotFluct)
885 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
886 factor*electronDensity/beta2) ;
889 }
while (loss < 0. || loss > 2.0*MeanLoss);
930 Tm = Tm-ipotFluct+e0 ;
941 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
1009 siga=std::sqrt(a3) ;
1030 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1031 ea = na*ipotFluct*alfa1;
1032 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));