Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PAIySection.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// G4PAIySection.cc -- class implementation file
29//
30// GEANT 4 class implementation file
31//
32// For information related to this code, please, contact
33// the Geant4 Collaboration.
34//
35// R&D: [email protected]
36//
37// History:
38//
39// 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class
40// 26.07.09, V.Ivanchenko added protection for mumerical exceptions for
41// low-density materials
42// 21.11.10 V. Grichine bug fixed in Initialise for reading sandia table from
43// material. Warning: the table is tuned for photo-effect not PAI model.
44// 23.06.13 V.Grichine arrays->G4DataVectors
45//
46
47#include "G4PAIySection.hh"
48
49#include "globals.hh"
50#include "G4PhysicalConstants.hh"
51#include "G4SystemOfUnits.hh"
52#include "G4ios.hh"
53#include "G4Poisson.hh"
54#include "G4Material.hh"
55#include "G4MaterialCutsCouple.hh"
56#include "G4SandiaTable.hh"
57#include "G4Exp.hh"
58#include "G4Log.hh"
59
60using namespace std;
61
62// Local class constants
63
64const G4double G4PAIySection::fDelta = 0.005; // energy shift from interval border
65const G4double G4PAIySection::fError = 0.005; // error in lin-log approximation
66
67const G4int G4PAIySection::fMaxSplineSize = 500; // Max size of output spline
68 // arrays
69
70//////////////////////////////////////////////////////////////////
71//
72// Constructor
73//
74
75G4PAIySection::G4PAIySection()
76{
77 fSandia = nullptr;
78 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
79 fIntervalNumber = fSplineNumber = 0;
80 fVerbose = 0;
81
82 betaBohr = fine_structure_const;
83 G4double cofBetaBohr = 4.0;
84 G4double betaBohr2 = fine_structure_const*fine_structure_const;
85 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
86
87 fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
88 fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
89 fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
90 fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
91 fDifPAIySection = G4DataVector(fMaxSplineSize,0.0);
92 fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
93 fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
94 fIntegralPAIySection = G4DataVector(fMaxSplineSize,0.0);
95 fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
96 fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
97 fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
98
99 for( G4int i = 0; i < 500; ++i )
100 {
101 for( G4int j = 0; j < 112; ++j ) { fPAItable[i][j] = 0.0; }
102 }
103}
104
105////////////////////////////////////////////////////////////////////////////
106//
107//
108
109G4double G4PAIySection::GetLorentzFactor(G4int j) const
110{
111 return fLorentzFactor[j];
112}
113
114////////////////////////////////////////////////////////////////////////
115//
116// Constructor with beta*gamma square value called from G4PAIModel
117
118void G4PAIySection::Initialize( const G4Material* material,
119 G4double maxEnergyTransfer,
120 G4double betaGammaSq,
121 G4SandiaTable* sandia)
122{
123 if(fVerbose > 0)
124 {
125 G4cout<<G4endl;
126 G4cout<<"G4PAIySection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
127 G4cout<<G4endl;
128 }
129 G4int i, j;
130
131 fSandia = sandia;
132 fIntervalNumber = sandia->GetMaxInterval();
133 fDensity = material->GetDensity();
134 fElectronDensity = material->GetElectronDensity();
135
136 // fIntervalNumber--;
137
138 if( fVerbose > 0 )
139 {
140 G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "
141 <<fIntervalNumber<< " (beta*gamma)^2= " << betaGammaSq << G4endl;
142 }
143 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
144 fA1 = G4DataVector(fIntervalNumber+2,0.0);
145 fA2 = G4DataVector(fIntervalNumber+2,0.0);
146 fA3 = G4DataVector(fIntervalNumber+2,0.0);
147 fA4 = G4DataVector(fIntervalNumber+2,0.0);
148
149 for( i = 1; i <= fIntervalNumber; ++i )
150 {
151 if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV )
152 {
153 fIntervalNumber--;
154 continue;
155 }
156 if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer )
157 || i >= fIntervalNumber )
158 {
159 fEnergyInterval[i] = maxEnergyTransfer;
160 fIntervalNumber = i;
161 break;
162 }
163 fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
164 fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
165 fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
166 fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
167 fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
168
169 if( fVerbose > 0 ) {
170 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
171 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
172 }
173 }
174 if( fVerbose > 0 ) {
175 G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "
176 <<fIntervalNumber<<G4endl;
177 }
178 if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
179 {
180 fIntervalNumber++;
181 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
182 }
183 if( fVerbose > 0 )
184 {
185 for( i = 1; i <= fIntervalNumber; ++i )
186 {
187 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
188 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
189 }
190 }
191 if( fVerbose > 0 ) {
192 G4cout<<"Now checking, if two borders are too close together"<<G4endl;
193 }
194 for( i = 1; i < fIntervalNumber; ++i )
195 {
196 if( fEnergyInterval[i+1]-fEnergyInterval[i] >
197 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) ) continue;
198 else
199 {
200 for( j = i; j < fIntervalNumber; j++ )
201 {
202 fEnergyInterval[j] = fEnergyInterval[j+1];
203 fA1[j] = fA1[j+1];
204 fA2[j] = fA2[j+1];
205 fA3[j] = fA3[j+1];
206 fA4[j] = fA4[j+1];
207 }
208 fIntervalNumber--;
209 }
210 }
211 if( fVerbose > 0 )
212 {
213 for( i = 1; i <= fIntervalNumber; ++i )
214 {
215 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
216 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
217 }
218 }
219 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
220
221 ComputeLowEnergyCof(material);
222
223 G4double betaGammaSqRef =
224 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
225
226 NormShift(betaGammaSqRef);
227 SplainPAI(betaGammaSqRef);
228
229 // Preparation of integral PAI cross section for input betaGammaSq
230
231 for( i = 1; i <= fSplineNumber; ++i )
232 {
233 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
234
235 if( fVerbose > 0 ) G4cout<<i<<"; dNdxPAI = "<<fDifPAIySection[i]<<G4endl;
236 }
237 IntegralPAIySection();
238}
239
240/////////////////////////////////////////////////////////////////////////
241//
242// Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4.
243//
244
245void G4PAIySection::ComputeLowEnergyCof(const G4Material* material)
246{
247 G4int i, numberOfElements = (G4int)material->GetNumberOfElements();
248 G4double sumZ = 0., sumCof = 0.;
249
250 static const G4double p0 = 1.20923e+00;
251 static const G4double p1 = 3.53256e-01;
252 static const G4double p2 = -1.45052e-03;
253
254 G4double* thisMaterialZ = new G4double[numberOfElements];
255 G4double* thisMaterialCof = new G4double[numberOfElements];
256
257 for( i = 0; i < numberOfElements; ++i )
258 {
259 thisMaterialZ[i] = material->GetElement(i)->GetZ();
260 sumZ += thisMaterialZ[i];
261 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
262 }
263 for( i = 0; i < numberOfElements; ++i )
264 {
265 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
266 }
267 fLowEnergyCof = sumCof;
268 delete [] thisMaterialZ;
269 delete [] thisMaterialCof;
270 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
271}
272
273/////////////////////////////////////////////////////////////////////////
274//
275// General control function for class G4PAIySection
276//
277
278void G4PAIySection::InitPAI()
279{
280 G4int i;
281 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
282 fLorentzFactor[fRefGammaNumber] - 1;
283
284 // Preparation of integral PAI cross section for reference gamma
285
286 NormShift(betaGammaSq);
287 SplainPAI(betaGammaSq);
288
289 IntegralPAIySection();
290 IntegralCerenkov();
291 IntegralPlasmon();
292
293 for( i = 0; i<= fSplineNumber; ++i)
294 {
295 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i];
296
297 if(i != 0) fPAItable[i][0] = fSplineEnergy[i];
298 }
299 fPAItable[0][0] = fSplineNumber;
300
301 for( G4int j = 1; j < 112; ++j) // for other gammas
302 {
303 if( j == fRefGammaNumber ) continue;
304
305 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
306
307 for(i = 1; i <= fSplineNumber; ++i)
308 {
309 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
310 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
311 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
312 }
313 IntegralPAIySection();
314 IntegralCerenkov();
315 IntegralPlasmon();
316
317 for(i = 0; i <= fSplineNumber; ++i)
318 {
319 fPAItable[i][j] = fIntegralPAIySection[i];
320 }
321 }
322}
323
324///////////////////////////////////////////////////////////////////////
325//
326// Shifting from borders to intervals Creation of first energy points
327//
328
329void G4PAIySection::NormShift(G4double betaGammaSq)
330{
331 G4int i, j;
332
333 for( i = 1; i <= fIntervalNumber-1; ++i)
334 {
335 for( j = 1; j <= 2; ++j)
336 {
337 fSplineNumber = (i-1)*2 + j;
338
339 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
340 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
341 // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
342 // <<fSplineEnergy[fSplineNumber]<<G4endl;
343 }
344 }
345 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
346
347 j = 1;
348
349 for(i=2;i<=fSplineNumber;++i)
350 {
351 if(fSplineEnergy[i]<fEnergyInterval[j+1])
352 {
353 fIntegralTerm[i] = fIntegralTerm[i-1] +
354 RutherfordIntegral(j,fSplineEnergy[i-1],
355 fSplineEnergy[i] );
356 }
357 else
358 {
359 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
360 fEnergyInterval[j+1] );
361 j++;
362 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
363 RutherfordIntegral(j,fEnergyInterval[j],
364 fSplineEnergy[i] );
365 }
366 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
367 }
368 static const G4double nfactor =
369 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
370 fNormalizationCof = nfactor*fElectronDensity/fIntegralTerm[fSplineNumber];
371
372 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
373
374 // Calculation of PAI differrential cross-section (1/(keV*cm))
375 // in the energy points near borders of energy intervals
376
377 for(G4int k=1; k<=fIntervalNumber-1; ++k)
378 {
379 for(j=1; j<=2; ++j)
380 {
381 i = (k-1)*2 + j;
382 fImPartDielectricConst[i] = fNormalizationCof*
383 ImPartDielectricConst(k,fSplineEnergy[i]);
384 fRePartDielectricConst[i] = fNormalizationCof*
385 RePartDielectricConst(fSplineEnergy[i]);
386 fIntegralTerm[i] *= fNormalizationCof;
387
388 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
389 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
390 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
391 }
392 }
393
394} // end of NormShift
395
396/////////////////////////////////////////////////////////////////////////
397//
398// Creation of new energy points as geometrical mean of existing
399// one, calculation PAI_cs for them, while the error of logarithmic
400// linear approximation would be smaller than 'fError'
401
402void G4PAIySection::SplainPAI(G4double betaGammaSq)
403{
404 G4int k = 1;
405 G4int i = 1;
406
407 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
408 {
409 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
410 {
411 k++; // Here next energy point is in next energy interval
412 ++i;
413 continue;
414 }
415 // Shifting of arrayes for inserting the geometrical
416 // average of 'i' and 'i+1' energy points to 'i+1' place
417 fSplineNumber++;
418
419 for(G4int j = fSplineNumber; j >= i+2; j-- )
420 {
421 fSplineEnergy[j] = fSplineEnergy[j-1];
422 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
423 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
424 fIntegralTerm[j] = fIntegralTerm[j-1];
425
426 fDifPAIySection[j] = fDifPAIySection[j-1];
427 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
428 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
429 }
430 G4double x1 = fSplineEnergy[i];
431 G4double x2 = fSplineEnergy[i+1];
432 G4double yy1 = fDifPAIySection[i];
433 G4double y2 = fDifPAIySection[i+1];
434
435 G4double en1 = sqrt(x1*x2);
436 fSplineEnergy[i+1] = en1;
437
438 // Calculation of logarithmic linear approximation
439 // in this (enr) energy point, which number is 'i+1' now
440
441 G4double a = log10(y2/yy1)/log10(x2/x1);
442 G4double b = log10(yy1) - a*log10(x1);
443 G4double y = a*log10(en1) + b;
444 y = pow(10.,y);
445
446 // Calculation of the PAI dif. cross-section at this point
447
448 fImPartDielectricConst[i+1] = fNormalizationCof*
449 ImPartDielectricConst(k,fSplineEnergy[i+1]);
450 fRePartDielectricConst[i+1] = fNormalizationCof*
451 RePartDielectricConst(fSplineEnergy[i+1]);
452 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
453 RutherfordIntegral(k,fSplineEnergy[i],
454 fSplineEnergy[i+1]);
455
456 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
457 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
458 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
459
460 // Condition for next division of this segment or to pass
461 // to higher energies
462
463 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
464
465 G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])
466 /(fSplineEnergy[i+1]+fSplineEnergy[i]);
467
468 if( x < 0 )
469 {
470 x = -x;
471 }
472 if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
473 {
474 continue; // next division
475 }
476 i += 2; // pass to next segment
477
478 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
479 } // close 'while'
480
481} // end of SplainPAI
482
483
484////////////////////////////////////////////////////////////////////
485//
486// Integration over electrons that could be considered
487// quasi-free at energy transfer of interest
488
489G4double G4PAIySection::RutherfordIntegral( G4int k,
490 G4double x1,
491 G4double x2 )
492{
493 G4double c1, c2, c3;
494 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
495 G4double x12 = x1*x2;
496 c1 = (x2 - x1)/x12;
497 c2 = (x2 - x1)*(x2 + x1)/(x12*x12);
498 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
499 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
500
501 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
502
503} // end of RutherfordIntegral
504
505
506/////////////////////////////////////////////////////////////////
507//
508// Imaginary part of dielectric constant
509// (G4int k - interval number, G4double en1 - energy point)
510
511G4double G4PAIySection::ImPartDielectricConst( G4int k ,
512 G4double energy1 )
513{
514 G4double energy2,energy3,energy4,result;
515
516 energy2 = energy1*energy1;
517 energy3 = energy2*energy1;
518 energy4 = energy3*energy1;
519
520 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
521 result *=hbarc/energy1;
522
523 return result;
524
525} // end of ImPartDielectricConst
526
527
528//////////////////////////////////////////////////////////////////////////////
529//
530// Real part of dielectric constant minus unit: epsilon_1 - 1
531// (G4double enb - energy point)
532//
533
534G4double G4PAIySection::RePartDielectricConst(G4double enb)
535{
536 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
537 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
538
539 x0 = enb;
540 result = 0;
541
542 for(G4int i=1;i<=fIntervalNumber-1;++i)
543 {
544 x1 = fEnergyInterval[i];
545 x2 = fEnergyInterval[i+1];
546 xx1 = x1 - x0;
547 xx2 = x2 - x0;
548 xx12 = xx2/xx1;
549
550 if(xx12<0)
551 {
552 xx12 = -xx12;
553 }
554 xln1 = log(x2/x1);
555 xln2 = log(xx12);
556 xln3 = log((x2 + x0)/(x1 + x0));
557 x02 = x0*x0;
558 x03 = x02*x0;
559 x04 = x03*x0;
560 x05 = x04*x0;
561 G4double x12 = x1*x2;
562 c1 = (x2 - x1)/x12;
563 c2 = (x2 - x1)*(x2 +x1)/(x12*x12);
564 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
565
566 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
567 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
568 result -= fA3[i]*c2/2/x02;
569 result -= fA4[i]*c3/3/x02;
570
571 cof1 = fA1[i]/x02 + fA3[i]/x04;
572 cof2 = fA2[i]/x03 + fA4[i]/x05;
573
574 result += 0.5*(cof1 +cof2)*xln2;
575 result += 0.5*(cof1 - cof2)*xln3;
576 }
577 result *= 2*hbarc/pi;
578
579 return result;
580
581} // end of RePartDielectricConst
582
583//////////////////////////////////////////////////////////////////////
584//
585// PAI differential cross-section in terms of
586// simplified Allison's equation
587//
588
589G4double G4PAIySection::DifPAIySection( G4int i ,
590 G4double betaGammaSq )
591{
592 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
593 //G4double beta, be4;
594 //G4double be4;
595 // G4double betaBohr2 = fine_structure_const*fine_structure_const;
596 // G4double betaBohr4 = betaBohr2*betaBohr2*4.0;
597 be2 = betaGammaSq/(1 + betaGammaSq);
598 //be4 = be2*be2;
599 beta = sqrt(be2);
600 cof = 1;
601 x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
602
603 if( betaGammaSq < 0.01 ) x2 = log(be2);
604 else
605 {
606 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
607 (1/betaGammaSq - fRePartDielectricConst[i]) +
608 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
609 }
610 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
611 {
612 x6=0;
613 }
614 else
615 {
616 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
617 x5 = -1 - fRePartDielectricConst[i] +
618 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
619 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
620
621 x7 = atan2(fImPartDielectricConst[i],x3);
622 x6 = x5 * x7;
623 }
624 // if(fImPartDielectricConst[i] == 0) x6 = 0;
625
626 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
627 // if( x4 < 0.0 ) x4 = 0.0;
628 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
629 fImPartDielectricConst[i]*fImPartDielectricConst[i];
630
631 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
632 result = std::max(result, 1.0e-8);
633 result *= fine_structure_const/be2/pi;
634 // low energy correction
635
636 G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
637
638 result *= (1 - exp(-beta/betaBohr/lowCof));
639
640 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr));
641 // result *= (1-exp(-be2/betaBohr2));
642 // result *= (1-exp(-be4/betaBohr4));
643 // if(fDensity >= 0.1)
644 if(x8 > 0.)
645 {
646 result /= x8;
647 }
648 return result;
649
650} // end of DifPAIySection
651
652//////////////////////////////////////////////////////////////////////////
653//
654// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
655
656G4double G4PAIySection::PAIdNdxCerenkov( G4int i ,
657 G4double betaGammaSq )
658{
659 G4double logarithm, x3, x5, argument, modul2, dNdxC;
660 G4double be2, be4;
661
662 //G4double cof = 1.0;
663
664 be2 = betaGammaSq/(1 + betaGammaSq);
665 be4 = be2*be2;
666
667 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
668 else
669 {
670 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
671 (1/betaGammaSq - fRePartDielectricConst[i]) +
672 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
673 logarithm += log(1+1.0/betaGammaSq);
674 }
675
676 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
677 {
678 argument = 0.0;
679 }
680 else
681 {
682 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
683 x5 = -1.0 - fRePartDielectricConst[i] +
684 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
685 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
686 if( x3 == 0.0 ) argument = 0.5*pi;
687 else argument = atan2(fImPartDielectricConst[i],x3);
688 argument *= x5 ;
689 }
690 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
691
692 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
693
694 dNdxC *= fine_structure_const/be2/pi;
695
696 dNdxC *= (1-exp(-be4/betaBohr4));
697
698 // if(fDensity >= 0.1)
699 // {
700 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
701 fImPartDielectricConst[i]*fImPartDielectricConst[i];
702 if(modul2 > 0.)
703 {
704 dNdxC /= modul2;
705 }
706 return dNdxC;
707
708} // end of PAIdNdxCerenkov
709
710//////////////////////////////////////////////////////////////////////////
711//
712// Calculation od dN/dx of collisions with creation of longitudinal EM
713// excitations (plasmons, delta-electrons)
714
715G4double G4PAIySection::PAIdNdxPlasmon( G4int i ,
716 G4double betaGammaSq )
717{
718 G4double cof, resonance, modul2, dNdxP;
719 G4double be2, be4;
720
721 cof = 1;
722
723 be2 = betaGammaSq/(1 + betaGammaSq);
724 be4 = be2*be2;
725
726 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
727 resonance *= fImPartDielectricConst[i]/hbarc;
728
729 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
730
731 dNdxP = std::max(dNdxP, 1.0e-8);
732
733 dNdxP *= fine_structure_const/be2/pi;
734 dNdxP *= (1-exp(-be4/betaBohr4));
735
736// if( fDensity >= 0.1 )
737// {
738 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
739 fImPartDielectricConst[i]*fImPartDielectricConst[i];
740 if(modul2 > 0.)
741 {
742 dNdxP /= modul2;
743 }
744 return dNdxP;
745
746} // end of PAIdNdxPlasmon
747
748////////////////////////////////////////////////////////////////////////
749//
750// Calculation of the PAI integral cross-section
751// fIntegralPAIySection[1] = specific primary ionisation, 1/cm
752// and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm
753
754void G4PAIySection::IntegralPAIySection()
755{
756 fIntegralPAIySection[fSplineNumber] = 0;
757 fIntegralPAIdEdx[fSplineNumber] = 0;
758 fIntegralPAIySection[0] = 0;
759 G4int k = fIntervalNumber -1;
760
761 for(G4int i = fSplineNumber-1; i >= 1; i--)
762 {
763 if(fSplineEnergy[i] >= fEnergyInterval[k])
764 {
765 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
766 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
767 }
768 else
769 {
770 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
771 SumOverBorder(i+1,fEnergyInterval[k]);
772 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
773 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
774 k--;
775 }
776 }
777} // end of IntegralPAIySection
778
779////////////////////////////////////////////////////////////////////////
780//
781// Calculation of the PAI Cerenkov integral cross-section
782// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
783// and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
784
785void G4PAIySection::IntegralCerenkov()
786{
787 G4int i, k;
788 fIntegralCerenkov[fSplineNumber] = 0;
789 fIntegralCerenkov[0] = 0;
790 k = fIntervalNumber -1;
791
792 for( i = fSplineNumber-1; i >= 1; i-- )
793 {
794 if(fSplineEnergy[i] >= fEnergyInterval[k])
795 {
796 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
797 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
798 }
799 else
800 {
801 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
802 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
803 k--;
804 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
805 }
806 }
807
808} // end of IntegralCerenkov
809
810////////////////////////////////////////////////////////////////////////
811//
812// Calculation of the PAI Plasmon integral cross-section
813// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
814// and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
815
816void G4PAIySection::IntegralPlasmon()
817{
818 fIntegralPlasmon[fSplineNumber] = 0;
819 fIntegralPlasmon[0] = 0;
820 G4int k = fIntervalNumber -1;
821 for(G4int i=fSplineNumber-1;i>=1;i--)
822 {
823 if(fSplineEnergy[i] >= fEnergyInterval[k])
824 {
825 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
826 }
827 else
828 {
829 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
830 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
831 k--;
832 }
833 }
834
835} // end of IntegralPlasmon
836
837//////////////////////////////////////////////////////////////////////
838//
839// Calculation the PAI integral cross-section inside
840// of interval of continuous values of photo-ionisation
841// cross-section. Parameter 'i' is the number of interval.
842
843G4double G4PAIySection::SumOverInterval( G4int i )
844{
845 G4double x0,x1,y0,yy1,a,b,c,result;
846
847 x0 = fSplineEnergy[i];
848 x1 = fSplineEnergy[i+1];
849
850 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
851
852 y0 = fDifPAIySection[i];
853 yy1 = fDifPAIySection[i+1];
854 //G4cout << "## x0= " << x0 << " x1= " << x1 << G4endl;
855 c = x1/x0;
856 //G4cout << "c= " << c << " y0= " << y0 << " yy1= " << yy1 << G4endl;
857 a = log10(yy1/y0)/log10(c);
858 //G4cout << "a= " << a << G4endl;
859 // b = log10(y0) - a*log10(x0);
860 b = y0/pow(x0,a);
861 a += 1;
862 if(a == 0)
863 {
864 result = b*log(x1/x0);
865 }
866 else
867 {
868 result = y0*(x1*pow(c,a-1) - x0)/a;
869 }
870 a++;
871 if(a == 0)
872 {
873 fIntegralPAIySection[0] += b*log(x1/x0);
874 }
875 else
876 {
877 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
878 }
879 return result;
880
881} // end of SumOverInterval
882
883/////////////////////////////////
884
885G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
886{
887 G4double x0,x1,y0,yy1,a,b,c,result;
888
889 x0 = fSplineEnergy[i];
890 x1 = fSplineEnergy[i+1];
891
892 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
893
894 y0 = fDifPAIySection[i];
895 yy1 = fDifPAIySection[i+1];
896 c = x1/x0;
897 a = log10(yy1/y0)/log10(c);
898 // b = log10(y0) - a*log10(x0);
899 b = y0/pow(x0,a);
900 a += 2;
901 if(a == 0)
902 {
903 result = b*log(x1/x0);
904 }
905 else
906 {
907 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
908 }
909 return result;
910
911} // end of SumOverInterval
912
913//////////////////////////////////////////////////////////////////////
914//
915// Calculation the PAI Cerenkov integral cross-section inside
916// of interval of continuous values of photo-ionisation Cerenkov
917// cross-section. Parameter 'i' is the number of interval.
918
919G4double G4PAIySection::SumOverInterCerenkov( G4int i )
920{
921 G4double x0,x1,y0,yy1,a,c,result;
922
923 x0 = fSplineEnergy[i];
924 x1 = fSplineEnergy[i+1];
925
926 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
927
928 y0 = fdNdxCerenkov[i];
929 yy1 = fdNdxCerenkov[i+1];
930 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
931 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
932
933 c = x1/x0;
934 a = log10(yy1/y0)/log10(c);
935 G4double b = 0.0;
936 if(a < 20.) b = y0/pow(x0,a);
937
938 a += 1.0;
939 if(a == 0) result = b*log(c);
940 else result = y0*(x1*pow(c,a-1) - x0)/a;
941 a += 1.0;
942
943 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
944 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
945 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
946 return result;
947
948} // end of SumOverInterCerenkov
949
950//////////////////////////////////////////////////////////////////////
951//
952// Calculation the PAI Plasmon integral cross-section inside
953// of interval of continuous values of photo-ionisation Plasmon
954// cross-section. Parameter 'i' is the number of interval.
955
956G4double G4PAIySection::SumOverInterPlasmon( G4int i )
957{
958 G4double x0,x1,y0,yy1,a,c,result;
959
960 x0 = fSplineEnergy[i];
961 x1 = fSplineEnergy[i+1];
962
963 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
964
965 y0 = fdNdxPlasmon[i];
966 yy1 = fdNdxPlasmon[i+1];
967 c = x1/x0;
968 a = log10(yy1/y0)/log10(c);
969
970 G4double b = 0.0;
971 if(a < 20.) b = y0/pow(x0,a);
972
973 a += 1.0;
974 if(a == 0) result = b*log(x1/x0);
975 else result = y0*(x1*pow(c,a-1) - x0)/a;
976 a += 1.0;
977
978 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
979 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
980
981 return result;
982
983} // end of SumOverInterPlasmon
984
985///////////////////////////////////////////////////////////////////////////////
986//
987// Integration of PAI cross-section for the case of
988// passing across border between intervals
989
990G4double G4PAIySection::SumOverBorder( G4int i ,
991 G4double en0 )
992{
993 G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
994
995 e0 = en0;
996 x0 = fSplineEnergy[i];
997 x1 = fSplineEnergy[i+1];
998 y0 = fDifPAIySection[i];
999 yy1 = fDifPAIySection[i+1];
1000
1001 //c = x1/x0;
1002 d = e0/x0;
1003 a = log10(yy1/y0)/log10(x1/x0);
1004
1005 G4double b = 0.0;
1006 if(a < 20.) b = y0/pow(x0,a);
1007
1008 a += 1;
1009 if(a == 0)
1010 {
1011 result = b*log(x0/e0);
1012 }
1013 else
1014 {
1015 result = y0*(x0 - e0*pow(d,a-1))/a;
1016 }
1017 a++;
1018 if(a == 0)
1019 {
1020 fIntegralPAIySection[0] += b*log(x0/e0);
1021 }
1022 else
1023 {
1024 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1025 }
1026 x0 = fSplineEnergy[i - 1];
1027 x1 = fSplineEnergy[i - 2];
1028 y0 = fDifPAIySection[i - 1];
1029 yy1 = fDifPAIySection[i - 2];
1030
1031 //c = x1/x0;
1032 d = e0/x0;
1033 a = log10(yy1/y0)/log10(x1/x0);
1034 // b0 = log10(y0) - a*log10(x0);
1035 b = y0/pow(x0,a);
1036 a += 1;
1037 if(a == 0)
1038 {
1039 result += b*log(e0/x0);
1040 }
1041 else
1042 {
1043 result += y0*(e0*pow(d,a-1) - x0)/a;
1044 }
1045 a++;
1046 if(a == 0)
1047 {
1048 fIntegralPAIySection[0] += b*log(e0/x0);
1049 }
1050 else
1051 {
1052 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1053 }
1054 return result;
1055
1056}
1057
1058///////////////////////////////////////////////////////////////////////
1059
1060G4double G4PAIySection::SumOverBorderdEdx( G4int i ,
1061 G4double en0 )
1062{
1063 G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
1064
1065 e0 = en0;
1066 x0 = fSplineEnergy[i];
1067 x1 = fSplineEnergy[i+1];
1068 y0 = fDifPAIySection[i];
1069 yy1 = fDifPAIySection[i+1];
1070
1071 //c = x1/x0;
1072 d = e0/x0;
1073 a = log10(yy1/y0)/log10(x1/x0);
1074
1075 G4double b = 0.0;
1076 if(a < 20.) b = y0/pow(x0,a);
1077
1078 a += 2;
1079 if(a == 0)
1080 {
1081 result = b*log(x0/e0);
1082 }
1083 else
1084 {
1085 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1086 }
1087 x0 = fSplineEnergy[i - 1];
1088 x1 = fSplineEnergy[i - 2];
1089 y0 = fDifPAIySection[i - 1];
1090 yy1 = fDifPAIySection[i - 2];
1091
1092 //c = x1/x0;
1093 d = e0/x0;
1094 a = log10(yy1/y0)/log10(x1/x0);
1095
1096 if(a < 20.) b = y0/pow(x0,a);
1097
1098 a += 2;
1099 if(a == 0)
1100 {
1101 result += b*log(e0/x0);
1102 }
1103 else
1104 {
1105 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1106 }
1107 return result;
1108
1109}
1110
1111///////////////////////////////////////////////////////////////////////////////
1112//
1113// Integration of Cerenkov cross-section for the case of
1114// passing across border between intervals
1115
1116G4double G4PAIySection::SumOverBordCerenkov( G4int i ,
1117 G4double en0 )
1118{
1119 G4double x0,x1,y0,yy1,a,e0,c,d,result;
1120
1121 e0 = en0;
1122 x0 = fSplineEnergy[i];
1123 x1 = fSplineEnergy[i+1];
1124 y0 = fdNdxCerenkov[i];
1125 yy1 = fdNdxCerenkov[i+1];
1126
1127 // G4cout<<G4endl;
1128 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1129 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1130 c = x1/x0;
1131 d = e0/x0;
1132 a = log10(yy1/y0)/log10(c);
1133
1134 G4double b = 0.0;
1135 if(a < 20.) b = y0/pow(x0,a);
1136
1137 a += 1.0;
1138 if( a == 0 ) result = b*log(x0/e0);
1139 else result = y0*(x0 - e0*pow(d,a-1))/a;
1140 a += 1.0;
1141
1142 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1143 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1144
1145 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1146
1147 x0 = fSplineEnergy[i - 1];
1148 x1 = fSplineEnergy[i - 2];
1149 y0 = fdNdxCerenkov[i - 1];
1150 yy1 = fdNdxCerenkov[i - 2];
1151
1152 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1153 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1154
1155 c = x1/x0;
1156 d = e0/x0;
1157 a = log10(yy1/y0)/log10(x1/x0);
1158
1159 // G4cout << "a= " << a << G4endl;
1160 if(a > 20.0) b = 0.0;
1161 else b = y0/pow(x0,a);
1162
1163 //G4cout << "b= " << b << G4endl;
1164
1165 a += 1.0;
1166 if( a == 0 ) result += b*log(e0/x0);
1167 else result += y0*(e0*pow(d,a-1) - x0 )/a;
1168 a += 1.0;
1169 //G4cout << "result= " << result << G4endl;
1170
1171 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
1172 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1173
1174 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1175
1176 return result;
1177
1178}
1179
1180///////////////////////////////////////////////////////////////////////////////
1181//
1182// Integration of Plasmon cross-section for the case of
1183// passing across border between intervals
1184
1185G4double G4PAIySection::SumOverBordPlasmon( G4int i ,
1186 G4double en0 )
1187{
1188 G4double x0,x1,y0,yy1,a,c,d,e0,result;
1189
1190 e0 = en0;
1191 x0 = fSplineEnergy[i];
1192 x1 = fSplineEnergy[i+1];
1193 y0 = fdNdxPlasmon[i];
1194 yy1 = fdNdxPlasmon[i+1];
1195
1196 c = x1/x0;
1197 d = e0/x0;
1198 a = log10(yy1/y0)/log10(c);
1199
1200 G4double b = 0.0;
1201 if(a < 20.) b = y0/pow(x0,a);
1202
1203 a += 1.0;
1204 if( a == 0 ) result = b*log(x0/e0);
1205 else result = y0*(x0 - e0*pow(d,a-1))/a;
1206 a += 1.0;
1207
1208 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1209 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1210
1211 x0 = fSplineEnergy[i - 1];
1212 x1 = fSplineEnergy[i - 2];
1213 y0 = fdNdxPlasmon[i - 1];
1214 yy1 = fdNdxPlasmon[i - 2];
1215
1216 c = x1/x0;
1217 d = e0/x0;
1218 a = log10(yy1/y0)/log10(c);
1219
1220 if(a < 20.) b = y0/pow(x0,a);
1221
1222 a += 1.0;
1223 if( a == 0 ) result += b*log(e0/x0);
1224 else result += y0*(e0*pow(d,a-1) - x0)/a;
1225 a += 1.0;
1226
1227 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
1228 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1229
1230 return result;
1231
1232}
1233
1234/////////////////////////////////////////////////////////////////////////
1235//
1236//
1237
1238G4double G4PAIySection::GetStepEnergyLoss( G4double step )
1239{
1240 G4int iTransfer ;
1241 G4long numOfCollisions;
1242 G4double loss = 0.0;
1243 G4double meanNumber, position;
1244
1245 // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl;
1246
1247
1248
1249 meanNumber = fIntegralPAIySection[1]*step;
1250 numOfCollisions = G4Poisson(meanNumber);
1251
1252 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1253
1254 while(numOfCollisions)
1255 {
1256 position = fIntegralPAIySection[1]*G4UniformRand();
1257
1258 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1259 {
1260 if( position >= fIntegralPAIySection[iTransfer] ) break;
1261 }
1262 loss += fSplineEnergy[iTransfer] ;
1263 numOfCollisions--;
1264 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1265 }
1266 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
1267
1268 return loss;
1269}
1270
1271/////////////////////////////////////////////////////////////////////////
1272//
1273//
1274
1275G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
1276{
1277 G4int iTransfer ;
1278 G4long numOfCollisions;
1279 G4double loss = 0.0;
1280 G4double meanNumber, position;
1281
1282 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1283
1284
1285
1286 meanNumber = fIntegralCerenkov[1]*step;
1287 numOfCollisions = G4Poisson(meanNumber);
1288
1289 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1290
1291 while(numOfCollisions)
1292 {
1293 position = fIntegralCerenkov[1]*G4UniformRand();
1294
1295 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1296 {
1297 if( position >= fIntegralCerenkov[iTransfer] ) break;
1298 }
1299 loss += fSplineEnergy[iTransfer] ;
1300 numOfCollisions--;
1301 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1302 }
1303 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1304
1305 return loss;
1306}
1307
1308/////////////////////////////////////////////////////////////////////////
1309//
1310//
1311
1312G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
1313{
1314 G4int iTransfer ;
1315 G4long numOfCollisions;
1316 G4double loss = 0.0;
1317 G4double meanNumber, position;
1318
1319 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1320
1321
1322
1323 meanNumber = fIntegralPlasmon[1]*step;
1324 numOfCollisions = G4Poisson(meanNumber);
1325
1326 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1327
1328 while(numOfCollisions)
1329 {
1330 position = fIntegralPlasmon[1]*G4UniformRand();
1331
1332 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1333 {
1334 if( position >= fIntegralPlasmon[iTransfer] ) break;
1335 }
1336 loss += fSplineEnergy[iTransfer] ;
1337 numOfCollisions--;
1338 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1339 }
1340 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
1341
1342 return loss;
1343}
1344
1345/////////////////////////////////////////////////////////////////////////////
1346//
1347
1348void G4PAIySection::CallError(G4int i, const G4String& methodName) const
1349{
1350 G4String head = "G4PAIySection::" + methodName + "()";
1351 G4ExceptionDescription ed;
1352 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
1353 G4Exception(head,"pai001",FatalException,ed);
1354}
1355
1356/////////////////////////////////////////////////////////////////////////////
1357//
1358// Init array of Lorentz factors
1359//
1360
1361G4int G4PAIySection::fNumberOfGammas = 111;
1362
1363const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1
1364{
13650.0,
13661.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
13671.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
13681.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
13691.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
13702.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
13713.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
13725.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
13738.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
13741.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
13752.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
13765.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
13771.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
13781.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
13793.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
13806.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
13811.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
13822.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
13834.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
13848.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
13851.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
13863.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
13875.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
13881.065799e+05
1389};
1390
1391///////////////////////////////////////////////////////////////////////
1392//
1393// The number of gamma for creation of spline (near ion-min , G ~ 4 )
1394//
1395
1396const G4int G4PAIySection::fRefGammaNumber = 29;
1397
1398//
1399// end of G4PAIySection implementation file
1400//
1401////////////////////////////////////////////////////////////////////////////
1402
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4double
double G4double
Definition: G4Types.hh:83
G4long
long G4long
Definition: G4Types.hh:87
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:175
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:197
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:181
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:212
G4PAIySection::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition: G4PAIySection.cc:1275
G4PAIySection::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition: G4PAIySection.cc:534
G4PAIySection::IntegralPlasmon
void IntegralPlasmon()
Definition: G4PAIySection.cc:816
G4PAIySection::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4PAIySection.cc:511
G4PAIySection::GetStepPlasmonLoss
G4double GetStepPlasmonLoss(G4double step)
Definition: G4PAIySection.cc:1312
G4PAIySection::SumOverBordPlasmon
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
Definition: G4PAIySection.cc:1185
G4PAIySection::NormShift
void NormShift(G4double betaGammaSq)
Definition: G4PAIySection.cc:329
G4PAIySection::ComputeLowEnergyCof
void ComputeLowEnergyCof(const G4Material *material)
Definition: G4PAIySection.cc:245
G4PAIySection::SumOverBorderdEdx
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
Definition: G4PAIySection.cc:1060
G4PAIySection::SumOverInterPlasmon
G4double SumOverInterPlasmon(G4int intervalNumber)
Definition: G4PAIySection.cc:956
G4PAIySection::IntegralCerenkov
void IntegralCerenkov()
Definition: G4PAIySection.cc:785
G4PAIySection::IntegralPAIySection
void IntegralPAIySection()
Definition: G4PAIySection.cc:754
G4PAIySection::GetStepEnergyLoss
G4double GetStepEnergyLoss(G4double step)
Definition: G4PAIySection.cc:1238
G4PAIySection::SumOverInterCerenkov
G4double SumOverInterCerenkov(G4int intervalNumber)
Definition: G4PAIySection.cc:919
G4PAIySection::SumOverInterval
G4double SumOverInterval(G4int intervalNumber)
Definition: G4PAIySection.cc:843
G4PAIySection::SumOverBorder
G4double SumOverBorder(G4int intervalNumber, G4double energy)
Definition: G4PAIySection.cc:990
G4PAIySection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIySection.cc:118
G4PAIySection::SumOverBordCerenkov
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
Definition: G4PAIySection.cc:1116
G4PAIySection::RutherfordIntegral
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
Definition: G4PAIySection.cc:489
G4PAIySection::DifPAIySection
G4double DifPAIySection(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIySection.cc:589
G4PAIySection::SumOverIntervaldEdx
G4double SumOverIntervaldEdx(G4int intervalNumber)
Definition: G4PAIySection.cc:885
G4PAIySection::PAIdNdxCerenkov
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIySection.cc:656
G4PAIySection::GetLorentzFactor
G4double GetLorentzFactor(G4int j) const
Definition: G4PAIySection.cc:109
G4PAIySection::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
Definition: G4PAIySection.cc:715
G4PAIySection::SplainPAI
void SplainPAI(G4double betaGammaSq)
Definition: G4PAIySection.cc:402
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.cc:723
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition: G4SandiaTable.cc:1066