| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class
30 // File name: G4PAIPhotData.cc
31 //
32 // Author: V.Grichine based on G4PAIMod
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, G4
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, lowes
67 fHighestKineticEnergy = tmax;
68
69 if(tmax < 10*fLowestKineticEnergy)
70 {
71 fHighestKineticEnergy = 10*fLowestKineticE
72 }
73 else if(tmax > highestTkin)
74 {
75 fHighestKineticEnergy = std::max(highestTk
76 }
77 fTotBin = (G4int)(nPerDecade*
78 std::log10(fHighestKineticEnergy/fLowe
79
80 fParticleEnergyVector = new G4PhysicsLogVect
81 fHighestKineticEnergy,
82 fTotBin);
83 if(0 < ver) {
84 G4cout << "### G4PAIPhotData: Nbins= " <<
85 << " Tmin(MeV)= " << fLowestKineticEnergy
86 << " Tmax(GeV)= " << fHighestKineticEnerg
87 << " tmin(keV)= " << tmin/keV << G4endl;
88 }
89 }
90
91 //////////////////////////////////////////////
92
93 G4PAIPhotData::~G4PAIPhotData()
94 {
95 //G4cout << "G4PAIPhotData::~G4PAIPhotData()
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()
122 }
123
124 //////////////////////////////////////////////
125
126 void G4PAIPhotData::Initialise(const G4Materia
127 G4double cut,
128 {
129 G4ProductionCutsTable* theCoupleTable=
130 G4ProductionCutsTable::GetProductionCu
131 G4int numOfCouples = (G4int)theCoupleTable->
132 G4int jMatCC;
133
134 for (jMatCC = 0; jMatCC < numOfCouples; ++jM
135 {
136 if( couple == theCoupleTable->GetMaterialC
137 }
138 if( jMatCC == numOfCouples && jMatCC > 0 ) -
139
140 const vector<G4double>* deltaCutInKineticEn
141 const vector<G4double>* photonCutInKineticE
142 G4double deltaCutInKineticEnergyNow = (*del
143 G4double photonCutInKineticEnergyNow = (*pho
144 /*
145 G4cout<<"G4PAIPhotData::Initialise: "<<"cut
146 <<deltaCutInKineticEnergyNow/keV<<" ke
147 <<photonCutInKineticEnergyNow/keV<<" keV"<<G
148 */
149 // if( deltaCutInKineticEnergyNow != cut ) d
150
151 auto dEdxCutVector =
152 new G4PhysicsLogVector(fLowestKineticEnerg
153 fHighestKineticEnergy,
154 fTotBin);
155
156 auto dNdxCutVector =
157 new G4PhysicsLogVector(fLowestKineticEnerg
158 fHighestKineticEnergy,
159 fTotBin);
160 auto dNdxCutPhotonVector =
161 new G4PhysicsLogVector(fLowestKineticEnerg
162 fHighestKineticEnergy,
163 fTotBin);
164 auto dNdxCutPlasmonVector =
165 new G4PhysicsLogVector(fLowestKineticEnerg
166 fHighestKineticEnergy,
167 fTotBin);
168
169 const G4Material* mat = couple->GetMaterial(
170 fSandia.Initialize(const_cast<G4Material*>(m
171
172 auto PAItransferTable = new G4PhysicsTable(f
173 auto PAIphotonTable = new G4PhysicsTable(fTo
174 auto PAIplasmonTable = new G4PhysicsTable(fT
175
176 auto PAIdEdxTable = new G4PhysicsTable(fTotB
177 auto dEdxMeanVector =
178 new G4PhysicsLogVector(fLowestKineticEnerg
179 fHighestKineticEnergy,
180 fTotBin);
181
182 // low energy Sandia interval
183 G4double Tmin = fSandia.GetSandiaMatTablePAI
184
185 // energy safety
186 const G4double deltaLow = 100.*eV;
187
188 for (G4int i = 0; i <= fTotBin; ++i)
189 {
190 G4double kinEnergy = fParticleEnergyVector
191 G4double Tmax = model->ComputeMaxEnergy(ki
192 G4double tau = kinEnergy/proton_mass_c2;
193 G4double bg2 = tau*( tau + 2. );
194
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin
196
197 fPAIxSection.Initialize( mat, Tmax, bg2, &
198
199 //G4cout << i << ". TransferMax(keV)= "<<
200 // << cut/keV << " E(MeV)= " << kinEn
201
202 G4int n = fPAIxSection.GetSplineSize();
203
204 auto transferVector = new G4PhysicsFreeVec
205 auto photonVector = new G4PhysicsFreeVec
206 auto plasmonVector = new G4PhysicsFreeVec
207
208 auto dEdxVector = new G4PhysicsFreeVec
209
210 for( G4int k = 0; k < n; k++ )
211 {
212 G4double t = fPAIxSection.GetSplineEnerg
213
214 transferVector->PutValue(k , t,
215 t*fPAIxSection.
216 photonVector->PutValue(k , t,
217 t*fPAIxSection.
218 plasmonVector->PutValue(k , t,
219 t*fPAIxSection.
220
221 dEdxVector->PutValue(k, t, fPAIxSection.
222 }
223 // G4cout << *transferVector << G4endl;
224
225 G4double ionloss = std::max(fPAIxSection.G
226 dEdxMeanVector->PutValue(i,ionloss);
227
228 G4double dNdxCut = transferVector->Value(d
229 G4double dNdxCutPhoton = photonVector->Val
230 G4double dNdxCutPlasmon = plasmonVector->V
231
232 G4double dEdxCut = dEdxVector->Value(cut)/
233 //G4cout << "i= " << i << " x= " << dNdxCu
234
235 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
236 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton =
237 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon
238
239 dNdxCutVector->PutValue(i, dNdxCut);
240 dNdxCutPhotonVector->PutValue(i, dNdxCutPh
241 dNdxCutPlasmonVector->PutValue(i, dNdxCutP
242
243 dEdxCutVector->PutValue(i, dEdxCut);
244
245 PAItransferTable->insertAt(i,transferVecto
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(dNdxCutPhotonV
261 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmo
262
263 fdEdxCutTable.push_back(dEdxCutVector);
264 }
265
266 //////////////////////////////////////////////
267
268 G4double G4PAIPhotData::DEDXPerVolume(G4int co
269 G4double cut) const
270 {
271 // VI: iPlace is the low edge index of the b
272 // iPlace is in interval from 0 to (N-1)
273 std::size_t iPlace = fParticleEnergyVector->
274 std::size_t nPlace = fParticleEnergyVector->
275
276 G4bool one = true;
277 if(scaledTkin >= fParticleEnergyVector->Ener
278 else if(scaledTkin > fParticleEnergyVector->
279 one = false;
280 }
281
282 // VI: apply interpolation of the vector
283 G4double dEdx = fdEdxTable[coupleIndex]->Val
284 G4double del = (*(fPAIdEdxBank[coupleIndex]
285 if(!one) {
286 G4double del2 = (*(fPAIdEdxBank[coupleInde
287 G4double E1 = fParticleEnergyVector->Energ
288 G4double E2 = fParticleEnergyVector->Energ
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 cou
305 G4double scaledTkin,
306 G4double tcut, G4double tmax) co
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->
315 std::size_t nPlace = fParticleEnergyVector->
316
317 G4bool one = true;
318
319 if( scaledTkin >= fParticleEnergyVector
320 else if( scaledTkin > fParticleEnergyVector-
321
322
323 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex]
324 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])
325
326 xscPh = xscPh2;
327 xscEl = xscEl2;
328
329 cross = xscPh + xscEl;
330
331 if( !one )
332 {
333 xscEl2 = (*fdNdxCutPlasmonTable[coupleInde
334
335 G4double E1 = fParticleEnergyVector->Energ
336 G4double E2 = fParticleEnergyVector->Energ
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
346
347 E1 = fParticleEnergyVector->Energy(iPlace)
348 E2 = fParticleEnergyVector->Energy(iPlace+
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 coupleInd
368 {
369 G4double cross, xscEl, xscEl2, xscPh, xscPh2
370 // iPlace is in interval from 0 to (N-1)
371
372 std::size_t iPlace = fParticleEnergyVector->
373 std::size_t nPlace = fParticleEnergyVector->
374
375 G4bool one = true;
376
377 if( scaledTkin >= fParticleEnergyVector
378 else if( scaledTkin > fParticleEnergyVector-
379
380
381 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex]
382 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])
383
384 xscPh = xscPh2;
385 xscEl = xscEl2;
386
387 cross = xscPh + xscEl;
388
389 if( !one )
390 {
391 xscEl2 = (*fdNdxCutPlasmonTable[coupleInde
392
393 G4double E1 = fParticleEnergyVector->Energ
394 G4double E2 = fParticleEnergyVector->Energ
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
404
405 E1 = fParticleEnergyVector->Energy(iPlace)
406 E2 = fParticleEnergyVector->Energy(iPlace+
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
426 }
427 return plRatio;
428 }
429
430 //////////////////////////////////////////////
431
432 G4double G4PAIPhotData::SampleAlongStepTransfe
433
434 G4double scaledTkin,
435 G4double stepFactor) const
436 {
437 G4double loss = 0.0;
438 G4double omega;
439 G4double position, E1, E2, W1, W2, W, dNdxCu
440
441 std::size_t iPlace = fParticleEnergyVector->
442 std::size_t nPlace = fParticleEnergyVector->
443
444 G4bool one = true;
445
446 if (scaledTkin >= fParticleEnergyVector-
447 else if(scaledTkin > fParticleEnergyVector->
448
449 G4PhysicsLogVector* vcut = fdNdxCutTable[cou
450 G4PhysicsVector* v1 = (*(fPAIxscBank[co
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)*s
458
459 meanNumber = meanN1;
460
461 // G4cout<<"iPlace = "<<iPlace<< " meanN1= "
462 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< "
463
464 dNdxCut2 = dNdxCut1;
465 W1 = 1.0;
466 W2 = 0.0;
467 if(!one)
468 {
469 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+
470 dNdxCut2 = (*vcut)[iPlace+1];
471 e2 = v2->Energy(0);
472
473 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)
474
475 E1 = fParticleEnergyVector->Energy(iPlace)
476 E2 = fParticleEnergyVector->Energy(iPlace+
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 << " me
483 // << " dNdxCut2= " << dNdxCut2 << G4en
484 }
485 if( meanNumber <= 0.0) return 0.0;
486
487 G4int numOfCollisions = (G4int)G4Poisson(mea
488
489 //G4cout << "N= " << numOfCollisions << G4en
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 - dNdxC
497 omega = GetEnergyTransfer(coupleIndex, iPl
498
499 //G4cout << "omega(keV)= " << omega/keV <<
500
501 if(!one)
502 {
503 position = dNdxCut2 + ((*v2)[0]/e2 - dNd
504 G4double omega2 = GetEnergyTransfer(coup
505 omega = omega*W1 + omega2*W2;
506 }
507 //G4cout << "omega(keV)= " << omega/keV <<
508
509 loss += omega;
510 if( loss > kinEnergy) { break; }
511 }
512
513 // G4cout<<"PAIPhotData AlongStepLoss = "<<l
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::SampleAlongStepPhotonT
525
526 G4double scaledTkin,
527 G4double stepFactor) const
528 {
529 G4double loss = 0.0;
530 G4double omega;
531 G4double position, E1, E2, W1, W2, W, dNdxCu
532
533 std::size_t iPlace = fParticleEnergyVector->
534 std::size_t nPlace = fParticleEnergyVector->
535
536 G4bool one = true;
537
538 if (scaledTkin >= fParticleEnergyVector-
539 else if(scaledTkin > fParticleEnergyVector->
540
541 G4PhysicsLogVector* vcut = fdNdxCutPhotonTab
542 G4PhysicsVector* v1 = (*(fPAIphotonBank
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)*s
550
551 meanNumber = meanN1;
552
553 // G4cout<<"iPlace = "<<iPlace<< " meanN1= "
554 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< "
555
556 dNdxCut2 = dNdxCut1;
557 W1 = 1.0;
558 W2 = 0.0;
559 if(!one)
560 {
561 v2 = (*(fPAIphotonBank[coupleIndex]))(iPla
562 dNdxCut2 = (*vcut)[iPlace+1];
563 e2 = v2->Energy(0);
564
565 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)
566
567 E1 = fParticleEnergyVector->Energy(iPlace)
568 E2 = fParticleEnergyVector->Energy(iPlace+
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 << " me
575 // << " dNdxCut2= " << dNdxCut2 << G4en
576 }
577 if( meanNumber <= 0.0) return 0.0;
578
579 G4int numOfCollisions = (G4int)G4Poisson(mea
580
581 //G4cout << "N= " << numOfCollisions << G4en
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 - dNdxC
589 omega = GetEnergyPhotonTransfer(coupleInde
590
591 //G4cout << "omega(keV)= " << omega/keV <<
592
593 if(!one)
594 {
595 position = dNdxCut2 + ((*v2)[0]/e2 - dNd
596 G4double omega2 = GetEnergyPhotonTransfe
597 omega = omega*W1 + omega2*W2;
598 }
599 //G4cout << "omega(keV)= " << omega/keV <<
600
601 loss += omega;
602 if( loss > kinEnergy) { break; }
603 }
604
605 // G4cout<<"PAIPhotData AlongStepLoss = "<<l
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::SampleAlongStepPlasmon
617
618 G4double scaledTkin,
619 G4double stepFactor) const
620 {
621 G4double loss = 0.0;
622 G4double omega;
623 G4double position, E1, E2, W1, W2, W, dNdxCu
624
625 std::size_t iPlace = fParticleEnergyVector->
626 std::size_t nPlace = fParticleEnergyVector->
627
628 G4bool one = true;
629
630 if (scaledTkin >= fParticleEnergyVector-
631 else if(scaledTkin > fParticleEnergyVector->
632
633 G4PhysicsLogVector* vcut = fdNdxCutPlasmonTa
634 G4PhysicsVector* v1 = (*(fPAIplasmonBan
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)*s
642
643 meanNumber = meanN1;
644
645 // G4cout<<"iPlace = "<<iPlace<< " meanN1= "
646 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< "
647
648 dNdxCut2 = dNdxCut1;
649 W1 = 1.0;
650 W2 = 0.0;
651 if(!one)
652 {
653 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPl
654 dNdxCut2 = (*vcut)[iPlace+1];
655 e2 = v2->Energy(0);
656
657 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)
658
659 E1 = fParticleEnergyVector->Energy(iPlace)
660 E2 = fParticleEnergyVector->Energy(iPlace+
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 << " me
667 // << " dNdxCut2= " << dNdxCut2 << G4en
668 }
669 if( meanNumber <= 0.0) return 0.0;
670
671 G4int numOfCollisions = (G4int)G4Poisson(mea
672
673 //G4cout << "N= " << numOfCollisions << G4en
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 - dNdxC
681 omega = GetEnergyPlasmonTransfer(coupleInd
682
683 //G4cout << "omega(keV)= " << omega/keV <<
684
685 if(!one)
686 {
687 position = dNdxCut2 + ((*v2)[0]/e2 - dNd
688 G4double omega2 = GetEnergyPlasmonTransf
689 omega = omega*W1 + omega2*W2;
690 }
691 //G4cout << "omega(keV)= " << omega/keV <<
692
693 loss += omega;
694 if( loss > kinEnergy) { break; }
695 }
696
697 // G4cout<<"PAIPhotData AlongStepLoss = "<<l
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
709 // according to passed scaled kinetic energy o
710
711 G4double G4PAIPhotData::SamplePostStepTransfer
712 G4double scaledTkin) const
713 {
714 //G4cout<<"G4PAIPhotData::GetPostStepTransfe
715 G4double transfer = 0.0;
716 G4double rand = G4UniformRand();
717
718 std::size_t nPlace = fParticleEnergyVector->
719
720 // std::size_t iTransfer, iTr1, iTr2;
721 G4double position, dNdxCut1, dNdxCut2, E1, E
722
723 G4PhysicsVector* cutv = fdNdxCutTable[couple
724
725 // Fermi plato, try from left
726 if( scaledTkin >= fParticleEnergyVector->Get
727 {
728 position = (*cutv)[nPlace]*rand;
729 transfer = GetEnergyTransfer(coupleIndex,
730 }
731 else if( scaledTkin <= fParticleEnergyVector
732 {
733 position = (*cutv)[0]*rand;
734 transfer = GetEnergyTransfer(coupleIndex,
735 }
736 else
737 {
738 std::size_t iPlace = fParticleEnergyVector
739
740 dNdxCut1 = (*cutv)[iPlace];
741 dNdxCut2 = (*cutv)[iPlace+1];
742
743 E1 = fParticleEnergyVector->Energy(iPlace)
744 E2 = fParticleEnergyVector->Energy(iPlace+
745 W = 1.0/(E2 - E1);
746 W1 = (E2 - scaledTkin)*W;
747 W2 = (scaledTkin - E1)*W;
748
749 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d
750 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= "
751
752 position = dNdxCut1*rand;
753 G4double tr1 = GetEnergyTransfer(coupleInd
754
755 position = dNdxCut2*rand;
756 G4double tr2 = GetEnergyTransfer(coupleInd
757
758 transfer = tr1*W1 + tr2*W2;
759 }
760 //G4cout<<"PAImodel PostStepTransfer = "<<tr
761 if(transfer < 0.0 ) { transfer = 0.0; }
762 return transfer;
763 }
764
765 //////////////////////////////////////////////
766
767 G4double G4PAIPhotData::SamplePostStepPhotonTr
768 G4double scaledTkin) const
769 {
770 //G4cout<<"G4PAIPhotData::GetPostStepTransfe
771 G4double transfer = 0.0;
772 G4double rand = G4UniformRand();
773
774 std::size_t nPlace = fParticleEnergyVector->
775
776 // std::size_t iTransfer, iTr1, iTr2;
777 G4double position, dNdxCut1, dNdxCut2, E1, E
778
779 G4PhysicsVector* cutv = fdNdxCutPhotonTable[
780
781 // Fermi plato, try from left
782
783 if( scaledTkin >= fParticleEnergyVector->Get
784 {
785 position = (*cutv)[nPlace]*rand;
786 transfer = GetEnergyPhotonTransfer(coupleI
787 }
788 else if( scaledTkin <= fParticleEnergyVector
789 {
790 position = (*cutv)[0]*rand;
791 transfer = GetEnergyPhotonTransfer(coupleI
792 }
793 else
794 {
795 std::size_t iPlace = fParticleEnergyVector
796
797 dNdxCut1 = (*cutv)[iPlace];
798 dNdxCut2 = (*cutv)[iPlace+1];
799
800 E1 = fParticleEnergyVector->Energy(iPlace)
801 E2 = fParticleEnergyVector->Energy(iPlace+
802 W = 1.0/(E2 - E1);
803 W1 = (E2 - scaledTkin)*W;
804 W2 = (scaledTkin - E1)*W;
805
806 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d
807 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= "
808
809 position = dNdxCut1*rand;
810
811 G4double tr1 = GetEnergyPhotonTransfer(cou
812
813 position = dNdxCut2*rand;
814 G4double tr2 = GetEnergyPhotonTransfer(cou
815
816 transfer = tr1*W1 + tr2*W2;
817 }
818 //G4cout<<"PAImodel PostStepTransfer = "<<tr
819 if(transfer < 0.0 ) { transfer = 0.0; }
820 return transfer;
821 }
822
823 //////////////////////////////////////////////
824
825 G4double G4PAIPhotData::SamplePostStepPlasmonT
826 G4double scaledTkin) const
827 {
828 //G4cout<<"G4PAIPhotData::GetPostStepTransfe
829 G4double transfer = 0.0;
830 G4double rand = G4UniformRand();
831
832 std::size_t nPlace = fParticleEnergyVector->
833
834 // std::size_t iTransfer, iTr1, iTr2;
835 G4double position, dNdxCut1, dNdxCut2, E1, E
836
837 G4PhysicsVector* cutv = fdNdxCutPlasmonTable
838
839 // Fermi plato, try from left
840 if( scaledTkin >= fParticleEnergyVector->Get
841 {
842 position = (*cutv)[nPlace]*rand;
843 transfer = GetEnergyPlasmonTransfer(couple
844 }
845 else if( scaledTkin <= fParticleEnergyVector
846 {
847 position = (*cutv)[0]*rand;
848 transfer = GetEnergyPlasmonTransfer(couple
849 }
850 else
851 {
852 std::size_t iPlace = fParticleEnergyVector
853
854 dNdxCut1 = (*cutv)[iPlace];
855 dNdxCut2 = (*cutv)[iPlace+1];
856
857 E1 = fParticleEnergyVector->Energy(iPlace)
858 E2 = fParticleEnergyVector->Energy(iPlace+
859 W = 1.0/(E2 - E1);
860 W1 = (E2 - scaledTkin)*W;
861 W2 = (scaledTkin - E1)*W;
862
863 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d
864 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= "
865
866 position = dNdxCut1*rand;
867 G4double tr1 = GetEnergyPlasmonTransfer(co
868
869 position = dNdxCut2*rand;
870 G4double tr2 = GetEnergyPlasmonTransfer(co
871
872 transfer = tr1*W1 + tr2*W2;
873 }
874 //G4cout<<"PAImodel PostStepPlasmonTransfer
875
876 if(transfer < 0.0 ) transfer = 0.0;
877
878 return transfer;
879 }
880
881 //////////////////////////////////////////////
882 //
883 // Returns PAI energy transfer according to pa
884 // indexes of particle kinetic enegry and rand
885
886 G4double G4PAIPhotData::GetEnergyTransfer(G4in
887 std::size_t iPlace,
888 G4double position) const
889 {
890 G4PhysicsVector* v = (*(fPAIxscBank[coupleIn
891 if(position*v->Energy(0) >= (*v)[0]) { retur
892
893 std::size_t iTransferMax = v->GetVectorLengt
894
895 std::size_t iTransfer;
896 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0),
897
898 for(iTransfer=1; iTransfer<=iTransferMax; ++
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= "
907 // << " x1= " << x1 << " x2= " << x2 << G4
908
909 energyTransfer = x1;
910 if ( x1 != x2 ) {
911 if ( y1 == y2 ) {
912 energyTransfer += (x2 - x1)*G4UniformRan
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/(positi
927 }
928 }
929 // G4cout << "x1(keV)= " << x1/keV << " x2(
930 // << " y1= " << y1 << " y2= " << y2 << "
931 // << " E(keV)= " << energyTransfer/keV <<
932 return energyTransfer;
933 }
934
935 //////////////////////////////////////////////
936
937 G4double G4PAIPhotData::GetEnergyPhotonTransfe
938 std::size_t iPlace,
939 G4double position) const
940 {
941 G4PhysicsVector* v = (*(fPAIphotonBank[coupl
942 if(position*v->Energy(0) >= (*v)[0]) return
943
944 std::size_t iTransferMax = v->GetVectorLengt
945
946 std::size_t iTransfer;
947 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0),
948
949 for(iTransfer=1; iTransfer<=iTransferMax; ++
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= "
959 // << " x1= " << x1 << " x2= " << x2 << G4
960
961 energyTransfer = x1;
962
963 if ( x1 != x2 )
964 {
965 if ( y1 == y2 )
966 {
967 energyTransfer += (x2 - x1)*G4UniformRan
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/(positi
987 }
988 }
989 // G4cout << "x1(keV)= " << x1/keV << " x2(
990 // << " y1= " << y1 << " y2= " << y2 << "
991 // << " E(keV)= " << energyTransfer/keV <<
992
993 return energyTransfer;
994 }
995
996 //////////////////////////////////////////////
997
998 G4double G4PAIPhotData::GetEnergyPlasmonTransf
999 std::size_t iPlace,
1000 G4double position) const
1001 {
1002 G4PhysicsVector* v = (*(fPAIplasmonBank[cou
1003
1004 if( position*v->Energy(0) >= (*v)[0] ) ret
1005
1006 std::size_t iTransferMax = v->GetVectorLeng
1007
1008 std::size_t iTransfer;
1009 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0)
1010
1011 for(iTransfer = 1; iTransfer <= iTransferMa
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= "
1021 // << " x1= " << x1 << " x2= " << x2 << G
1022
1023 energyTransfer = x1;
1024
1025 if ( x1 != x2 )
1026 {
1027 if ( y1 == y2 )
1028 {
1029 energyTransfer += (x2 - x1)*G4UniformRa
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/(posit
1051 }
1052 }
1053 // G4cout << "x1(keV)= " << x1/keV << " x2
1054 // << " y1= " << y1 << " y2= " << y2 << "
1055 // << " E(keV)= " << energyTransfer/keV <
1056
1057 return energyTransfer;
1058 }
1059
1060 //
1061 //
1062 /////////////////////////////////////////////
1063
1064