Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PAIPhotData.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: G4PAIPhotData.cc
31//
32// Author: V.Grichine based on G4PAIModelData MT
33//
34// Creation date: 07.10.2013
35//
36// Modifications:
37//
38
39#include "G4PAIPhotData.hh"
40#include "G4PAIPhotModel.hh"
41
42#include "G4SystemOfUnits.hh"
44
45#include "G4PhysicsLogVector.hh"
47#include "G4PhysicsTable.hh"
50#include "G4SandiaTable.hh"
51#include "Randomize.hh"
52#include "G4Poisson.hh"
53
54////////////////////////////////////////////////////////////////////////
55
56using namespace std;
57
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}
90
91////////////////////////////////////////////////////////////////////////////
92
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}
123
124///////////////////////////////////////////////////////////////////////////////
125
127 G4double cut, G4PAIPhotModel* model)
128{
129 G4ProductionCutsTable* theCoupleTable=
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}
268
269//////////////////////////////////////////////////////////////////////////////
270
272 G4double cut) const
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}
303
304/////////////////////////////////////////////////////////////////////////
305
308 G4double scaledTkin,
309 G4double tcut, G4double tmax) const
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}
366
367/////////////////////////////////////////////////////////////////////////
368
370G4PAIPhotData::GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
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}
432
433///////////////////////////////////////////////////////////////////////
434
436 G4double kinEnergy,
437 G4double scaledTkin,
438 G4double stepFactor) const
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}
524
525////////////////////////////////////////////////////////////////////////
526
528 G4double kinEnergy,
529 G4double scaledTkin,
530 G4double stepFactor) const
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}
616
617//////////////////////////////////////////////////////////////////
618
620 G4double kinEnergy,
621 G4double scaledTkin,
622 G4double stepFactor) const
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}
708
709///////////////////////////////////////////////////////////////////////
710//
711// Returns post step PAI energy transfer > cut electron energy
712// according to passed scaled kinetic energy of particle
713
715 G4double scaledTkin) const
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}
767
768////////////////////////////////////////////////////////////////////////
769
771 G4double scaledTkin) const
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}
825
826//////////////////////////////////////////////////////////////////////////
827
829 G4double scaledTkin) const
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}
883
884///////////////////////////////////////////////////////////////////////
885//
886// Returns PAI energy transfer according to passed
887// indexes of particle kinetic enegry and random x-section
888
889G4double G4PAIPhotData::GetEnergyTransfer(G4int coupleIndex,
890 std::size_t iPlace,
891 G4double position) const
892{
893 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
894 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
895
896 std::size_t iTransferMax = v->GetVectorLength() - 1;
897
898 std::size_t iTransfer;
899 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
900
901 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
902 x2 = v->Energy(iTransfer);
903 y2 = (*v)[iTransfer]/x2;
904 if(position >= y2) { break; }
905 }
906
907 x1 = v->Energy(iTransfer-1);
908 y1 = (*v)[iTransfer-1]/x1;
909 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
910 // << " x1= " << x1 << " x2= " << x2 << G4endl;
911
912 energyTransfer = x1;
913 if ( x1 != x2 ) {
914 if ( y1 == y2 ) {
915 energyTransfer += (x2 - x1)*G4UniformRand();
916 } else {
917 if(x1*1.1 < x2) {
918 const G4int nbins = 5;
919 G4double del = (x2 - x1)/G4int(nbins);
920 x2 = x1;
921 for(G4int i=1; i<=nbins; ++i) {
922 x2 += del;
923 y2 = v->Value(x2)/x2;
924 if(position >= y2) { break; }
925 x1 = x2;
926 y1 = y2;
927 }
928 }
929 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
930 }
931 }
932 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
933 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
934 // << " E(keV)= " << energyTransfer/keV << G4endl;
935 return energyTransfer;
936}
937
938/////////////////////////////////////////////////////////////////
939
940G4double G4PAIPhotData::GetEnergyPhotonTransfer(G4int coupleIndex,
941 std::size_t iPlace,
942 G4double position) const
943{
944 G4PhysicsVector* v = (*(fPAIphotonBank[coupleIndex]))(iPlace);
945 if(position*v->Energy(0) >= (*v)[0]) return v->Energy(0);
946
947 std::size_t iTransferMax = v->GetVectorLength() - 1;
948
949 std::size_t iTransfer;
950 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
951
952 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
953 {
954 x2 = v->Energy(iTransfer);
955 y2 = (*v)[iTransfer]/x2;
956 if(position >= y2) break;
957 }
958 x1 = v->Energy(iTransfer-1);
959 y1 = (*v)[iTransfer-1]/x1;
960
961 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
962 // << " x1= " << x1 << " x2= " << x2 << G4endl;
963
964 energyTransfer = x1;
965
966 if ( x1 != x2 )
967 {
968 if ( y1 == y2 )
969 {
970 energyTransfer += (x2 - x1)*G4UniformRand();
971 }
972 else
973 {
974 if( x1*1.1 < x2 )
975 {
976 const G4int nbins = 5;
977 G4double del = (x2 - x1)/G4int(nbins);
978 x2 = x1;
979
980 for(G4int i=1; i<=nbins; ++i)
981 {
982 x2 += del;
983 y2 = v->Value(x2)/x2;
984 if(position >= y2) { break; }
985 x1 = x2;
986 y1 = y2;
987 }
988 }
989 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
990 }
991 }
992 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
993 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
994 // << " E(keV)= " << energyTransfer/keV << G4endl;
995
996 return energyTransfer;
997}
998
999/////////////////////////////////////////////////////////////////////////
1000
1001G4double G4PAIPhotData::GetEnergyPlasmonTransfer(G4int coupleIndex,
1002 std::size_t iPlace,
1003 G4double position) const
1004{
1005 G4PhysicsVector* v = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
1006
1007 if( position*v->Energy(0) >= (*v)[0] ) return v->Energy(0);
1008
1009 std::size_t iTransferMax = v->GetVectorLength() - 1;
1010
1011 std::size_t iTransfer;
1012 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1013
1014 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1015 {
1016 x2 = v->Energy(iTransfer);
1017 y2 = (*v)[iTransfer]/x2;
1018 if(position >= y2) break;
1019 }
1020 x1 = v->Energy(iTransfer-1);
1021 y1 = (*v)[iTransfer-1]/x1;
1022
1023 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
1024 // << " x1= " << x1 << " x2= " << x2 << G4endl;
1025
1026 energyTransfer = x1;
1027
1028 if ( x1 != x2 )
1029 {
1030 if ( y1 == y2 )
1031 {
1032 energyTransfer += (x2 - x1)*G4UniformRand();
1033 }
1034 else
1035 {
1036 if(x1*1.1 < x2)
1037 {
1038 const G4int nbins = 5;
1039 G4double del = (x2 - x1)/G4int(nbins);
1040 x2 = x1;
1041
1042 for( G4int i = 1; i <= nbins; ++i )
1043 {
1044 x2 += del;
1045 y2 = v->Value(x2)/x2;
1046
1047 if(position >= y2) break;
1048
1049 x1 = x2;
1050 y1 = y2;
1051 }
1052 }
1053 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1054 }
1055 }
1056 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
1057 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
1058 // << " E(keV)= " << energyTransfer/keV << G4endl;
1059
1060 return energyTransfer;
1061}
1062
1063//
1064//
1065//////////////////////////////////////////////////////////////////////
1066
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
@ idxG4ElectronCut
@ idxG4GammaCut
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 CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double ComputeMaxEnergy(G4double scaledEnergy)
G4double GetIntegralCerenkov(G4int i) const
G4int GetSplineSize() const
G4double GetIntegralPAIdEdx(G4int i) const
G4double GetIntegralPlasmon(G4int i) const
G4double GetSplineEnergy(G4int i) const
G4double GetIntegralPAIxSection(G4int i) const
G4double GetMeanEnergyLoss() 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
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
static G4ProductionCutsTable * GetProductionCutsTable()
void Initialize(const G4Material *)
G4double GetSandiaMatTablePAI(G4int, G4int) const
#define W
Definition: crc32.c:84