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
G4PAIPhotData Class Reference

#include <G4PAIPhotData.hh>

Public Member Functions

 G4PAIPhotData (G4double tmin, G4double tmax, G4int verbose)
 
 ~G4PAIPhotData ()
 
void Initialise (const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
 
G4double DEDXPerVolume (G4int coupleIndex, G4double scaledTkin, G4double cut) const
 
G4double CrossSectionPerVolume (G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
 
G4double GetPlasmonRatio (G4int coupleIndex, G4double scaledTkin) const
 
G4double SampleAlongStepTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SampleAlongStepPhotonTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SampleAlongStepPlasmonTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SamplePostStepTransfer (G4int coupleIndex, G4double scaledTkin) const
 
G4double SamplePostStepPhotonTransfer (G4int coupleIndex, G4double scaledTkin) const
 
G4double SamplePostStepPlasmonTransfer (G4int coupleIndex, G4double scaledTkin) const
 
G4PAIPhotData & operator= (const G4PAIPhotData &right)=delete
 
 G4PAIPhotData (const G4PAIPhotData &)=delete
 

Detailed Description

Definition at line 65 of file G4PAIPhotData.hh.

Constructor & Destructor Documentation

◆ G4PAIPhotData() [1/2]

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

Definition at line 58 of file G4PAIPhotData.cc.

59{
60 const G4int nPerDecade = 10;
61 const G4double lowestTkin = 50*keV;
62 const G4double highestTkin = 10*TeV;
63
64 // fPAIxSection.SetVerbose(ver);
65
66 fLowestKineticEnergy = std::max(tmin, lowestTkin);
67 fHighestKineticEnergy = tmax;
68
69 if(tmax < 10*fLowestKineticEnergy)
70 {
71 fHighestKineticEnergy = 10*fLowestKineticEnergy;
72 }
73 else if(tmax > highestTkin)
74 {
75 fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
76 }
77 fTotBin = (G4int)(nPerDecade*
78 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
79
80 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
81 fHighestKineticEnergy,
82 fTotBin);
83 if(0 < ver) {
84 G4cout << "### G4PAIPhotData: Nbins= " << fTotBin
85 << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
86 << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
87 << " tmin(keV)= " << tmin/keV << G4endl;
88 }
89}
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

◆ ~G4PAIPhotData()

G4PAIPhotData::~G4PAIPhotData ( )

Definition at line 93 of file G4PAIPhotData.cc.

94{
95 //G4cout << "G4PAIPhotData::~G4PAIPhotData() " << this << G4endl;
96 std::size_t n = fPAIxscBank.size();
97 if(0 < n)
98 {
99 for(std::size_t i=0; i<n; ++i)
100 {
101 if(fPAIxscBank[i])
102 {
103 fPAIxscBank[i]->clearAndDestroy();
104 delete fPAIxscBank[i];
105 fPAIxscBank[i] = nullptr;
106 }
107 if(fPAIdEdxBank[i])
108 {
109 fPAIdEdxBank[i]->clearAndDestroy();
110 delete fPAIdEdxBank[i];
111 fPAIdEdxBank[i] = nullptr;
112 }
113 delete fdEdxTable[i];
114 delete fdNdxCutTable[i];
115 fdEdxTable[i] = nullptr;
116 fdNdxCutTable[i] = nullptr;
117 }
118 }
119 delete fParticleEnergyVector;
120 fParticleEnergyVector = nullptr;
121 //G4cout << "G4PAIPhotData::~G4PAIPhotData() done for " << this << G4endl;
122}

◆ G4PAIPhotData() [2/2]

G4PAIPhotData::G4PAIPhotData ( const G4PAIPhotData & )
delete

Member Function Documentation

◆ CrossSectionPerVolume()

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

Definition at line 304 of file G4PAIPhotData.cc.

307{
308 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
309
310 cross=tcut+tmax;
311
312 // iPlace is in interval from 0 to (N-1)
313
314 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
315 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
316
317 G4bool one = true;
318
319 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
320 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
321
322
323 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
324 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
325
326 xscPh = xscPh2;
327 xscEl = xscEl2;
328
329 cross = xscPh + xscEl;
330
331 if( !one )
332 {
333 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
334
335 G4double E1 = fParticleEnergyVector->Energy(iPlace);
336 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
337
338 G4double W = 1.0/(E2 - E1);
339 G4double W1 = (E2 - scaledTkin)*W;
340 G4double W2 = (scaledTkin - E1)*W;
341
342 xscEl *= W1;
343 xscEl += W2*xscEl2;
344
345 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
346
347 E1 = fParticleEnergyVector->Energy(iPlace);
348 E2 = fParticleEnergyVector->Energy(iPlace+1);
349
350 W = 1.0/(E2 - E1);
351 W1 = (E2 - scaledTkin)*W;
352 W2 = (scaledTkin - E1)*W;
353
354 xscPh *= W1;
355 xscPh += W2*xscPh2;
356
357 cross = xscEl + xscPh;
358 }
359 if( cross < 0.0) cross = 0.0;
360
361 return cross;
362}
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 G4PAIPhotModel::CrossSectionPerVolume().

◆ DEDXPerVolume()

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

Definition at line 268 of file G4PAIPhotData.cc.

270{
271 // VI: iPlace is the low edge index of the bin
272 // iPlace is in interval from 0 to (N-1)
273 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
274 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
275
276 G4bool one = true;
277 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
278 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
279 one = false;
280 }
281
282 // VI: apply interpolation of the vector
283 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
284 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
285 if(!one) {
286 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
287 G4double E1 = fParticleEnergyVector->Energy(iPlace);
288 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
289 G4double W = 1.0/(E2 - E1);
290 G4double W1 = (E2 - scaledTkin)*W;
291 G4double W2 = (scaledTkin - E1)*W;
292 del *= W1;
293 del += W2*del2;
294 }
295 dEdx -= del;
296
297 if( dEdx < 0.) { dEdx = 0.; }
298 return dEdx;
299}

Referenced by G4PAIPhotModel::ComputeDEDXPerVolume().

◆ GetPlasmonRatio()

G4double G4PAIPhotData::GetPlasmonRatio ( G4int coupleIndex,
G4double scaledTkin ) const

Definition at line 367 of file G4PAIPhotData.cc.

368{
369 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
370 // iPlace is in interval from 0 to (N-1)
371
372 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
373 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
374
375 G4bool one = true;
376
377 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
378 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
379
380
381 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
382 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
383
384 xscPh = xscPh2;
385 xscEl = xscEl2;
386
387 cross = xscPh + xscEl;
388
389 if( !one )
390 {
391 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
392
393 G4double E1 = fParticleEnergyVector->Energy(iPlace);
394 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
395
396 G4double W = 1.0/(E2 - E1);
397 G4double W1 = (E2 - scaledTkin)*W;
398 G4double W2 = (scaledTkin - E1)*W;
399
400 xscEl *= W1;
401 xscEl += W2*xscEl2;
402
403 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
404
405 E1 = fParticleEnergyVector->Energy(iPlace);
406 E2 = fParticleEnergyVector->Energy(iPlace+1);
407
408 W = 1.0/(E2 - E1);
409 W1 = (E2 - scaledTkin)*W;
410 W2 = (scaledTkin - E1)*W;
411
412 xscPh *= W1;
413 xscPh += W2*xscPh2;
414
415 cross = xscEl + xscPh;
416 }
417 if( cross <= 0.0)
418 {
419 plRatio = 2.0;
420 }
421 else
422 {
423 plRatio = xscEl/cross;
424
425 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
426 }
427 return plRatio;
428}

Referenced by G4PAIPhotModel::SampleSecondaries().

◆ Initialise()

void G4PAIPhotData::Initialise ( const G4MaterialCutsCouple * couple,
G4double cut,
G4PAIPhotModel * model )

Definition at line 126 of file G4PAIPhotData.cc.

128{
129 G4ProductionCutsTable* theCoupleTable=
130 G4ProductionCutsTable::GetProductionCutsTable();
131 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
132 G4int jMatCC;
133
134 for (jMatCC = 0; jMatCC < numOfCouples; ++jMatCC )
135 {
136 if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
137 }
138 if( jMatCC == numOfCouples && jMatCC > 0 ) --jMatCC;
139
140 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
141 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
142 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
143 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
144 /*
145 G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
146 <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
147 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
148 */
149 // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??
150
151 auto dEdxCutVector =
152 new G4PhysicsLogVector(fLowestKineticEnergy,
153 fHighestKineticEnergy,
154 fTotBin);
155
156 auto dNdxCutVector =
157 new G4PhysicsLogVector(fLowestKineticEnergy,
158 fHighestKineticEnergy,
159 fTotBin);
160 auto dNdxCutPhotonVector =
161 new G4PhysicsLogVector(fLowestKineticEnergy,
162 fHighestKineticEnergy,
163 fTotBin);
164 auto dNdxCutPlasmonVector =
165 new G4PhysicsLogVector(fLowestKineticEnergy,
166 fHighestKineticEnergy,
167 fTotBin);
168
169 const G4Material* mat = couple->GetMaterial();
170 fSandia.Initialize(const_cast<G4Material*>(mat));
171
172 auto PAItransferTable = new G4PhysicsTable(fTotBin+1);
173 auto PAIphotonTable = new G4PhysicsTable(fTotBin+1);
174 auto PAIplasmonTable = new G4PhysicsTable(fTotBin+1);
175
176 auto PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
177 auto dEdxMeanVector =
178 new G4PhysicsLogVector(fLowestKineticEnergy,
179 fHighestKineticEnergy,
180 fTotBin);
181
182 // low energy Sandia interval
183 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
184
185 // energy safety
186 const G4double deltaLow = 100.*eV;
187
188 for (G4int i = 0; i <= fTotBin; ++i)
189 {
190 G4double kinEnergy = fParticleEnergyVector->Energy(i);
191 G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
192 G4double tau = kinEnergy/proton_mass_c2;
193 G4double bg2 = tau*( tau + 2. );
194
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
196
197 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
198
199 //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
200 // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
201
202 G4int n = fPAIxSection.GetSplineSize();
203
204 auto transferVector = new G4PhysicsFreeVector(n);
205 auto photonVector = new G4PhysicsFreeVector(n);
206 auto plasmonVector = new G4PhysicsFreeVector(n);
207
208 auto dEdxVector = new G4PhysicsFreeVector(n);
209
210 for( G4int k = 0; k < n; k++ )
211 {
212 G4double t = fPAIxSection.GetSplineEnergy(k+1);
213
214 transferVector->PutValue(k , t,
215 t*fPAIxSection.GetIntegralPAIxSection(k+1));
216 photonVector->PutValue(k , t,
217 t*fPAIxSection.GetIntegralCerenkov(k+1));
218 plasmonVector->PutValue(k , t,
219 t*fPAIxSection.GetIntegralPlasmon(k+1));
220
221 dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
222 }
223 // G4cout << *transferVector << G4endl;
224
225 G4double ionloss = std::max(fPAIxSection.GetMeanEnergyLoss(), 0.0);// total <dE/dx>
226 dEdxMeanVector->PutValue(i,ionloss);
227
228 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
229 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
230 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
231
232 G4double dEdxCut = dEdxVector->Value(cut)/cut;
233 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
234
235 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
236 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
237 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
238
239 dNdxCutVector->PutValue(i, dNdxCut);
240 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
241 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
242
243 dEdxCutVector->PutValue(i, dEdxCut);
244
245 PAItransferTable->insertAt(i,transferVector);
246 PAIphotonTable->insertAt(i,photonVector);
247 PAIplasmonTable->insertAt(i,plasmonVector);
248 PAIdEdxTable->insertAt(i,dEdxVector);
249
250 } // end of Tkin loop
251
252 fPAIxscBank.push_back(PAItransferTable);
253 fPAIphotonBank.push_back(PAIphotonTable);
254 fPAIplasmonBank.push_back(PAIplasmonTable);
255
256 fPAIdEdxBank.push_back(PAIdEdxTable);
257 fdEdxTable.push_back(dEdxMeanVector);
258
259 fdNdxCutTable.push_back(dNdxCutVector);
260 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
261 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
262
263 fdEdxCutTable.push_back(dEdxCutVector);
264}
idxG4ElectronCut
@ idxG4ElectronCut
Definition G4ProductionCuts.hh:48
idxG4GammaCut
@ idxG4GammaCut
Definition G4ProductionCuts.hh:47
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4PAIPhotModel::ComputeMaxEnergy
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition G4PAIPhotModel.hh:160
G4PAIxSection::GetIntegralCerenkov
G4double GetIntegralCerenkov(G4int i) const
Definition G4PAIxSection.hh:298
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition G4PAIxSection.hh:170
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition G4PAIxSection.hh:292
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition G4PAIxSection.hh:310
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition G4PAIxSection.hh:280
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition G4PAIxSection.hh:286
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition G4PAIxSection.hh:182
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition G4PAIxSection.cc:606
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
Definition G4ProductionCutsTable.hh:243
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4SandiaTable::Initialize
void Initialize(const G4Material *)
Definition G4SandiaTable.cc:642
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition G4SandiaTable.cc:949

Referenced by G4PAIPhotModel::Initialise().

◆ operator=()

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

◆ SampleAlongStepPhotonTransfer()

G4double G4PAIPhotData::SampleAlongStepPhotonTransfer ( G4int coupleIndex,
G4double kinEnergy,
G4double scaledTkin,
G4double stepFactor ) const

Definition at line 524 of file G4PAIPhotData.cc.

528{
529 G4double loss = 0.0;
530 G4double omega;
531 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
532
533 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
534 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
535
536 G4bool one = true;
537
538 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
539 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
540
541 G4PhysicsLogVector* vcut = fdNdxCutPhotonTable[coupleIndex];
542 G4PhysicsVector* v1 = (*(fPAIphotonBank[coupleIndex]))(iPlace);
543 G4PhysicsVector* v2 = nullptr;
544
545 dNdxCut1 = (*vcut)[iPlace];
546 G4double e1 = v1->Energy(0);
547 G4double e2 = e1;
548
549 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
550
551 meanNumber = meanN1;
552
553 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
554 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
555
556 dNdxCut2 = dNdxCut1;
557 W1 = 1.0;
558 W2 = 0.0;
559 if(!one)
560 {
561 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
562 dNdxCut2 = (*vcut)[iPlace+1];
563 e2 = v2->Energy(0);
564
565 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
566
567 E1 = fParticleEnergyVector->Energy(iPlace);
568 E2 = fParticleEnergyVector->Energy(iPlace+1);
569 W = 1.0/(E2 - E1);
570 W1 = (E2 - scaledTkin)*W;
571 W2 = (scaledTkin - E1)*W;
572 meanNumber = W1*meanN1 + W2*meanN2;
573
574 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
575 // << " dNdxCut2= " << dNdxCut2 << G4endl;
576 }
577 if( meanNumber <= 0.0) return 0.0;
578
579 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
580
581 //G4cout << "N= " << numOfCollisions << G4endl;
582
583 if( 0 == numOfCollisions) return 0.0;
584
585 for(G4int i=0; i< numOfCollisions; ++i)
586 {
587 G4double rand = G4UniformRand();
588 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
589 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
590
591 //G4cout << "omega(keV)= " << omega/keV << G4endl;
592
593 if(!one)
594 {
595 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
596 G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
597 omega = omega*W1 + omega2*W2;
598 }
599 //G4cout << "omega(keV)= " << omega/keV << G4endl;
600
601 loss += omega;
602 if( loss > kinEnergy) { break; }
603 }
604
605 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
606 //<<step/mm<<" mm"<<G4endl;
607
608 if ( loss > kinEnergy) loss = kinEnergy;
609 else if( loss < 0.) loss = 0.;
610
611 return loss;
612}
G4Poisson
G4long G4Poisson(G4double mean)
Definition G4Poisson.hh:50
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52

Referenced by G4PAIPhotModel::SampleFluctuations().

◆ SampleAlongStepPlasmonTransfer()

G4double G4PAIPhotData::SampleAlongStepPlasmonTransfer ( G4int coupleIndex,
G4double kinEnergy,
G4double scaledTkin,
G4double stepFactor ) const

Definition at line 616 of file G4PAIPhotData.cc.

620{
621 G4double loss = 0.0;
622 G4double omega;
623 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
624
625 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
626 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
627
628 G4bool one = true;
629
630 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
631 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
632
633 G4PhysicsLogVector* vcut = fdNdxCutPlasmonTable[coupleIndex];
634 G4PhysicsVector* v1 = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
635 G4PhysicsVector* v2 = nullptr;
636
637 dNdxCut1 = (*vcut)[iPlace];
638 G4double e1 = v1->Energy(0);
639 G4double e2 = e1;
640
641 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
642
643 meanNumber = meanN1;
644
645 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
646 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
647
648 dNdxCut2 = dNdxCut1;
649 W1 = 1.0;
650 W2 = 0.0;
651 if(!one)
652 {
653 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
654 dNdxCut2 = (*vcut)[iPlace+1];
655 e2 = v2->Energy(0);
656
657 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
658
659 E1 = fParticleEnergyVector->Energy(iPlace);
660 E2 = fParticleEnergyVector->Energy(iPlace+1);
661 W = 1.0/(E2 - E1);
662 W1 = (E2 - scaledTkin)*W;
663 W2 = (scaledTkin - E1)*W;
664 meanNumber = W1*meanN1 + W2*meanN2;
665
666 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
667 // << " dNdxCut2= " << dNdxCut2 << G4endl;
668 }
669 if( meanNumber <= 0.0) return 0.0;
670
671 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
672
673 //G4cout << "N= " << numOfCollisions << G4endl;
674
675 if( 0 == numOfCollisions) return 0.0;
676
677 for(G4int i=0; i< numOfCollisions; ++i)
678 {
679 G4double rand = G4UniformRand();
680 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
681 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
682
683 //G4cout << "omega(keV)= " << omega/keV << G4endl;
684
685 if(!one)
686 {
687 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
688 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
689 omega = omega*W1 + omega2*W2;
690 }
691 //G4cout << "omega(keV)= " << omega/keV << G4endl;
692
693 loss += omega;
694 if( loss > kinEnergy) { break; }
695 }
696
697 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
698 //<<step/mm<<" mm"<<G4endl;
699
700 if ( loss > kinEnergy) loss = kinEnergy;
701 else if( loss < 0.) loss = 0.;
702
703 return loss;
704}

Referenced by G4PAIPhotModel::SampleFluctuations().

◆ SampleAlongStepTransfer()

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

Definition at line 432 of file G4PAIPhotData.cc.

436{
437 G4double loss = 0.0;
438 G4double omega;
439 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
440
441 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
442 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
443
444 G4bool one = true;
445
446 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
447 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
448
449 G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
450 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
451 G4PhysicsVector* v2 = nullptr;
452
453 dNdxCut1 = (*vcut)[iPlace];
454 G4double e1 = v1->Energy(0);
455 G4double e2 = e1;
456
457 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
458
459 meanNumber = meanN1;
460
461 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
462 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
463
464 dNdxCut2 = dNdxCut1;
465 W1 = 1.0;
466 W2 = 0.0;
467 if(!one)
468 {
469 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
470 dNdxCut2 = (*vcut)[iPlace+1];
471 e2 = v2->Energy(0);
472
473 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
474
475 E1 = fParticleEnergyVector->Energy(iPlace);
476 E2 = fParticleEnergyVector->Energy(iPlace+1);
477 W = 1.0/(E2 - E1);
478 W1 = (E2 - scaledTkin)*W;
479 W2 = (scaledTkin - E1)*W;
480 meanNumber = W1*meanN1 + W2*meanN2;
481
482 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
483 // << " dNdxCut2= " << dNdxCut2 << G4endl;
484 }
485 if( meanNumber <= 0.0) return 0.0;
486
487 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
488
489 //G4cout << "N= " << numOfCollisions << G4endl;
490
491 if( 0 == numOfCollisions) return 0.0;
492
493 for(G4int i=0; i< numOfCollisions; ++i)
494 {
495 G4double rand = G4UniformRand();
496 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
497 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
498
499 //G4cout << "omega(keV)= " << omega/keV << G4endl;
500
501 if(!one)
502 {
503 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
504 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
505 omega = omega*W1 + omega2*W2;
506 }
507 //G4cout << "omega(keV)= " << omega/keV << G4endl;
508
509 loss += omega;
510 if( loss > kinEnergy) { break; }
511 }
512
513 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
514 //<<step/mm<<" mm"<<G4endl;
515
516 if ( loss > kinEnergy) loss = kinEnergy;
517 else if( loss < 0.) loss = 0.;
518
519 return loss;
520}

◆ SamplePostStepPhotonTransfer()

G4double G4PAIPhotData::SamplePostStepPhotonTransfer ( G4int coupleIndex,
G4double scaledTkin ) const

Definition at line 767 of file G4PAIPhotData.cc.

769{
770 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
771 G4double transfer = 0.0;
772 G4double rand = G4UniformRand();
773
774 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
775
776 // std::size_t iTransfer, iTr1, iTr2;
777 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
778
779 G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];
780
781 // Fermi plato, try from left
782
783 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
784 {
785 position = (*cutv)[nPlace]*rand;
786 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
787 }
788 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
789 {
790 position = (*cutv)[0]*rand;
791 transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
792 }
793 else
794 {
795 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
796
797 dNdxCut1 = (*cutv)[iPlace];
798 dNdxCut2 = (*cutv)[iPlace+1];
799
800 E1 = fParticleEnergyVector->Energy(iPlace);
801 E2 = fParticleEnergyVector->Energy(iPlace+1);
802 W = 1.0/(E2 - E1);
803 W1 = (E2 - scaledTkin)*W;
804 W2 = (scaledTkin - E1)*W;
805
806 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
807 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
808
809 position = dNdxCut1*rand;
810
811 G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
812
813 position = dNdxCut2*rand;
814 G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
815
816 transfer = tr1*W1 + tr2*W2;
817 }
818 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
819 if(transfer < 0.0 ) { transfer = 0.0; }
820 return transfer;
821}
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const

