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);
176 if(ionloss < 0.0) ionloss = 0.0;
178 dEdxMeanVector->PutValue(i,ionloss);
180 PAItransferTable->insertAt(i,transferVector);
181 PAIdEdxTable->insertAt(i,dEdxVector);
184 fPAIxscBank.push_back(PAItransferTable);
185 fPAIdEdxBank.push_back(PAIdEdxTable);
188 fdEdxTable.push_back(dEdxMeanVector);
198 std::size_t iPlace(0);
199 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin, iPlace);
207 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
208 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
213 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
216 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
228 dEdx = std::max(dEdx, 0.);
242 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
246 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
247 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
254 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
256 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
258 cross = (cross2-cross1);
261 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
262 - (*table)(iPlace+1)->Value(tmax)/tmax;
273 cross = std::max(cross, 0.0);
288 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
292 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
293 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
310 meanN11 = (*v1)[0]/e1;
311 meanN12 = v1->
Value(e2)/e2;
312 meanNumber = (meanN11 - meanN12)*stepFactor;
320 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
325 meanN21 = (*v2)[0]/e1;
326 meanN22 = v2->
Value(e2)/e2;
330 W1 = (E2 - scaledTkin)*
W;
331 W2 = (scaledTkin - E1)*
W;
333 meanNumber += (meanN21 - meanN22)*stepFactor*W2;
337 if(meanNumber < 0.0) {
return loss; }
342 if(0 == numOfCollisions) {
return loss; }
345 for(
G4int i=0; i< numOfCollisions; ++i) {
347 position = meanN12 + (meanN11 - meanN12)*rand;
348 omega = GetEnergyTransfer(coupleIndex, iPlace,
position);
351 position = meanN22 + (meanN21 - meanN22)*rand;
352 omega2 = GetEnergyTransfer(coupleIndex, iPlace+1,
position);
359 if(loss > kinEnergy) {
break; }
363 if(loss > kinEnergy) { loss = kinEnergy; }
364 else if(loss < 0.) { loss = 0.; }
384 std::size_t iPlace = fParticleEnergyVector->
FindBin(scaledTkin, 0);
387 if(scaledTkin >= fParticleEnergyVector->
Energy(nPlace)) { iPlace = nPlace; }
388 else if(scaledTkin > fParticleEnergyVector->
Energy(0)) {
396 if(emax < emin) {
return transfer; }
407 transfer = GetEnergyTransfer(coupleIndex, iPlace,
position);
415 emin = std::max(tmin, v2->
Energy(0));
418 dNdx1 = v2->
Value(emin)/emin;
419 dNdx2 = v2->
Value(emax)/emax;
433 position = dNdx2 + (dNdx1 - dNdx2)*rand;
444 transfer = std::max(transfer, 0.0);
462 std::size_t iTransfer;
463 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
466 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
467 x2 = v->
Energy(iTransfer);
468 y2 = (*v)[iTransfer]/x2;
470 if(iTransfer == iTransferMax) {
return v->
GetMaxEnergy(); }
473 x1 = v->
Energy(iTransfer-1);
474 y1 = (*v)[iTransfer-1]/x1;
486 const G4int nbins = 5;
489 for(
G4int i=1; i<=nbins; ++i) {
491 y2 = v->
Value(x2)/x2;
501 energyTransfer = (y2 - y1)*x1*x2/(
position*(x1 - x2) - y1*x1 + y2*x2);
507 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