59 const G4int nPerDecade = 10;
65 fLowestKineticEnergy = std::max(tmin, lowestTkin);
66 fHighestKineticEnergy = tmax;
67 if(tmax < 10*fLowestKineticEnergy) {
68 fHighestKineticEnergy = 10*fLowestKineticEnergy;
69 }
else if(tmax > highestTkin) {
70 fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
72 fTotBin = (
G4int)(nPerDecade*
73 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
76 fHighestKineticEnergy,
79 G4cout <<
"### G4PAIModelData: Nbins= " << fTotBin
80 <<
" Tlowest(keV)= " << lowestTkin/keV
81 <<
" Tmin(keV)= " << fLowestKineticEnergy/keV
82 <<
" Tmax(GeV)= " << fHighestKineticEnergy/GeV
91 std::size_t n = fPAIxscBank.size();
93 for(std::size_t i=0; i<n; ++i) {
95 fPAIxscBank[i]->clearAndDestroy();
96 delete fPAIxscBank[i];
99 fPAIdEdxBank[i]->clearAndDestroy();
100 delete fPAIdEdxBank[i];
102 delete fdEdxTable[i];
105 delete fParticleEnergyVector;
118 auto dEdxMeanVector =
120 fHighestKineticEnergy,
126 static const G4double deltaLow = 100.*eV;
128 for (
G4int i = 0; i <= fTotBin; ++i) {
132 G4double tau = kinEnergy/proton_mass_c2;
135 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
137 fPAIySection.
Initialize(mat, Tmax, bg2, &fSandia);
144 for(
G4int k = 0; k < n; ++k) {
158 for(
G4int k = kmin; k < n; ++k)
166 transferVector->PutValue(k, t, t*tr);
175 dEdxMeanVector->PutValue(i,ionloss);
177 PAItransferTable->insertAt(i,transferVector);
178 PAIdEdxTable->insertAt(i,dEdxVector);
181 fPAIxscBank.push_back(PAItransferTable);
182 fPAIdEdxBank.push_back(PAIdEdxTable);
185 fdEdxTable.push_back(dEdxMeanVector);
195 std::size_t iPlace(0);
196 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin, iPlace);
204 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
205 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
210 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
213 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
225 dEdx = std::max(dEdx, 0.);
239 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
243 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
244 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
251 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
253 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
255 cross = (cross2-cross1);
258 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
259 - (*table)(iPlace+1)->Value(tmax)/tmax;
270 cross = std::max(cross, 0.0);
285 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
289 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
290 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
307 meanN11 = (*v1)[0]/e1;
308 meanN12 = v1->
Value(e2)/e2;
309 meanNumber = (meanN11 - meanN12)*stepFactor;
317 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
322 meanN21 = (*v2)[0]/e1;
323 meanN22 = v2->
Value(e2)/e2;
327 W1 = (E2 - scaledTkin)*
W;
328 W2 = (scaledTkin - E1)*
W;
330 meanNumber += (meanN21 - meanN22)*stepFactor*W2;
334 if(meanNumber < 0.0) {
return loss; }
339 if(0 == numOfCollisions) {
return loss; }
342 for(
G4int i=0; i< numOfCollisions; ++i) {
344 position = meanN12 + (meanN11 - meanN12)*rand;
345 omega = GetEnergyTransfer(coupleIndex, iPlace,
position);
348 position = meanN22 + (meanN21 - meanN22)*rand;
349 omega2 = GetEnergyTransfer(coupleIndex, iPlace+1,
position);
356 if(loss > kinEnergy) {
break; }
360 if(loss > kinEnergy) { loss = kinEnergy; }
361 else if(loss < 0.) { loss = 0.; }
381 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
384 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
385 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
393 if(emax < emin) {
return transfer; }
404 transfer = GetEnergyTransfer(coupleIndex, iPlace,
position);
412 emin = std::max(tmin, v2->
Energy(0));
415 dNdx1 = v2->
Value(emin)/emin;
416 dNdx2 = v2->
Value(emax)/emax;
430 position = dNdx2 + (dNdx1 - dNdx2)*rand;
441 transfer = std::max(transfer, 0.0);
459 std::size_t iTransfer;
460 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
463 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
464 x2 = v->
Energy(iTransfer);
465 y2 = (*v)[iTransfer]/x2;
467 if(iTransfer == iTransferMax) {
return v->
GetMaxEnergy(); }
470 x1 = v->
Energy(iTransfer-1);
471 y1 = (*v)[iTransfer-1]/x1;
483 const G4int nbins = 5;
486 for(
G4int i=1; i<=nbins; ++i) {
488 y2 = v->
Value(x2)/x2;
498 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
504 return energyTransfer;
G4long G4Poisson(G4double mean)
G4GLOB_DLL std::ostream G4cout
const G4Material * GetMaterial() const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
void Initialise(const G4MaterialCutsCouple *, G4PAIModel *)
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin, G4double tmin, G4double tmax) const
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double tmax, G4double stepFactor) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double ComputeMaxEnergy(G4double scaledEnergy)
G4double GetSplineEnergy(G4int i) const
G4double GetIntegralPAIdEdx(G4int i) const
G4double GetMeanEnergyLoss() const
G4double GetIntegralPAIySection(G4int i) const
G4int GetSplineSize() const
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double GetMaxEnergy() const
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
std::size_t GetVectorLength() const
std::size_t FindBin(const G4double energy, std::size_t idx) const
void Initialize(const G4Material *)
G4double GetSandiaMatTablePAI(G4int, G4int) const