Referenced by G4PAIPhotModel::SampleSecondaries().

◆ SamplePostStepPlasmonTransfer()

G4double G4PAIPhotData::SamplePostStepPlasmonTransfer ( G4int coupleIndex,
G4double scaledTkin ) const

Definition at line 825 of file G4PAIPhotData.cc.

827{
828 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
829 G4double transfer = 0.0;
830 G4double rand = G4UniformRand();
831
832 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
833
834 // std::size_t iTransfer, iTr1, iTr2;
835 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
836
837 G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];
838
839 // Fermi plato, try from left
840 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
841 {
842 position = (*cutv)[nPlace]*rand;
843 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
844 }
845 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
846 {
847 position = (*cutv)[0]*rand;
848 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
849 }
850 else
851 {
852 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
853
854 dNdxCut1 = (*cutv)[iPlace];
855 dNdxCut2 = (*cutv)[iPlace+1];
856
857 E1 = fParticleEnergyVector->Energy(iPlace);
858 E2 = fParticleEnergyVector->Energy(iPlace+1);
859 W = 1.0/(E2 - E1);
860 W1 = (E2 - scaledTkin)*W;
861 W2 = (scaledTkin - E1)*W;
862
863 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
864 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
865
866 position = dNdxCut1*rand;
867 G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
868
869 position = dNdxCut2*rand;
870 G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
871
872 transfer = tr1*W1 + tr2*W2;
873 }
874 //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl;
875
876 if(transfer < 0.0 ) transfer = 0.0;
877
878 return transfer;
879}

