| Geant4 Cross Reference |
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"
43 #include "G4PhysicalConstants.hh"
44
45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsFreeVector.hh"
47 #include "G4PhysicsTable.hh"
48 #include "G4MaterialCutsCouple.hh"
49 #include "G4ProductionCutsTable.hh"
50 #include "G4SandiaTable.hh"
51 #include "Randomize.hh"
52 #include "G4Poisson.hh"
53
54 ////////////////////////////////////////////////////////////////////////
55
56 using namespace std;
57
58 G4PAIPhotData::G4PAIPhotData(G4double tmin, G4double tmax, G4int ver)
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
93 G4PAIPhotData::~G4PAIPhotData()
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
126 void G4PAIPhotData::Initialise(const G4MaterialCutsCouple* couple,
127 G4double cut, G4PAIPhotModel* model)
128 {
129 G4ProductionCutsTable* theCoupleTable=
130 G4ProductionCutsTable::GetProductionCutsTable();
131 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
132 G4int jMatCC;
133
134 for (jMatCC = 0; jMatCC < numOfCouples; ++jMatCC )
135 {
136 if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
137 }
138 if( jMatCC == numOfCouples && jMatCC > 0 ) --jMatCC;
139
140 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
141 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
142 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
143 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
144 /*
145 G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
146 <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
147 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
148 */
149 // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??
150
151 auto dEdxCutVector =
152 new G4PhysicsLogVector(fLowestKineticEnergy,
153 fHighestKineticEnergy,
154 fTotBin);
155
156 auto dNdxCutVector =
157 new G4PhysicsLogVector(fLowestKineticEnergy,
158 fHighestKineticEnergy,
159 fTotBin);
160 auto dNdxCutPhotonVector =
161 new G4PhysicsLogVector(fLowestKineticEnergy,
162 fHighestKineticEnergy,
163 fTotBin);
164 auto dNdxCutPlasmonVector =
165 new G4PhysicsLogVector(fLowestKineticEnergy,
166 fHighestKineticEnergy,
167 fTotBin);
168
169 const G4Material* mat = couple->GetMaterial();
170 fSandia.Initialize(const_cast<G4Material*>(mat));
171
172 auto PAItransferTable = new G4PhysicsTable(fTotBin+1);
173 auto PAIphotonTable = new G4PhysicsTable(fTotBin+1);
174 auto PAIplasmonTable = new G4PhysicsTable(fTotBin+1);
175
176 auto PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
177 auto dEdxMeanVector =
178 new G4PhysicsLogVector(fLowestKineticEnergy,
179 fHighestKineticEnergy,
180 fTotBin);
181
182 // low energy Sandia interval
183 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
184
185 // energy safety
186 const G4double deltaLow = 100.*eV;
187
188 for (G4int i = 0; i <= fTotBin; ++i)
189 {
190 G4double kinEnergy = fParticleEnergyVector->Energy(i);
191 G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
192 G4double tau = kinEnergy/proton_mass_c2;
193 G4double bg2 = tau*( tau + 2. );
194
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
196
197 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
198
199 //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
200 // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
201
202 G4int n = fPAIxSection.GetSplineSize();
203
204 auto transferVector = new G4PhysicsFreeVector(n);
205 auto photonVector = new G4PhysicsFreeVector(n);
206 auto plasmonVector = new G4PhysicsFreeVector(n);
207
208 auto dEdxVector = new G4PhysicsFreeVector(n);
209
210 for( G4int k = 0; k < n; k++ )
211 {
212 G4double t = fPAIxSection.GetSplineEnergy(k+1);
213
214 transferVector->PutValue(k , t,
215 t*fPAIxSection.GetIntegralPAIxSection(k+1));
216 photonVector->PutValue(k , t,
217 t*fPAIxSection.GetIntegralCerenkov(k+1));
218 plasmonVector->PutValue(k , t,
219 t*fPAIxSection.GetIntegralPlasmon(k+1));
220
221 dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
222 }
223 // G4cout << *transferVector << G4endl;
224
225 G4double ionloss = std::max(fPAIxSection.GetMeanEnergyLoss(), 0.0);// total <dE/dx>
226 dEdxMeanVector->PutValue(i,ionloss);
227
228 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
229 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
230 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
231
232 G4double dEdxCut = dEdxVector->Value(cut)/cut;
233 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
234
235 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
236 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
237 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
238
239 dNdxCutVector->PutValue(i, dNdxCut);
240 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
241 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
242
243 dEdxCutVector->PutValue(i, dEdxCut);
244
245 PAItransferTable->insertAt(i,transferVector);
246 PAIphotonTable->insertAt(i,photonVector);
247 PAIplasmonTable->insertAt(i,plasmonVector);
248 PAIdEdxTable->insertAt(i,dEdxVector);
249
250 } // end of Tkin loop
251
252 fPAIxscBank.push_back(PAItransferTable);
253 fPAIphotonBank.push_back(PAIphotonTable);
254 fPAIplasmonBank.push_back(PAIplasmonTable);
255
256 fPAIdEdxBank.push_back(PAIdEdxTable);
257 fdEdxTable.push_back(dEdxMeanVector);
258
259 fdNdxCutTable.push_back(dNdxCutVector);
260 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
261 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
262
263 fdEdxCutTable.push_back(dEdxCutVector);
264 }
265
266 //////////////////////////////////////////////////////////////////////////////
267
268 G4double G4PAIPhotData::DEDXPerVolume(G4int coupleIndex, G4double scaledTkin,
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
303 G4double
304 G4PAIPhotData::CrossSectionPerVolume(G4int coupleIndex,
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
366 G4double
367 G4PAIPhotData::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
432 G4double G4PAIPhotData::SampleAlongStepTransfer(G4int coupleIndex,
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
524 G4double G4PAIPhotData::SampleAlongStepPhotonTransfer(G4int coupleIndex,
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
616 G4double G4PAIPhotData::SampleAlongStepPlasmonTransfer(G4int coupleIndex,
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
711 G4double G4PAIPhotData::SamplePostStepTransfer(G4int coupleIndex,
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
767 G4double G4PAIPhotData::SamplePostStepPhotonTransfer(G4int coupleIndex,
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
825 G4double G4PAIPhotData::SamplePostStepPlasmonTransfer(G4int coupleIndex,
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
886 G4double 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
937 G4double 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
998 G4double 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
1064