Geant4 11.2.2
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:67
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 232 of file G4PAIModelData.cc.

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

Referenced by G4PAIModel::CrossSectionPerVolume().

◆ DEDXPerVolume()

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

Definition at line 190 of file G4PAIModelData.cc.

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

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 = std::max(fPAIySection.GetMeanEnergyLoss(), 0.0);// total <dE/dx>
175 dEdxMeanVector->PutValue(i,ionloss);
176
177 PAItransferTable->insertAt(i,transferVector);
178 PAIdEdxTable->insertAt(i,dEdxVector);
179
180 } // end of Tkin loop`
181 fPAIxscBank.push_back(PAItransferTable);
182 fPAIdEdxBank.push_back(PAIdEdxTable);
183 //G4cout << "dEdxMeanVector: " << G4endl;
184 //G4cout << *dEdxMeanVector << G4endl;
185 fdEdxTable.push_back(dEdxMeanVector);
186}
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:642
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition G4SandiaTable.cc:949

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 276 of file G4PAIModelData.cc.

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

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

Referenced by G4PAIModel::SampleSecondaries().

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