Referenced by G4PAIPhotModel::SampleSecondaries().

◆ SamplePostStepTransfer()

G4double G4PAIPhotData::SamplePostStepTransfer ( G4int coupleIndex,
G4double scaledTkin ) const

Definition at line 711 of file G4PAIPhotData.cc.

713{
714 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
715 G4double transfer = 0.0;
716 G4double rand = G4UniformRand();
717
718 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
719
720 // std::size_t iTransfer, iTr1, iTr2;
721 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
722
723 G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
724
725 // Fermi plato, try from left
726 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
727 {
728 position = (*cutv)[nPlace]*rand;
729 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
730 }
731 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
732 {
733 position = (*cutv)[0]*rand;
734 transfer = GetEnergyTransfer(coupleIndex, 0, position);
735 }
736 else
737 {
738 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
739
740 dNdxCut1 = (*cutv)[iPlace];
741 dNdxCut2 = (*cutv)[iPlace+1];
742
743 E1 = fParticleEnergyVector->Energy(iPlace);
744 E2 = fParticleEnergyVector->Energy(iPlace+1);
745 W = 1.0/(E2 - E1);
746 W1 = (E2 - scaledTkin)*W;
747 W2 = (scaledTkin - E1)*W;
748
749 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
750 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
751
752 position = dNdxCut1*rand;
753 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
754
755 position = dNdxCut2*rand;
756 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
757
758 transfer = tr1*W1 + tr2*W2;
759 }
760 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
761 if(transfer < 0.0 ) { transfer = 0.0; }
762 return transfer;
763}
The documentation for this class was generated from the following files: