Geant4 11.2.2
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 = 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}
265
266//////////////////////////////////////////////////////////////////////////////
267
269 G4double cut) const
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}
300
301/////////////////////////////////////////////////////////////////////////
302
305 G4double scaledTkin,
306 G4double tcut, G4double tmax) const
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}
363
364/////////////////////////////////////////////////////////////////////////
365
367G4PAIPhotData::GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
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}
429
430///////////////////////////////////////////////////////////////////////
431
433 G4double kinEnergy,
434 G4double scaledTkin,
435 G4double stepFactor) const
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}
521
522////////////////////////////////////////////////////////////////////////
523
525 G4double kinEnergy,
526 G4double scaledTkin,
527 G4double stepFactor) const
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}
613
614//////////////////////////////////////////////////////////////////
615
617 G4double kinEnergy,
618 G4double scaledTkin,
619 G4double stepFactor) const
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}
705
706///////////////////////////////////////////////////////////////////////
707//
708// Returns post step PAI energy transfer > cut electron energy
709// according to passed scaled kinetic energy of particle
710
712 G4double scaledTkin) const
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}
764
765////////////////////////////////////////////////////////////////////////
766
768 G4double scaledTkin) const
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}
822
823//////////////////////////////////////////////////////////////////////////
824
826 G4double scaledTkin) const
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}
880
881///////////////////////////////////////////////////////////////////////
882//
883// Returns PAI energy transfer according to passed
884// indexes of particle kinetic enegry and random x-section
885
886G4double G4PAIPhotData::GetEnergyTransfer(G4int coupleIndex,
887 std::size_t iPlace,
888 G4double position) const
889{
890 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
891 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
892
893 std::size_t iTransferMax = v->GetVectorLength() - 1;
894
895 std::size_t iTransfer;
896 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
897
898 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
899 x2 = v->Energy(iTransfer);
900 y2 = (*v)[iTransfer]/x2;
901 if(position >= y2) { break; }
902 }
903
904 x1 = v->Energy(iTransfer-1);
905 y1 = (*v)[iTransfer-1]/x1;
906 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
907 // << " x1= " << x1 << " x2= " << x2 << G4endl;
908
909 energyTransfer = x1;
910 if ( x1 != x2 ) {
911 if ( y1 == y2 ) {
912 energyTransfer += (x2 - x1)*G4UniformRand();
913 } else {
914 if(x1*1.1 < x2) {
915 const G4int nbins = 5;
916 G4double del = (x2 - x1)/G4int(nbins);
917 x2 = x1;
918 for(G4int i=1; i<=nbins; ++i) {
919 x2 += del;
920 y2 = v->Value(x2)/x2;
921 if(position >= y2) { break; }
922 x1 = x2;
923 y1 = y2;
924 }
925 }
926 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
927 }
928 }
929 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
930 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
931 // << " E(keV)= " << energyTransfer/keV << G4endl;
932 return energyTransfer;
933}
934
935/////////////////////////////////////////////////////////////////
936
937G4double G4PAIPhotData::GetEnergyPhotonTransfer(G4int coupleIndex,
938 std::size_t iPlace,
939 G4double position) const
940{
941 G4PhysicsVector* v = (*(fPAIphotonBank[coupleIndex]))(iPlace);
942 if(position*v->Energy(0) >= (*v)[0]) return v->Energy(0);
943
944 std::size_t iTransferMax = v->GetVectorLength() - 1;
945
946 std::size_t iTransfer;
947 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
948
949 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
950 {
951 x2 = v->Energy(iTransfer);
952 y2 = (*v)[iTransfer]/x2;
953 if(position >= y2) break;
954 }
955 x1 = v->Energy(iTransfer-1);
956 y1 = (*v)[iTransfer-1]/x1;
957
958 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
959 // << " x1= " << x1 << " x2= " << x2 << G4endl;
960
961 energyTransfer = x1;
962
963 if ( x1 != x2 )
964 {
965 if ( y1 == y2 )
966 {
967 energyTransfer += (x2 - x1)*G4UniformRand();
968 }
969 else
970 {
971 if( x1*1.1 < x2 )
972 {
973 const G4int nbins = 5;
974 G4double del = (x2 - x1)/G4int(nbins);
975 x2 = x1;
976
977 for(G4int i=1; i<=nbins; ++i)
978 {
979 x2 += del;
980 y2 = v->Value(x2)/x2;
981 if(position >= y2) { break; }
982 x1 = x2;
983 y1 = y2;
984 }
985 }
986 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
987 }
988 }
989 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
990 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
991 // << " E(keV)= " << energyTransfer/keV << G4endl;
992
993 return energyTransfer;
994}
995
996/////////////////////////////////////////////////////////////////////////
997
998G4double G4PAIPhotData::GetEnergyPlasmonTransfer(G4int coupleIndex,
999 std::size_t iPlace,
1000 G4double position) const
1001{
1002 G4PhysicsVector* v = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
1003
1004 if( position*v->Energy(0) >= (*v)[0] ) return v->Energy(0);
1005
1006 std::size_t iTransferMax = v->GetVectorLength() - 1;
1007
1008 std::size_t iTransfer;
1009 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1010
1011 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1012 {
1013 x2 = v->Energy(iTransfer);
1014 y2 = (*v)[iTransfer]/x2;
1015 if(position >= y2) break;
1016 }
1017 x1 = v->Energy(iTransfer-1);
1018 y1 = (*v)[iTransfer-1]/x1;
1019
1020 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
1021 // << " x1= " << x1 << " x2= " << x2 << G4endl;
1022
1023 energyTransfer = x1;
1024
1025 if ( x1 != x2 )
1026 {
1027 if ( y1 == y2 )
1028 {
1029 energyTransfer += (x2 - x1)*G4UniformRand();
1030 }
1031 else
1032 {
1033 if(x1*1.1 < x2)
1034 {
1035 const G4int nbins = 5;
1036 G4double del = (x2 - x1)/G4int(nbins);
1037 x2 = x1;
1038
1039 for( G4int i = 1; i <= nbins; ++i )
1040 {
1041 x2 += del;
1042 y2 = v->Value(x2)/x2;
1043
1044 if(position >= y2) break;
1045
1046 x1 = x2;
1047 y1 = y2;
1048 }
1049 }
1050 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1051 }
1052 }
1053 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
1054 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
1055 // << " E(keV)= " << energyTransfer/keV << G4endl;
1056
1057 return energyTransfer;
1058}
1059
1060//
1061//
1062//////////////////////////////////////////////////////////////////////
1063
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:67
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:85