Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PAIModelData.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class
30// File name: G4PAIModelData.cc
31//
32// Author: V. Ivanchenko based on V.Grichine code of G4PAIModel
33//
34// Creation date: 16.08.2013
35//
36// Modifications:
37//
38
39#include "G4PAIModelData.hh"
40#include "G4PAIModel.hh"
41
42#include "G4SystemOfUnits.hh"
44
45#include "G4PhysicsLogVector.hh"
47#include "G4PhysicsTable.hh"
49#include "G4SandiaTable.hh"
50#include "Randomize.hh"
51#include "G4Poisson.hh"
52
53////////////////////////////////////////////////////////////////////////
54
55using namespace std;
56
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}
86
87////////////////////////////////////////////////////////////////////////////
88
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}
107
108///////////////////////////////////////////////////////////////////////////////
109
111 G4PAIModel* model)
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}
190
191//////////////////////////////////////////////////////////////////////////////
192
194 G4double cut) const
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}
231
232/////////////////////////////////////////////////////////////////////////
233
236 G4double scaledTkin,
237 G4double tcut, G4double tmax) const
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}
276
277///////////////////////////////////////////////////////////////////////
278
280 G4double kinEnergy,
281 G4double scaledTkin,
282 G4double tmax,
283 G4double stepFactor) const
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}
367
368///////////////////////////////////////////////////////////////////////
369//
370// Returns post step PAI energy transfer > cut electron energy
371// according to passed scaled kinetic energy of particle
372
374 G4double scaledTkin,
375 G4double tmin,
376 G4double tmax) const
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}
447
448///////////////////////////////////////////////////////////////////////
449//
450// Returns PAI energy transfer according to passed
451// indexes of particle kinetic enegry and random x-section
452
453G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
454 std::size_t iPlace,
455 G4double position) const
456{
457 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
458 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
459
460 std::size_t iTransferMax = v->GetVectorLength() - 1;
461
462 std::size_t iTransfer;
463 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
464
465 //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl;
466 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
467 x2 = v->Energy(iTransfer);
468 y2 = (*v)[iTransfer]/x2;
469 if(position >= y2) { break; }
470 if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); }
471 }
472
473 x1 = v->Energy(iTransfer-1);
474 y1 = (*v)[iTransfer-1]/x1;
475 /*
476 G4cout << "i= " << iTransfer << " imax= " << iTransferMax
477 << " x1= " << x1 << " x2= " << x2
478 << " y1= " << y1 << " y2= " << y2 << G4endl;
479 */
480 energyTransfer = x1;
481 if ( x1 != x2 ) {
482 if ( y1 == y2 ) {
483 energyTransfer += (x2 - x1)*G4UniformRand();
484 } else {
485 if(x1*1.1 < x2) {
486 const G4int nbins = 5;
487 G4double del = (x2 - x1)/G4int(nbins);
488 x2 = x1;
489 for(G4int i=1; i<=nbins; ++i) {
490 x2 += del;
491 y2 = v->Value(x2)/x2;
492 if(position >= y2) {
493 break;
494 }
495 x1 = x2;
496 y1 = y2;
497 }
498 }
499 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
500 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl;
501 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
502 }
503 }
504 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
505 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
506 // << " E(keV)= " << energyTransfer/keV << G4endl;
507 return energyTransfer;
508}
509
510//////////////////////////////////////////////////////////////////////
511
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
const G4Material * GetMaterial() const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
void Initialise(const G4MaterialCutsCouple *, G4PAIModel *)
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin, G4double tmin, G4double tmax) const
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double tmax, G4double stepFactor) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIModel.hh:166
G4double GetSplineEnergy(G4int i) const
void SetVerbose(G4int v)
G4double GetIntegralPAIdEdx(G4int i) const
G4double GetMeanEnergyLoss() const
G4double GetIntegralPAIySection(G4int i) const
G4int GetSplineSize() const
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
G4double GetMaxEnergy() const
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
std::size_t GetVectorLength() const
std::size_t FindBin(const G4double energy, std::size_t idx) const
void Initialize(const G4Material *)
G4double GetSandiaMatTablePAI(G4int, G4int) const
#define W
Definition: crc32.c:84