60 const G4int nPerDecade = 10;
64 fPAIySection.SetVerbose(ver);
66 fLowestKineticEnergy =
std::max(tmin, lowestTkin);
67 fHighestKineticEnergy = tmax;
68 if(tmax < 10*fLowestKineticEnergy) {
69 fHighestKineticEnergy = 10*fLowestKineticEnergy;
70 }
else if(tmax > highestTkin) {
71 fHighestKineticEnergy =
std::max(highestTkin, 10*fLowestKineticEnergy);
73 fTotBin = (
G4int)(nPerDecade*
74 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
77 fHighestKineticEnergy,
80 G4cout <<
"### G4PAIModelData: Nbins= " << fTotBin
81 <<
" Tlowest(MeV)= " << fLowestKineticEnergy/
MeV
82 <<
" Tmin(keV)= " << tmin/
keV
83 <<
" Tmax(GeV)= " << fHighestKineticEnergy/
GeV
92 size_t n = fPAIxscBank.size();
94 for(
size_t i=0; i<
n; ++i) {
96 fPAIxscBank[i]->clearAndDestroy();
97 delete fPAIxscBank[i];
100 fPAIdEdxBank[i]->clearAndDestroy();
101 delete fPAIdEdxBank[i];
103 delete fdEdxTable[i];
104 delete fdNdxCutTable[i];
105 delete fdEdxCutTable[i];
108 delete fParticleEnergyVector;
117 fSandia.Initialize(const_cast<G4Material*>(mat));
124 fHighestKineticEnergy,
129 fHighestKineticEnergy,
134 fHighestKineticEnergy,
138 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
143 for (
G4int i = 0; i <= fTotBin; ++i) {
145 G4double kinEnergy = fParticleEnergyVector->Energy(i);
147 G4double tau = kinEnergy/proton_mass_c2;
150 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
152 fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
157 G4int n = fPAIySection.GetSplineSize();
159 for(
G4int k = 0; k <
n; ++k) {
160 if(fPAIySection.GetIntegralPAIySection(k+1) <= 0.0) {
173 for(
G4int k = kmin; k <
n; ++k)
175 G4double t = fPAIySection.GetSplineEnergy(k+1);
176 tr = fPAIySection.GetIntegralPAIySection(k+1);
181 transferVector->
PutValue(k , t, t*tr);
182 dEdxVector->
PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
186 G4double ionloss = fPAIySection.GetMeanEnergyLoss();
188 if(ionloss < 0.0) ionloss = 0.0;
190 dEdxMeanVector->
PutValue(i,ionloss);
195 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
196 if(dEdxCut < 0.0) { dEdxCut = 0.0; }
197 dNdxCutVector->
PutValue(i, dNdxCut);
198 dEdxCutVector->
PutValue(i, dEdxCut);
200 PAItransferTable->
insertAt(i,transferVector);
201 PAIdEdxTable->
insertAt(i,dEdxVector);
209 fPAIxscBank.push_back(PAItransferTable);
210 fPAIdEdxBank.push_back(PAIdEdxTable);
219 fdEdxTable.push_back(dEdxMeanVector);
220 fdNdxCutTable.push_back(dNdxCutVector);
221 fdEdxCutTable.push_back(dEdxCutVector);
231 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
232 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
235 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
236 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
241 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
242 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
244 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
245 G4double E1 = fParticleEnergyVector->Energy(iPlace);
246 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
255 if( dEdx < 0.) { dEdx = 0.; }
269 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
270 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
273 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
274 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
281 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
283 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
285 cross = (cross2-cross1);
288 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
289 - (*table)(iPlace+1)->Value(tmax)/tmax;
291 G4double E1 = fParticleEnergyVector->Energy(iPlace);
292 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
300 if( cross < 0.0) { cross = 0.0; }
315 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
316 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
319 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
320 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
327 dNdxCut1 = (*vcut)[iPlace];
331 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
341 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
342 dNdxCut2 = (*vcut)[iPlace+1];
344 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
345 E1 = fParticleEnergyVector->Energy(iPlace);
346 E2 = fParticleEnergyVector->Energy(iPlace+1);
348 W1 = (E2 - scaledTkin)*W;
349 W2 = (scaledTkin - E1)*W;
350 meanNumber = W1*meanN1 + W2*meanN2;
354 if(meanNumber < 0.0) {
return 0.0; }
360 if(0 == numOfCollisions) {
return 0.0; }
362 for(
G4int i=0; i< numOfCollisions; ++i) {
364 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
365 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
368 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
369 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
370 omega = omega*W1 + omega2*W2;
375 if(loss > kinEnergy) {
break; }
380 if(loss > kinEnergy) { loss = kinEnergy; }
381 else if(loss < 0.) { loss = 0.; }
399 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
400 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
403 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
404 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
415 G4double dNdxCut1 = (*vcut)[iPlace];
424 position = dNdxCutM + (dNdxCut1 - dNdxCutM)*rand;
425 transfer = GetEnergyTransfer(coupleIndex, iPlace, position);
433 G4double dNdxCut2 = (*vcut)[iPlace+1];
441 E1 = fParticleEnergyVector->Energy(iPlace);
442 E2 = fParticleEnergyVector->Energy(iPlace+1);
444 W1 = (E2 - scaledTkin)*W;
445 W2 = (scaledTkin - E1)*W;
451 position = dNdxCutM + (dNdxCut1 - dNdxCutM)*rand;
454 tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
459 position = dNdxCutM2 + (dNdxCut2 - dNdxCutM2)*rand;
460 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
464 transfer = tr1*W1 + tr2*W2;
468 if(transfer < 0.0 ) { transfer = 0.0; }
482 if(position*v->
Energy(0) >= (*v)[0]) {
return v->
Energy(0); }
487 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
490 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
491 x2 = v->
Energy(iTransfer);
492 y2 = (*v)[iTransfer]/x2;
493 if(position >= y2) {
break; }
494 if(iTransfer == iTransferMax) {
return v->
GetMaxEnergy(); }
497 x1 = v->
Energy(iTransfer-1);
498 y1 = (*v)[iTransfer-1]/x1;
509 const G4int nbins = 5;
512 for(
G4int i=1; i<=nbins; ++i) {
514 y2 = v->
Value(x2)/x2;
515 if(position >= y2) {
break; }
522 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
528 return energyTransfer;
G4long G4Poisson(G4double mean)
G4double GetMaxEnergy() const
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
G4double GetEnergyTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
size_t GetVectorLength() const
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin, G4double scaledTmax) const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4GLOB_DLL std::ostream G4cout
void PutValue(size_t index, G4double theValue)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIModel *)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void insertAt(size_t, G4PhysicsVector *)
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double ComputeMaxEnergy(G4double scaledEnergy)
const G4Material * GetMaterial() const