Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4PAIModelData Class Reference

#include <G4PAIModelData.hh>

Public Member Functions

 G4PAIModelData (G4double tmin, G4double tmax, G4int verbose)
 
 ~G4PAIModelData ()
 
void Initialise (const G4MaterialCutsCouple *, G4PAIModel *)
 
G4double DEDXPerVolume (G4int coupleIndex, G4double scaledTkin, G4double cut) const
 
G4double CrossSectionPerVolume (G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
 
G4double SampleAlongStepTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double tmax, G4double stepFactor) const
 
G4double SamplePostStepTransfer (G4int coupleIndex, G4double scaledTkin, G4double tmin, G4double tmax) const
 
G4PAIModelData & operator= (const G4PAIModelData &right)=delete
 
 G4PAIModelData (const G4PAIModelData &)=delete
 

Detailed Description

Definition at line 67 of file G4PAIModelData.hh.

Constructor & Destructor Documentation

◆ G4PAIModelData() [1/2]

G4PAIModelData::G4PAIModelData ( G4double  tmin,
G4double  tmax,
G4int  verbose 
)
explicit

Definition at line 57 of file G4PAIModelData.cc.

58{
59 const G4int nPerDecade = 10;
60 const G4double lowestTkin = 50*keV;
61 const G4double highestTkin = 10*TeV;
62
63 fPAIySection.SetVerbose(ver);
64
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);
71 }
72 fTotBin = (G4int)(nPerDecade*
73 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
74
75 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
76 fHighestKineticEnergy,
77 fTotBin);
78 if(0 < ver) {
79 G4cout << "### G4PAIModelData: Nbins= " << fTotBin
80 << " Tlowest(keV)= " << lowestTkin/keV
81 << " Tmin(keV)= " << fLowestKineticEnergy/keV
82 << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
83 << G4endl;
84 }
85}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIySection::SetVerbose
void SetVerbose(G4int v)
Definition: G4PAIySection.hh:139

◆ ~G4PAIModelData()

G4PAIModelData::~G4PAIModelData ( )

Definition at line 89 of file G4PAIModelData.cc.

90{
91 std::size_t n = fPAIxscBank.size();
92 if(0 < n) {
93 for(std::size_t i=0; i<n; ++i) {
94 if(fPAIxscBank[i]) {
95 fPAIxscBank[i]->clearAndDestroy();
96 delete fPAIxscBank[i];
97 }
98 if(fPAIdEdxBank[i]) {
99 fPAIdEdxBank[i]->clearAndDestroy();
100 delete fPAIdEdxBank[i];
101 }
102 delete fdEdxTable[i];
103 }
104 }
105 delete fParticleEnergyVector;
106}

◆ G4PAIModelData() [2/2]

G4PAIModelData::G4PAIModelData ( const G4PAIModelData &  )
delete

Member Function Documentation

◆ CrossSectionPerVolume()

G4double G4PAIModelData::CrossSectionPerVolume ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  tcut,
G4double  tmax 
) const

Definition at line 235 of file G4PAIModelData.cc.

238{
239 G4double cross, cross1, cross2;
240
241 // iPlace is in interval from 0 to (N-1)
242 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
243 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
244
245 G4bool one = true;
246 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
247 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
248 one = false;
249 }
250 G4PhysicsTable* table = fPAIxscBank[coupleIndex];
251
252 //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
253 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
254 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
255 //G4cout<<"cross1 = "<<cross1<<G4endl;
256 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
257 //G4cout<<"cross2 = "<<cross2<<G4endl;
258 cross = (cross2-cross1);
259 //G4cout<<"cross = "<<cross<<G4endl;
260 if(!one) {
261 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
262 - (*table)(iPlace+1)->Value(tmax)/tmax;
263
264 G4double E1 = fParticleEnergyVector->Energy(iPlace);
265 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
266 G4double W = 1.0/(E2 - E1);
267 G4double W1 = (E2 - scaledTkin)*W;
268 G4double W2 = (scaledTkin - E1)*W;
269 cross *= W1;
270 cross += W2*cross2;
271 }
272
273 cross = std::max(cross, 0.0);
274 return cross;
275}
G4bool
bool G4bool
Definition: G4Types.hh:86
G4PhysicsVector::Energy
G4double Energy(const std::size_t index) const
G4PhysicsVector::GetVectorLength
std::size_t GetVectorLength() const
G4PhysicsVector::FindBin
std::size_t FindBin(const G4double energy, std::size_t idx) const
Definition: G4PhysicsVector.cc:173
W
#define W
Definition: crc32.c:84

