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 <<
" Tmin(MeV)= " << fLowestKineticEnergy/
MeV
82 <<
" Tmax(GeV)= " << fHighestKineticEnergy/
GeV
91 size_t n = fPAIxscBank.size();
93 for(
size_t i=0; i<
n; ++i) {
95 delete fPAIxscBank[i];
98 delete fPAIdEdxBank[i];
110 fSandia.Initialize(const_cast<G4Material*>(mat));
117 fHighestKineticEnergy,
122 fHighestKineticEnergy,
127 fHighestKineticEnergy,
131 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
136 for (
G4int i = 0; i <= fTotBin; ++i) {
138 G4double kinEnergy = fParticleEnergyVector->Energy(i);
143 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
145 fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
150 G4int n = fPAIySection.GetSplineSize();
154 for(
G4int k = 0; k <
n; k++ )
156 G4double t = fPAIySection.GetSplineEnergy(k+1);
158 t*fPAIySection.GetIntegralPAIySection(k+1));
159 dEdxVector->
PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
163 G4double ionloss = fPAIySection.GetMeanEnergyLoss();
165 if(ionloss < 0.0) ionloss = 0.0;
167 dEdxMeanVector->
PutValue(i,ionloss);
172 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
173 dNdxCutVector->
PutValue(i, dNdxCut);
174 dEdxCutVector->
PutValue(i, dEdxCut);
176 PAItransferTable->
insertAt(i,transferVector);
177 PAIdEdxTable->
insertAt(i,dEdxVector);
181 fPAIxscBank.push_back(PAItransferTable);
182 fPAIdEdxBank.push_back(PAIdEdxTable);
183 fdEdxTable.push_back(dEdxMeanVector);
185 fdNdxCutTable.push_back(dNdxCutVector);
186 fdEdxCutTable.push_back(dEdxCutVector);
196 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
197 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
200 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
201 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
206 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
207 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
209 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
210 G4double E1 = fParticleEnergyVector->Energy(iPlace);
211 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
220 if( dEdx < 0.) { dEdx = 0.; }
234 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
235 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
238 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
239 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
246 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
248 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
250 cross = (cross2-cross1);
253 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
254 - (*table)(iPlace+1)->Value(tmax)/tmax;
256 G4double E1 = fParticleEnergyVector->Energy(iPlace);
257 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
265 if( cross < 0.0) { cross = 0.0; }
280 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
281 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
284 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
285 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
292 dNdxCut1 = (*vcut)[iPlace];
296 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
306 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
307 dNdxCut2 = (*vcut)[iPlace+1];
309 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
310 E1 = fParticleEnergyVector->Energy(iPlace);
311 E2 = fParticleEnergyVector->Energy(iPlace+1);
313 W1 = (E2 - scaledTkin)*W;
314 W2 = (scaledTkin - E1)*W;
315 meanNumber = W1*meanN1 + W2*meanN2;
319 if(meanNumber < 0.0) {
return 0.0; }
325 if(0 == numOfCollisions) {
return 0.0; }
327 for(
G4int i=0; i< numOfCollisions; ++i) {
329 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
330 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
333 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
334 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
335 omega = omega*W1 + omega2*W2;
340 if(loss > kinEnergy) {
break; }
345 if(loss > kinEnergy) { loss = kinEnergy; }
346 else if(loss < 0.) { loss = 0.; }
362 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
370 if(scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
372 position = (*cutv)[nPlace]*rand;
373 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
375 else if(scaledTkin <= fParticleEnergyVector->Energy(0))
377 position = (*cutv)[0]*rand;
378 transfer = GetEnergyTransfer(coupleIndex, 0, position);
382 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
384 dNdxCut1 = (*cutv)[iPlace];
385 dNdxCut2 = (*cutv)[iPlace+1];
387 E1 = fParticleEnergyVector->Energy(iPlace);
388 E2 = fParticleEnergyVector->Energy(iPlace+1);
390 W1 = (E2 - scaledTkin)*W;
391 W2 = (scaledTkin - E1)*W;
396 position = dNdxCut1*rand;
397 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
399 position = dNdxCut2*rand;
400 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
402 transfer = tr1*W1 + tr2*W2;
405 if(transfer < 0.0 ) { transfer = 0.0; }
414 G4double G4PAIModelData::GetEnergyTransfer(
G4int coupleIndex,
419 if(position*v->
Energy(0) >= (*v)[0]) {
return v->
Energy(0); }
426 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
428 y2 = (*v)[iTransfer]/
x2;
429 if(position >=
y2) {
break; }
433 y1 = (*v)[iTransfer-1]/
x1;
443 const G4int nbins = 5;
446 for(
G4int i=1; i<=nbins; ++i) {
449 if(position >=
y2) {
break; }
460 return energyTransfer;
G4long G4Poisson(G4double mean)
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
size_t GetVectorLength() const
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4GLOB_DLL std::ostream G4cout
void PutValue(size_t index, G4double theValue)
const XML_Char XML_Content * model
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 *)
void insertAt(size_t, G4PhysicsVector *)
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double ComputeMaxEnergy(G4double scaledEnergy)
const G4Material * GetMaterial() const