Geant4 10.7.0
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
 

Detailed Description

Definition at line 67 of file G4PAIModelData.hh.

Constructor & Destructor Documentation

◆ G4PAIModelData()

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:142

◆ ~G4PAIModelData()

G4PAIModelData::~G4PAIModelData ( )

Definition at line 89 of file G4PAIModelData.cc.

90{
91 size_t n = fPAIxscBank.size();
92 if(0 < n) {
93 for(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}

Member Function Documentation

◆ CrossSectionPerVolume()

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

Definition at line 244 of file G4PAIModelData.cc.

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

Referenced by G4PAIModel::CrossSectionPerVolume().

◆ DEDXPerVolume()

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

Definition at line 202 of file G4PAIModelData.cc.

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

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 G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
117 G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
118 G4PhysicsLogVector* 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 G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
154 G4PhysicsFreeVector* 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 //transferVector->SetSpline(true);
184 //transferVector->FillSecondDerivatives();
185 //dEdxVector->SetSpline(true);
186 //dEdxVector->FillSecondDerivatives();
187
188 } // end of Tkin loop`
189 fPAIxscBank.push_back(PAItransferTable);
190 fPAIdEdxBank.push_back(PAIdEdxTable);
191 //G4cout << "dEdxMeanVector: " << G4endl;
192 //G4cout << *dEdxMeanVector << G4endl;
193 /*
194 dEdxMeanVector->SetSpline(true);
195 dEdxMeanVector->FillSecondDerivatives();
196 */
197 fdEdxTable.push_back(dEdxMeanVector);
198}
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4PAIModel::ComputeMaxEnergy
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIModel.hh:164
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:127
G4PAIySection::GetIntegralPAIySection
G4double GetIntegralPAIySection(G4int i) const
Definition: G4PAIySection.hh:209
G4PAIySection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIySection.hh:117
G4PAIySection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIySection.cc:118
G4PhysicsFreeVector::PutValue
void PutValue(std::size_t index, G4double energy, G4double dValue)
Definition: G4PhysicsFreeVector.hh:77
G4PhysicsTable::insertAt
void insertAt(std::size_t, G4PhysicsVector *)
G4PhysicsVector::PutValue
void PutValue(std::size_t index, G4double theValue)
G4SandiaTable::Initialize
void Initialize(const G4Material *)
Definition: G4SandiaTable.cc:705
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1031

Referenced by G4PAIModel::Initialise().

◆ SampleAlongStepTransfer()

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

Definition at line 288 of file G4PAIModelData.cc.

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

Referenced by G4PAIModel::SampleFluctuations().

◆ SamplePostStepTransfer()

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

Definition at line 382 of file G4PAIModelData.cc.

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

Referenced by G4PAIModel::SampleSecondaries().

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