Referenced by G4PAIModel::CrossSectionPerVolume().

◆ DEDXPerVolume()

G4double G4PAIModelData::DEDXPerVolume ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  cut 
) const

Definition at line 193 of file G4PAIModelData.cc.

195{
196 // VI: iPlace is the low edge index of the bin
197 // iPlace is in interval from 0 to (N-1)
198 std::size_t iPlace(0);
199 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin, iPlace);
200 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
201 /*
202 G4cout << "G4PAIModelData::DEDXPerVolume: coupleIdx= " << coupleIndex
203 << " Tscaled= " << scaledTkin << " cut= " << cut
204 << " iPlace= " << iPlace << " nPlace= " << nPlace << G4endl;
205 */
206 G4bool one = true;
207 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
208 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
209 one = false;
210 }
211
212 // VI: apply interpolation of the vector
213 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
214 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
215 if(!one) {
216 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
217 G4double E1 = fParticleEnergyVector->Energy(iPlace);
218 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
219 G4double W = 1.0/(E2 - E1);
220 G4double W1 = (E2 - scaledTkin)*W;
221 G4double W2 = (scaledTkin - E1)*W;
222 del *= W1;
223 del += W2*del2;
224 }
225 dEdx -= del;
226 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
227
228 dEdx = std::max(dEdx, 0.);
229 return dEdx;
230}

Referenced by G4PAIModel::ComputeDEDXPerVolume().

◆ Initialise()

void G4PAIModelData::Initialise ( const G4MaterialCutsCouple *  couple,
G4PAIModel *  model 
)

Definition at line 110 of file G4PAIModelData.cc.

112{
113 const G4Material* mat = couple->GetMaterial();
114 fSandia.Initialize(const_cast<G4Material*>(mat));
115
116 auto PAItransferTable = new G4PhysicsTable(fTotBin+1);
117 auto PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
118 auto dEdxMeanVector =
119 new G4PhysicsLogVector(fLowestKineticEnergy,
120 fHighestKineticEnergy,
121 fTotBin);
122 // low energy Sandia interval
123 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
124
125 // energy safety
126 static const G4double deltaLow = 100.*eV;
127
128 for (G4int i = 0; i <= fTotBin; ++i) {
129
130 G4double kinEnergy = fParticleEnergyVector->Energy(i);
131 G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
132 G4double tau = kinEnergy/proton_mass_c2;
133 G4double bg2 = tau*( tau + 2. );
134
135 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
136
137 fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
138
139 //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV
140 // << " E(MeV)= " << kinEnergy/MeV << G4endl;
141
142 G4int n = fPAIySection.GetSplineSize();
143 G4int kmin = 0;
144 for(G4int k = 0; k < n; ++k) {
145 if(fPAIySection.GetIntegralPAIySection(k+1) <= 0.0) {
146 kmin = k;
147 } else {
148 break;
149 }
150 }
151 n -= kmin;
152
153 auto transferVector = new G4PhysicsFreeVector(n);
154 auto dEdxVector = new G4PhysicsFreeVector(n);
155
156 //G4double tr0 = 0.0;
157 G4double tr = 0.0;
158 for(G4int k = kmin; k < n; ++k)
159 {
160 G4double t = fPAIySection.GetSplineEnergy(k+1);
161 tr = fPAIySection.GetIntegralPAIySection(k+1);
162 //if(tr >= tr0) { tr0 = tr; }
163 //else { G4cout << "G4PAIModelData::Initialise Warning: Ekin(MeV)= "
164 // << t/MeV << " IntegralTransfer= " << tr
165 // << " < " << tr0 << G4endl; }
166 transferVector->PutValue(k, t, t*tr);
167 dEdxVector->PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
168 }
169 //G4cout << "TransferVector:" << G4endl;
170 //G4cout << *transferVector << G4endl;
171 //G4cout << "DEDXVector:" << G4endl;
172 //G4cout << *dEdxVector << G4endl;
173
174 G4double ionloss = fPAIySection.GetMeanEnergyLoss();// total <dE/dx>
175
176 if(ionloss < 0.0) ionloss = 0.0;
177
178 dEdxMeanVector->PutValue(i,ionloss);
179
180 PAItransferTable->insertAt(i,transferVector);
181 PAIdEdxTable->insertAt(i,dEdxVector);
182
183 } // end of Tkin loop`
184 fPAIxscBank.push_back(PAItransferTable);
185 fPAIdEdxBank.push_back(PAIdEdxTable);
186 //G4cout << "dEdxMeanVector: " << G4endl;
187 //G4cout << *dEdxMeanVector << G4endl;
188 fdEdxTable.push_back(dEdxMeanVector);
189}
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4PAIModel::ComputeMaxEnergy
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIModel.hh:166
G4PAIySection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIySection.hh:203
G4PAIySection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIySection.hh:215
G4PAIySection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIySection.hh:124
G4PAIySection::GetIntegralPAIySection
G4double GetIntegralPAIySection(G4int i) const
Definition: G4PAIySection.hh:209
G4PAIySection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIySection.hh:114
G4PAIySection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIySection.cc:118
G4SandiaTable::Initialize
void Initialize(const G4Material *)
Definition: G4SandiaTable.cc:715
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1066

