Geant4 11.1.1
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:57
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 307 of file G4PAIPhotData.cc.

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

Referenced by G4PAIPhotModel::CrossSectionPerVolume().

◆ DEDXPerVolume()

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

Definition at line 271 of file G4PAIPhotData.cc.

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

Referenced by G4PAIPhotModel::ComputeDEDXPerVolume().

◆ GetPlasmonRatio()

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

Definition at line 370 of file G4PAIPhotData.cc.

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

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 = fPAIxSection.GetMeanEnergyLoss();// total <dE/dx>
226
227 if(ionloss < 0.0) ionloss = 0.0;
228
229 dEdxMeanVector->PutValue(i,ionloss);
230
231 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
232 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
233 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
234
235 G4double dEdxCut = dEdxVector->Value(cut)/cut;
236 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
237
238 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
239 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
240 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
241
242 dNdxCutVector->PutValue(i, dNdxCut);
243 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
244 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
245
246 dEdxCutVector->PutValue(i, dEdxCut);
247
248 PAItransferTable->insertAt(i,transferVector);
249 PAIphotonTable->insertAt(i,photonVector);
250 PAIplasmonTable->insertAt(i,plasmonVector);
251 PAIdEdxTable->insertAt(i,dEdxVector);
252
253 } // end of Tkin loop
254
255 fPAIxscBank.push_back(PAItransferTable);
256 fPAIphotonBank.push_back(PAIphotonTable);
257 fPAIplasmonBank.push_back(PAIplasmonTable);
258
259 fPAIdEdxBank.push_back(PAIdEdxTable);
260 fdEdxTable.push_back(dEdxMeanVector);
261
262 fdNdxCutTable.push_back(dNdxCutVector);
263 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
264 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
265
266 fdEdxCutTable.push_back(dEdxCutVector);
267}
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:301
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIxSection.hh:171
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIxSection.hh:295
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition: G4PAIxSection.hh:313
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIxSection.hh:283
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition: G4PAIxSection.hh:289
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIxSection.hh:183
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:246
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition: G4ProductionCutsTable.hh:239
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:233
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
G4SandiaTable::Initialize
void Initialize(const G4Material *)
Definition: G4SandiaTable.cc:715
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1066

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 527 of file G4PAIPhotData.cc.

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

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

Referenced by G4PAIPhotModel::SampleFluctuations().

◆ SampleAlongStepTransfer()

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

Definition at line 435 of file G4PAIPhotData.cc.

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

◆ SamplePostStepPhotonTransfer()

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

Definition at line 770 of file G4PAIPhotData.cc.

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

Referenced by G4PAIPhotModel::SampleSecondaries().

◆ SamplePostStepPlasmonTransfer()

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

Definition at line 828 of file G4PAIPhotData.cc.

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

Referenced by G4PAIPhotModel::SampleSecondaries().

◆ SamplePostStepTransfer()

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

Definition at line 714 of file G4PAIPhotData.cc.

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