Geant4 11.2.2
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 = 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}
187
188//////////////////////////////////////////////////////////////////////////////
189
191 G4double cut) const
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}
228
229/////////////////////////////////////////////////////////////////////////
230
233 G4double scaledTkin,
234 G4double tcut, G4double tmax) const
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}
273
274///////////////////////////////////////////////////////////////////////
275
277 G4double kinEnergy,
278 G4double scaledTkin,
279 G4double tmax,
280 G4double stepFactor) const
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}
364
365///////////////////////////////////////////////////////////////////////
366//
367// Returns post step PAI energy transfer > cut electron energy
368// according to passed scaled kinetic energy of particle
369
371 G4double scaledTkin,
372 G4double tmin,
373 G4double tmax) const
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}
444
445///////////////////////////////////////////////////////////////////////
446//
447// Returns PAI energy transfer according to passed
448// indexes of particle kinetic enegry and random x-section
449
450G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
451 std::size_t iPlace,
452 G4double position) const
453{
454 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
455 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
456
457 std::size_t iTransferMax = v->GetVectorLength() - 1;
458
459 std::size_t iTransfer;
460 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
461
462 //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl;
463 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
464 x2 = v->Energy(iTransfer);
465 y2 = (*v)[iTransfer]/x2;
466 if(position >= y2) { break; }
467 if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); }
468 }
469
470 x1 = v->Energy(iTransfer-1);
471 y1 = (*v)[iTransfer-1]/x1;
472 /*
473 G4cout << "i= " << iTransfer << " imax= " << iTransferMax
474 << " x1= " << x1 << " x2= " << x2
475 << " y1= " << y1 << " y2= " << y2 << G4endl;
476 */
477 energyTransfer = x1;
478 if ( x1 != x2 ) {
479 if ( y1 == y2 ) {
480 energyTransfer += (x2 - x1)*G4UniformRand();
481 } else {
482 if(x1*1.1 < x2) {
483 const G4int nbins = 5;
484 G4double del = (x2 - x1)/G4int(nbins);
485 x2 = x1;
486 for(G4int i=1; i<=nbins; ++i) {
487 x2 += del;
488 y2 = v->Value(x2)/x2;
489 if(position >= y2) {
490 break;
491 }
492 x1 = x2;
493 y1 = y2;
494 }
495 }
496 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
497 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl;
498 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
499 }
500 }
501 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
502 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
503 // << " E(keV)= " << energyTransfer/keV << G4endl;
504 return energyTransfer;
505}
506
507//////////////////////////////////////////////////////////////////////
508
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:67
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)
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:85