Referenced by G4PAIModel::Initialise().

◆ operator=()

G4PAIModelData & G4PAIModelData::operator= ( const G4PAIModelData &  right)
delete

◆ SampleAlongStepTransfer()

G4double G4PAIModelData::SampleAlongStepTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  tmax,
G4double  stepFactor 
) const

Definition at line 279 of file G4PAIModelData.cc.

284{
285 //G4cout << "=== G4PAIModelData::SampleAlongStepTransfer" << G4endl;
286 G4double loss = 0.0;
287
288 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
289 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
290
291 G4bool one = true;
292 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
293 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
294 one = false;
295 }
296
297 G4double meanNumber = 0.0;
298 G4double meanN11 = 0.0;
299 G4double meanN12 = 0.0;
300 G4double meanN21 = 0.0;
301 G4double meanN22 = 0.0;
302
303 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
304 G4PhysicsVector* v2 = nullptr;
305
306 G4double e1 = v1->Energy(0);
307 G4double e2 = std::min(tmax, v1->GetMaxEnergy());
308
309 if(e2 >= e1) {
310 meanN11 = (*v1)[0]/e1;
311 meanN12 = v1->Value(e2)/e2;
312 meanNumber = (meanN11 - meanN12)*stepFactor;
313 }
314 //G4cout<<"iPlace = "<<iPlace<< " meanN11= " << meanN11
315 // << " meanN12= " << meanN12 << G4endl;
316
317 G4double W1 = 1.0;
318 G4double W2 = 0.0;
319 if(!one) {
320 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
321
322 e1 = v2->Energy(0);
323 e2 = std::min(tmax, v2->GetMaxEnergy());
324 if(e2 >= e1) {
325 meanN21 = (*v2)[0]/e1;
326 meanN22 = v2->Value(e2)/e2;
327 G4double E1 = fParticleEnergyVector->Energy(iPlace);
328 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
329 G4double W = 1.0/(E2 - E1);
330 W1 = (E2 - scaledTkin)*W;
331 W2 = (scaledTkin - E1)*W;
332 meanNumber *= W1;
333 meanNumber += (meanN21 - meanN22)*stepFactor*W2;
334 }
335 }
336
337 if(meanNumber < 0.0) { return loss; }
338 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
339
340 //G4cout << "meanNumber= " << meanNumber << " N= " << numOfCollisions << G4endl;
341
342 if(0 == numOfCollisions) { return loss; }
343
344 G4double position, omega, omega2;
345 for(G4int i=0; i< numOfCollisions; ++i) {
346 G4double rand = G4UniformRand();
347 position = meanN12 + (meanN11 - meanN12)*rand;
348 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
349 //G4cout << "omega(keV)= " << omega/keV << G4endl;
350 if(!one) {
351 position = meanN22 + (meanN21 - meanN22)*rand;
352 omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
353 omega *= W1;
354 omega += omega2*W2;
355 }
356 //G4cout << "omega(keV)= " << omega/keV << G4endl;
357
358 loss += omega;
359 if(loss > kinEnergy) { break; }
360 }
361
362 //G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
363 if(loss > kinEnergy) { loss = kinEnergy; }
364 else if(loss < 0.) { loss = 0.; }
365 return loss;
366}
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::Value
G4double Value(const G4double energy, std::size_t &lastidx) const

Referenced by G4PAIModel::SampleFluctuations().

◆ SamplePostStepTransfer()

G4double G4PAIModelData::SamplePostStepTransfer ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  tmin,
G4double  tmax 
) const

Definition at line 373 of file G4PAIModelData.cc.

377{
378 //G4cout<<"=== G4PAIModelData::SamplePostStepTransfer idx= "<< coupleIndex
379 // << " Tkin= " << scaledTkin << " Tmax= " << tmax << G4endl;
380 G4double transfer = 0.0;
381 G4double rand = G4UniformRand();
382
383 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
384 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
385
386 G4bool one = true;
387 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
388 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
389 one = false;
390 }
391 G4PhysicsTable* table = fPAIxscBank[coupleIndex];
392 G4PhysicsVector* v1 = (*table)[iPlace];
393
394 G4double emin = std::max(tmin, v1->Energy(0));
395 G4double emax = std::min(tmax, v1->GetMaxEnergy());
396 if(emax < emin) { return transfer; }
397
398 G4double dNdx1 = v1->Value(emin)/emin;
399 G4double dNdx2 = v1->Value(emax)/emax;
400 /*
401 G4cout << "iPlace= " << iPlace << " nPlace= " << nPlace
402 << " emin= " << emin << " emax= " << emax
403 << " dNdx1= " << dNdx1 << " dNdx2= " << dNdx2
404 << " one: " << one << G4endl;
405 */
406 G4double position = dNdx2 + (dNdx1 - dNdx2)*rand;
407 transfer = GetEnergyTransfer(coupleIndex, iPlace, position);
408
409 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
410 // << " position= " << position << G4endl;
411
412 if(!one) {
413
414 G4PhysicsVector* v2 = (*table)[iPlace+1];
415 emin = std::max(tmin, v2->Energy(0));
416 emax = std::min(tmax, v2->GetMaxEnergy());
417 if(emin <= emax) {
418 dNdx1 = v2->Value(emin)/emin;
419 dNdx2 = v2->Value(emax)/emax;
420
421 //G4cout << " emax2= " << emax
422 // << " dNdx2= " << dNdx2 << " dNdx1= " << dNdx1 << G4endl;
423
424 G4double E1 = fParticleEnergyVector->Energy(iPlace);
425 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
426 G4double W = 1.0/(E2 - E1);
427 G4double W1 = (E2 - scaledTkin)*W;
428 G4double W2 = (scaledTkin - E1)*W;
429
430 //G4cout<< "E1= " << E1 << " E2= " << E2 <<" iPlace= " << iPlace
431 // << " W1= " << W1 << " W2= " << W2 <<G4endl;
432
433 position = dNdx2 + (dNdx1 - dNdx2)*rand;
434 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
435
436 //G4cout<<"PAImodel PostStepTransfer1 = "<<tr2/keV<<" keV"
437 // << " position= " << position << G4endl;
438 transfer *= W1;
439 transfer += tr2*W2;
440 }
441 }
442 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
443 // << " position= " << position << G4endl;
444 transfer = std::max(transfer, 0.0);
445 return transfer;
446}

Referenced by G4PAIModel::SampleSecondaries().

The documentation for this class was generated from the following files: