Geant4 11.2.2
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, G4double energy1 )
512{
513 G4double energy2,energy3,energy4,result;
514
515 energy2 = energy1*energy1;
516 energy3 = energy2*energy1;
517 energy4 = energy3*energy1;
518
519 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
520 result *=hbarc/energy1;
521
522 return result;
523
524} // end of ImPartDielectricConst
525
526
527//////////////////////////////////////////////////////////////////////////////
528//
529// Real part of dielectric constant minus unit: epsilon_1 - 1
530// (G4double enb - energy point)
531//
532
533G4double G4PAIySection::RePartDielectricConst(G4double enb)
534{
535 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
536 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
537
538 x0 = enb;
539 result = 0;
540
541 for(G4int i=1;i<=fIntervalNumber-1;++i)
542 {
543 x1 = fEnergyInterval[i];
544 x2 = fEnergyInterval[i+1];
545 xx1 = x1 - x0;
546 xx2 = x2 - x0;
547 xx12 = xx2/xx1;
548
549 if(xx12<0.)
550 {
551 xx12 = -xx12;
552 }
553 xln1 = log(x2/x1);
554 xln2 = log(xx12);
555 xln3 = log((x2 + x0)/(x1 + x0));
556 x02 = x0*x0;
557 x03 = x02*x0;
558 x04 = x03*x0;
559 x05 = x04*x0;
560 G4double x12 = x1*x2;
561 c1 = (x2 - x1)/x12;
562 c2 = (x2 - x1)*(x2 +x1)/(x12*x12);
563 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/(x12*x12*x12);
564
565 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
566 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
567 result -= fA3[i]*c2/2/x02;
568 result -= fA4[i]*c3/3/x02;
569
570 cof1 = fA1[i]/x02 + fA3[i]/x04;
571 cof2 = fA2[i]/x03 + fA4[i]/x05;
572
573 result += 0.5*(cof1 +cof2)*xln2;
574 result += 0.5*(cof1 - cof2)*xln3;
575 }
576 result *= 2*hbarc/pi;
577
578 return result;
579
580} // end of RePartDielectricConst
581
582//////////////////////////////////////////////////////////////////////
583//
584// PAI differential cross-section in terms of
585// simplified Allison's equation
586//
587
588G4double G4PAIySection::DifPAIySection( G4int i ,
589 G4double betaGammaSq )
590{
591 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
592 be2 = betaGammaSq/(1 + betaGammaSq);
593 beta = std::sqrt(be2);
594 cof = 1;
595 x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
596
597 if( betaGammaSq < 0.01 ) x2 = log(be2);
598 else
599 {
600 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
601 (1/betaGammaSq - fRePartDielectricConst[i]) +
602 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
603 }
604 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
605 {
606 x6=0;
607 }
608 else
609 {
610 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
611 x5 = -1 - fRePartDielectricConst[i] +
612 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
613 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
614
615 x7 = std::atan2(fImPartDielectricConst[i],x3);
616 x6 = x5 * x7;
617 }
618 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
619 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
620 fImPartDielectricConst[i]*fImPartDielectricConst[i];
621
622 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
623 result = std::max(result, 1.0e-8);
624 result *= fine_structure_const/(be2*pi);
625 // low energy correction
626
627 G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
628
629 result *= (1 - std::exp(-beta/(betaBohr*lowCof)));
630 if(x8 > 0.)
631 {
632 result /= x8;
633 }
634 return result;
635
636} // end of DifPAIySection
637
638//////////////////////////////////////////////////////////////////////////
639//
640// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
641
642G4double G4PAIySection::PAIdNdxCerenkov( G4int i, G4double betaGammaSq )
643{
644 G4double logarithm, x3, x5, argument, modul2, dNdxC;
645 G4double be2, be4;
646
647 be2 = betaGammaSq/(1 + betaGammaSq);
648 be4 = be2*be2;
649
650 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
651 else
652 {
653 logarithm = -std::log( (1/betaGammaSq - fRePartDielectricConst[i])*
654 (1/betaGammaSq - fRePartDielectricConst[i]) +
655 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
656 logarithm += std::log(1+1.0/betaGammaSq);
657 }
658
659 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
660 {
661 argument = 0.0;
662 }
663 else
664 {
665 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
666 x5 = -1.0 - fRePartDielectricConst[i] +
667 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
668 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
669 if( x3 == 0.0 ) argument = 0.5*pi;
670 else argument = std::atan2(fImPartDielectricConst[i],x3);
671 argument *= x5 ;
672 }
673 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
674
675 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
676
677 dNdxC *= fine_structure_const/be2/pi;
678
679 dNdxC *= (1 - std::exp(-be4/betaBohr4));
680
681 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
682 fImPartDielectricConst[i]*fImPartDielectricConst[i];
683 if(modul2 > 0.)
684 {
685 dNdxC /= modul2;
686 }
687 return dNdxC;
688
689} // end of PAIdNdxCerenkov
690
691//////////////////////////////////////////////////////////////////////////
692//
693// Calculation od dN/dx of collisions with creation of longitudinal EM
694// excitations (plasmons, delta-electrons)
695
696G4double G4PAIySection::PAIdNdxPlasmon( G4int i, G4double betaGammaSq )
697{
698 G4double cof, resonance, modul2, dNdxP;
699 G4double be2, be4;
700
701 cof = 1;
702
703 be2 = betaGammaSq/(1 + betaGammaSq);
704 be4 = be2*be2;
705
706 resonance = std::log(2*electron_mass_c2*be2/fSplineEnergy[i]);
707 resonance *= fImPartDielectricConst[i]/hbarc;
708
709 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
710
711 dNdxP = std::max(dNdxP, 1.0e-8);
712
713 dNdxP *= fine_structure_const/be2/pi;
714 dNdxP *= (1 - std::exp(-be4/betaBohr4));
715
716 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
717 fImPartDielectricConst[i]*fImPartDielectricConst[i];
718 if(modul2 > 0.)
719 {
720 dNdxP /= modul2;
721 }
722 return dNdxP;
723
724} // end of PAIdNdxPlasmon
725
726////////////////////////////////////////////////////////////////////////
727//
728// Calculation of the PAI integral cross-section
729// fIntegralPAIySection[1] = specific primary ionisation, 1/cm
730// and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm
731
732void G4PAIySection::IntegralPAIySection()
733{
734 fIntegralPAIySection[fSplineNumber] = 0;
735 fIntegralPAIdEdx[fSplineNumber] = 0;
736 fIntegralPAIySection[0] = 0;
737 G4int k = fIntervalNumber -1;
738
739 for(G4int i = fSplineNumber-1; i >= 1; i--)
740 {
741 if(fSplineEnergy[i] >= fEnergyInterval[k])
742 {
743 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
744 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
745 }
746 else
747 {
748 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
749 SumOverBorder(i+1,fEnergyInterval[k]);
750 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
751 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
752 k--;
753 }
754 }
755} // end of IntegralPAIySection
756
757////////////////////////////////////////////////////////////////////////
758//
759// Calculation of the PAI Cerenkov integral cross-section
760// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
761// and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
762
763void G4PAIySection::IntegralCerenkov()
764{
765 G4int i, k;
766 fIntegralCerenkov[fSplineNumber] = 0;
767 fIntegralCerenkov[0] = 0;
768 k = fIntervalNumber -1;
769
770 for( i = fSplineNumber-1; i >= 1; i-- )
771 {
772 if(fSplineEnergy[i] >= fEnergyInterval[k])
773 {
774 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
775 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
776 }
777 else
778 {
779 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
780 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
781 k--;
782 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
783 }
784 }
785
786} // end of IntegralCerenkov
787
788////////////////////////////////////////////////////////////////////////
789//
790// Calculation of the PAI Plasmon integral cross-section
791// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
792// and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
793
794void G4PAIySection::IntegralPlasmon()
795{
796 fIntegralPlasmon[fSplineNumber] = 0;
797 fIntegralPlasmon[0] = 0;
798 G4int k = fIntervalNumber -1;
799 for(G4int i=fSplineNumber-1;i>=1;i--)
800 {
801 if(fSplineEnergy[i] >= fEnergyInterval[k])
802 {
803 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
804 }
805 else
806 {
807 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
808 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
809 k--;
810 }
811 }
812} // end of IntegralPlasmon
813
814//////////////////////////////////////////////////////////////////////
815//
816// Calculation the PAI integral cross-section inside
817// of interval of continuous values of photo-ionisation
818// cross-section. Parameter 'i' is the number of interval.
819
820G4double G4PAIySection::SumOverInterval( G4int i )
821{
822 G4double x0,x1,y0,yy1,a,b,c,result;
823
824 x0 = fSplineEnergy[i];
825 x1 = fSplineEnergy[i+1];
826
827 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
828
829 y0 = fDifPAIySection[i];
830 yy1 = fDifPAIySection[i+1];
831 //G4cout << "## x0= " << x0 << " x1= " << x1 << G4endl;
832 c = x1/x0;
833 //G4cout << "c= " << c << " y0= " << y0 << " yy1= " << yy1 << G4endl;
834 a = log10(yy1/y0)/log10(c);
835 //G4cout << "a= " << a << G4endl;
836
837 b = 0.0;
838 if(a < 20.) b = y0/pow(x0,a);
839
840 a += 1;
841 if(a == 0)
842 {
843 result = b*log(x1/x0);
844 }
845 else
846 {
847 result = y0*(x1*pow(c,a-1) - x0)/a;
848 }
849 a++;
850 if(a == 0)
851 {
852 fIntegralPAIySection[0] += b*log(x1/x0);
853 }
854 else
855 {
856 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
857 }
858 return result;
859
860} // end of SumOverInterval
861
862/////////////////////////////////
863
864G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
865{
866 G4double x0,x1,y0,yy1,a,b,c,result;
867
868 x0 = fSplineEnergy[i];
869 x1 = fSplineEnergy[i+1];
870
871 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
872
873 y0 = fDifPAIySection[i];
874 yy1 = fDifPAIySection[i+1];
875 c = x1/x0;
876 a = log10(yy1/y0)/log10(c);
877
878 b = 0.0;
879 if(a < 20.) b = y0/pow(x0,a);
880
881 a += 2;
882 if(a == 0)
883 {
884 result = b*log(x1/x0);
885 }
886 else
887 {
888 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
889 }
890 return result;
891
892} // end of SumOverInterval
893
894//////////////////////////////////////////////////////////////////////
895//
896// Calculation the PAI Cerenkov integral cross-section inside
897// of interval of continuous values of photo-ionisation Cerenkov
898// cross-section. Parameter 'i' is the number of interval.
899
900G4double G4PAIySection::SumOverInterCerenkov( G4int i )
901{
902 G4double x0,x1,y0,yy1,a,c,result;
903
904 x0 = fSplineEnergy[i];
905 x1 = fSplineEnergy[i+1];
906
907 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
908
909 y0 = fdNdxCerenkov[i];
910 yy1 = fdNdxCerenkov[i+1];
911 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
912 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
913
914 c = x1/x0;
915 a = log10(yy1/y0)/log10(c);
916 G4double b = 0.0;
917 if(a < 20.) b = y0/pow(x0,a);
918
919 a += 1.0;
920 if(a == 0) result = b*log(c);
921 else result = y0*(x1*pow(c,a-1) - x0)/a;
922 a += 1.0;
923
924 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
925 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
926 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
927 return result;
928
929} // end of SumOverInterCerenkov
930
931//////////////////////////////////////////////////////////////////////
932//
933// Calculation the PAI Plasmon integral cross-section inside
934// of interval of continuous values of photo-ionisation Plasmon
935// cross-section. Parameter 'i' is the number of interval.
936
937G4double G4PAIySection::SumOverInterPlasmon( G4int i )
938{
939 G4double x0,x1,y0,yy1,a,c,result;
940
941 x0 = fSplineEnergy[i];
942 x1 = fSplineEnergy[i+1];
943
944 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
945
946 y0 = fdNdxPlasmon[i];
947 yy1 = fdNdxPlasmon[i+1];
948 c = x1/x0;
949 a = log10(yy1/y0)/log10(c);
950
951 G4double b = 0.0;
952 if(a < 20.) b = y0/pow(x0,a);
953
954 a += 1.0;
955 if(a == 0) result = b*log(x1/x0);
956 else result = y0*(x1*pow(c,a-1) - x0)/a;
957 a += 1.0;
958
959 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
960 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
961
962 return result;
963
964} // end of SumOverInterPlasmon
965
966///////////////////////////////////////////////////////////////////////////////
967//
968// Integration of PAI cross-section for the case of
969// passing across border between intervals
970
971G4double G4PAIySection::SumOverBorder( G4int i ,
972 G4double en0 )
973{
974 G4double x0,x1,y0,yy1,a,d,e0,result;
975
976 e0 = en0;
977 x0 = fSplineEnergy[i];
978 x1 = fSplineEnergy[i+1];
979 y0 = fDifPAIySection[i];
980 yy1 = fDifPAIySection[i+1];
981
982 d = e0/x0;
983 a = log10(yy1/y0)/log10(x1/x0);
984
985 G4double b = 0.0;
986 if(a < 20.) b = y0/pow(x0,a);
987
988 a += 1;
989 if(a == 0)
990 {
991 result = b*log(x0/e0);
992 }
993 else
994 {
995 result = y0*(x0 - e0*pow(d,a-1))/a;
996 }
997 a++;
998 if(a == 0)
999 {
1000 fIntegralPAIySection[0] += b*log(x0/e0);
1001 }
1002 else
1003 {
1004 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1005 }
1006 x0 = fSplineEnergy[i - 1];
1007 x1 = fSplineEnergy[i - 2];
1008 y0 = fDifPAIySection[i - 1];
1009 yy1 = fDifPAIySection[i - 2];
1010
1011 //c = x1/x0;
1012 d = e0/x0;
1013 a = log10(yy1/y0)/log10(x1/x0);
1014
1015 b = 0.0;
1016 if(a < 20.) b = y0/pow(x0,a);
1017
1018 a += 1;
1019 if(a == 0)
1020 {
1021 result += b*log(e0/x0);
1022 }
1023 else
1024 {
1025 result += y0*(e0*pow(d,a-1) - x0)/a;
1026 }
1027 a++;
1028 if(a == 0)
1029 {
1030 fIntegralPAIySection[0] += b*log(e0/x0);
1031 }
1032 else
1033 {
1034 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1035 }
1036 return result;
1037
1038}
1039
1040///////////////////////////////////////////////////////////////////////
1041
1042G4double G4PAIySection::SumOverBorderdEdx( G4int i ,
1043 G4double en0 )
1044{
1045 G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result;
1046
1047 e0 = en0;
1048 x0 = fSplineEnergy[i];
1049 x1 = fSplineEnergy[i+1];
1050 y0 = fDifPAIySection[i];
1051 yy1 = fDifPAIySection[i+1];
1052
1053 d = e0/x0;
1054 a = log10(yy1/y0)/log10(x1/x0);
1055
1056 G4double b = 0.0;
1057 if(a < 20.) b = y0/pow(x0,a);
1058
1059 a += 2;
1060 if(a == 0)
1061 {
1062 result = b*log(x0/e0);
1063 }
1064 else
1065 {
1066 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1067 }
1068 x0 = fSplineEnergy[i - 1];
1069 x1 = fSplineEnergy[i - 2];
1070 y0 = fDifPAIySection[i - 1];
1071 yy1 = fDifPAIySection[i - 2];
1072
1073 d = e0/x0;
1074 a = log10(yy1/y0)/log10(x1/x0);
1075
1076 b = 0.0;
1077 if(a < 20.) b = y0/pow(x0,a);
1078
1079 a += 2;
1080 if(a == 0)
1081 {
1082 result += b*log(e0/x0);
1083 }
1084 else
1085 {
1086 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1087 }
1088 return result;
1089}
1090
1091///////////////////////////////////////////////////////////////////////////////
1092//
1093// Integration of Cerenkov cross-section for the case of
1094// passing across border between intervals
1095
1096G4double G4PAIySection::SumOverBordCerenkov( G4int i ,
1097 G4double en0 )
1098{
1099 G4double x0,x1,y0,yy1,a,e0,c,d,result;
1100
1101 e0 = en0;
1102 x0 = fSplineEnergy[i];
1103 x1 = fSplineEnergy[i+1];
1104 y0 = fdNdxCerenkov[i];
1105 yy1 = fdNdxCerenkov[i+1];
1106
1107 // G4cout<<G4endl;
1108 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1109 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1110 c = x1/x0;
1111 d = e0/x0;
1112 a = log10(yy1/y0)/log10(c);
1113
1114 G4double b = 0.0;
1115 if(a < 20.) b = y0/pow(x0,a);
1116
1117 a += 1.0;
1118 if( a == 0 ) result = b*log(x0/e0);
1119 else result = y0*(x0 - e0*pow(d,a-1))/a;
1120 a += 1.0;
1121
1122 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1123 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1124
1125 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1126
1127 x0 = fSplineEnergy[i - 1];
1128 x1 = fSplineEnergy[i - 2];
1129 y0 = fdNdxCerenkov[i - 1];
1130 yy1 = fdNdxCerenkov[i - 2];
1131
1132 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1133 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1134
1135 c = x1/x0;
1136 d = e0/x0;
1137 a = log10(yy1/y0)/log10(x1/x0);
1138
1139 // G4cout << "a= " << a << G4endl;
1140 if(a > 20.0) b = 0.0;
1141 else b = y0/pow(x0,a);
1142
1143 //G4cout << "b= " << b << G4endl;
1144
1145 a += 1.0;
1146 if( a == 0 ) result += b*log(e0/x0);
1147 else result += y0*(e0*pow(d,a-1) - x0 )/a;
1148 a += 1.0;
1149 //G4cout << "result= " << result << G4endl;
1150
1151 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
1152 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1153
1154 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1155
1156 return result;
1157}
1158
1159///////////////////////////////////////////////////////////////////////////////
1160//
1161// Integration of Plasmon cross-section for the case of
1162// passing across border between intervals
1163
1164G4double G4PAIySection::SumOverBordPlasmon( G4int i ,
1165 G4double en0 )
1166{
1167 G4double x0,x1,y0,yy1,a,c,d,e0,result;
1168
1169 e0 = en0;
1170 x0 = fSplineEnergy[i];
1171 x1 = fSplineEnergy[i+1];
1172 y0 = fdNdxPlasmon[i];
1173 yy1 = fdNdxPlasmon[i+1];
1174
1175 c = x1/x0;
1176 d = e0/x0;
1177 a = log10(yy1/y0)/log10(c);
1178
1179 G4double b = 0.0;
1180 if(a < 20.) b = y0/pow(x0,a);
1181
1182 a += 1.0;
1183 if( a == 0 ) result = b*log(x0/e0);
1184 else result = y0*(x0 - e0*pow(d,a-1))/a;
1185 a += 1.0;
1186
1187 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1188 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1189
1190 x0 = fSplineEnergy[i - 1];
1191 x1 = fSplineEnergy[i - 2];
1192 y0 = fdNdxPlasmon[i - 1];
1193 yy1 = fdNdxPlasmon[i - 2];
1194
1195 c = x1/x0;
1196 d = e0/x0;
1197 a = log10(yy1/y0)/log10(c);
1198
1199 if(a < 20.) b = y0/pow(x0,a);
1200
1201 a += 1.0;
1202 if( a == 0 ) result += b*log(e0/x0);
1203 else result += y0*(e0*pow(d,a-1) - x0)/a;
1204 a += 1.0;
1205
1206 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
1207 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1208
1209 return result;
1210
1211}
1212
1213/////////////////////////////////////////////////////////////////////////
1214//
1215//
1216
1217G4double G4PAIySection::GetStepEnergyLoss( G4double step )
1218{
1219 G4int iTransfer ;
1220 G4long numOfCollisions;
1221 G4double loss = 0.0;
1222 G4double meanNumber, position;
1223
1224 // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl;
1225
1226
1227
1228 meanNumber = fIntegralPAIySection[1]*step;
1229 numOfCollisions = G4Poisson(meanNumber);
1230
1231 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1232
1233 while(numOfCollisions)
1234 {
1235 position = fIntegralPAIySection[1]*G4UniformRand();
1236
1237 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1238 {
1239 if( position >= fIntegralPAIySection[iTransfer] ) break;
1240 }
1241 loss += fSplineEnergy[iTransfer] ;
1242 numOfCollisions--;
1243 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1244 }
1245 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
1246
1247 return loss;
1248}
1249
1250/////////////////////////////////////////////////////////////////////////
1251//
1252//
1253
1254G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
1255{
1256 G4int iTransfer ;
1257 G4long numOfCollisions;
1258 G4double loss = 0.0;
1259 G4double meanNumber, position;
1260
1261 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1262
1263
1264
1265 meanNumber = fIntegralCerenkov[1]*step;
1266 numOfCollisions = G4Poisson(meanNumber);
1267
1268 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1269
1270 while(numOfCollisions)
1271 {
1272 position = fIntegralCerenkov[1]*G4UniformRand();
1273
1274 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1275 {
1276 if( position >= fIntegralCerenkov[iTransfer] ) break;
1277 }
1278 loss += fSplineEnergy[iTransfer] ;
1279 numOfCollisions--;
1280 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1281 }
1282 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1283
1284 return loss;
1285}
1286
1287/////////////////////////////////////////////////////////////////////////
1288//
1289//
1290
1291G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
1292{
1293 G4int iTransfer ;
1294 G4long numOfCollisions;
1295 G4double loss = 0.0;
1296 G4double meanNumber, position;
1297
1298 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl;
1299
1300
1301
1302 meanNumber = fIntegralPlasmon[1]*step;
1303 numOfCollisions = G4Poisson(meanNumber);
1304
1305 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1306
1307 while(numOfCollisions)
1308 {
1309 position = fIntegralPlasmon[1]*G4UniformRand();
1310
1311 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
1312 {
1313 if( position >= fIntegralPlasmon[iTransfer] ) break;
1314 }
1315 loss += fSplineEnergy[iTransfer] ;
1316 numOfCollisions--;
1317 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1318 }
1319 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
1320
1321 return loss;
1322}
1323
1324/////////////////////////////////////////////////////////////////////////////
1325//
1326
1327void G4PAIySection::CallError(G4int i, const G4String& methodName) const
1328{
1329 G4String head = "G4PAIySection::" + methodName + "()";
1330 G4ExceptionDescription ed;
1331 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber;
1332 G4Exception(head,"pai001",FatalException,ed);
1333}
1334
1335/////////////////////////////////////////////////////////////////////////////
1336//
1337// Init array of Lorentz factors
1338//
1339
1340G4int G4PAIySection::fNumberOfGammas = 111;
1341
1342const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1
1343{
13440.0,
13451.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
13461.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
13471.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
13481.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
13492.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
13503.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
13515.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
13528.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
13531.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
13542.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
13555.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
13561.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
13571.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
13583.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
13596.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
13601.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
13612.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
13624.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
13638.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
13641.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
13653.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
13665.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
13671.065799e+05
1368};
1369
1370///////////////////////////////////////////////////////////////////////
1371//
1372// The number of gamma for creation of spline (near ion-min , G ~ 4 )
1373//
1374
1375const G4int G4PAIySection::fRefGammaNumber = 29;
1376
1377//
1378// end of G4PAIySection implementation file
1379//
1380////////////////////////////////////////////////////////////////////////////
1381
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
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:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4Element::GetZ
G4double GetZ() const
Definition G4Element.hh:119
G4Material::GetDensity
G4double GetDensity() const
Definition G4Material.hh:174
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition G4Material.hh:192
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition G4Material.hh:203
G4Material::GetNumberOfElements
std::size_t GetNumberOfElements() const
Definition G4Material.hh:180
G4PAIySection::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition G4PAIySection.cc:1254
G4PAIySection::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition G4PAIySection.cc:533
G4PAIySection::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition G4PAIySection.cc:511
G4PAIySection::GetStepPlasmonLoss
G4double GetStepPlasmonLoss(G4double step)
Definition G4PAIySection.cc:1291
G4PAIySection::SumOverBordPlasmon
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
Definition G4PAIySection.cc:1164
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:1042
G4PAIySection::SumOverInterPlasmon
G4double SumOverInterPlasmon(G4int intervalNumber)
Definition G4PAIySection.cc:937
G4PAIySection::IntegralCerenkov
void IntegralCerenkov()
Definition G4PAIySection.cc:763
G4PAIySection::IntegralPAIySection
void IntegralPAIySection()
Definition G4PAIySection.cc:732
G4PAIySection::GetStepEnergyLoss
G4double GetStepEnergyLoss(G4double step)
Definition G4PAIySection.cc:1217
G4PAIySection::SumOverInterCerenkov
G4double SumOverInterCerenkov(G4int intervalNumber)
Definition G4PAIySection.cc:900
G4PAIySection::SumOverInterval
G4double SumOverInterval(G4int intervalNumber)
Definition G4PAIySection.cc:820
G4PAIySection::SumOverBorder
G4double SumOverBorder(G4int intervalNumber, G4double energy)
Definition G4PAIySection.cc:971
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:1096
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:588
G4PAIySection::SumOverIntervaldEdx
G4double SumOverIntervaldEdx(G4int intervalNumber)
Definition G4PAIySection.cc:864
G4PAIySection::PAIdNdxCerenkov
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
Definition G4PAIySection.cc:642
G4PAIySection::GetLorentzFactor
G4double GetLorentzFactor(G4int j) const
Definition G4PAIySection.cc:109
G4PAIySection::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
Definition G4PAIySection.cc:696
G4PAIySection::SplainPAI
void SplainPAI(G4double betaGammaSq)
Definition G4PAIySection.cc:402
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition G4SandiaTable.cc:650
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int) const
Definition G4SandiaTable.cc:949