Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4InitXscPAI.cc
이 파일의 문서화 페이지로 가기
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 // $Id: G4InitXscPAI.cc 96934 2016-05-18 09:10:41Z gcosmo $
27 //
28 //
29 // G4InitXscPAI.cc -- class implementation file
30 //
31 // GEANT 4 class implementation file
32 //
33 // For information related to this code, please, contact
34 // the Geant4 Collaboration.
35 //
36 // R&D: Vladimir.Grichine@cern.ch
37 //
38 // History:
39 //
40 
41 
42 
43 #include "G4InitXscPAI.hh"
44 
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ios.hh"
49 #include "G4Poisson.hh"
50 #include "G4Integrator.hh"
51 #include "G4Material.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4SandiaTable.hh"
54 
55 
56 
57 // Local class constants
58 
59 const G4double G4InitXscPAI::fDelta = 0.005 ; // energy shift from interval border
60 const G4int G4InitXscPAI::fPAIbin = 100 ; // size of energy transfer vectors
61 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border
62 
64 //
65 // Constructor
66 //
67 
68 using namespace std;
69 
71  : fPAIxscVector(nullptr),
72  fPAIdEdxVector(nullptr),
73  fPAIphotonVector(nullptr),
74  fPAIelectronVector(nullptr),
75  fChCosSqVector(nullptr),
76  fChWidthVector(nullptr)
77 {
78  G4int i, j, matIndex;
79 
80  fDensity = matCC->GetMaterial()->GetDensity();
82  matIndex = matCC->GetMaterial()->GetIndex();
83 
84  fSandia = new G4SandiaTable(matIndex);
86 
88 
89  for (i = 0; i < fIntervalNumber; i++)
90  {
91  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
92  }
93  for (i = 0; i < fIntervalNumber; i++)
94  {
95  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
96 
97  for(j = 1; j < 5 ; j++)
98  {
99  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
100  }
101  }
103  Normalisation();
104  fBetaGammaSq = fTmax = 0.0;
106 }
107 
108 
109 
110 
112 //
113 // Destructor
114 
116 {
117  if(fPAIxscVector) delete fPAIxscVector;
118  if(fPAIdEdxVector) delete fPAIdEdxVector;
121  if(fChCosSqVector) delete fChCosSqVector;
122  if(fChWidthVector) delete fChWidthVector;
123  delete fSandia;
124  delete fMatSandiaMatrix;
125 }
126 
128 //
129 // Kill close intervals, recalculate fIntervalNumber
130 
132 {
133  G4int i, j, k;
134  G4double energy1, energy2;
135 
136  for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
137  {
138  energy1 = (*(*fMatSandiaMatrix)[i])[0];
139  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
140 
141  if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ;
142  else
143  {
144  for(j = i; j < fIntervalNumber-1; j++)
145  {
146  for( k = 0; k < 5; k++ )
147  {
148  (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
149  }
150  }
151  fIntervalNumber-- ;
152  i-- ;
153  }
154  }
155 
156 }
157 
159 //
160 // Kill close intervals, recalculate fIntervalNumber
161 
163 {
164  G4int i, j;
165  G4double energy1, energy2, /*delta,*/ cof; // , shift;
166 
167  energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
168  energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
169 
170 
171  cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
172 
173  for( i = fIntervalNumber-2; i >= 0; i-- )
174  {
175  energy1 = (*(*fMatSandiaMatrix)[i])[0];
176  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
177 
178  cof += RutherfordIntegral(i,energy1,energy2);
179  // G4cout<<"norm. cof = "<<cof<<G4endl;
180  }
183  //delta = fNormalizationCof - cof;
184  fNormalizationCof /= cof;
185  // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
186  // <<"; at delta ="<<delta<<G4endl ;
187 
188  for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
189  {
190  for(j = 1; j < 5 ; j++)
191  {
192  (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
193  }
194  }
195  /*
196  if(delta > 0) // shift the first energy interval
197  {
198  for(i=1;i<100;i++)
199  {
200  energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
201  energy2 = (*(*fMatSandiaMatrix)[0])[0];
202  shift = RutherfordIntegral(0,energy1,energy2);
203  G4cout<<shift<<"\t";
204  if(shift >= delta) break;
205  }
206  (*(*fMatSandiaMatrix)[0])[0] = energy1;
207  cof += shift;
208  }
209  else if(delta < 0)
210  {
211  for(i=1;i<100;i++)
212  {
213  energy1 = (*(*fMatSandiaMatrix)[0])[0];
214  energy2 = (*(*fMatSandiaMatrix)[0])[0] +
215  ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
216  shift = RutherfordIntegral(0,energy1,energy2);
217  if( shift >= std::abs(delta) ) break;
218  }
219  (*(*fMatSandiaMatrix)[0])[0] = energy2;
220  cof -= shift;
221  }
222  G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
223  <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ;
224  */
225 }
226 
228 //
229 // Integration over electrons that could be considered
230 // quasi-free at energy transfer of interest
231 
233  G4double x1,
234  G4double x2 )
235 {
236  G4double c1, c2, c3, a1, a2, a3, a4 ;
237 
238  a1 = (*(*fMatSandiaMatrix)[k])[1];
239  a2 = (*(*fMatSandiaMatrix)[k])[2];
240  a3 = (*(*fMatSandiaMatrix)[k])[3];
241  a4 = (*(*fMatSandiaMatrix)[k])[4];
242  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
243  c1 = (x2 - x1)/x1/x2 ;
244  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
245  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
246  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
247 
248  return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
249 
250 } // end of RutherfordIntegral
251 
253 //
254 // Integrate photo-absorption cross-section from I1 up to omega
255 
257 {
258  G4int i;
259  G4double energy1, energy2, result = 0.;
260 
261  for( i = 0; i <= fIntervalTmax; i++ )
262  {
263  if(i == fIntervalTmax)
264  {
265  energy1 = (*(*fMatSandiaMatrix)[i])[0];
266  result += RutherfordIntegral(i,energy1,omega);
267  }
268  else
269  {
270  if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
271  {
272  energy1 = (*(*fMatSandiaMatrix)[i])[0];
273  result += RutherfordIntegral(i,energy1,omega);
274  break;
275  }
276  else
277  {
278  energy1 = (*(*fMatSandiaMatrix)[i])[0];
279  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
280  result += RutherfordIntegral(i,energy1,energy2);
281  }
282  }
283  // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
284  }
285  return result;
286 }
287 
288 
290 //
291 // Imaginary part of dielectric constant
292 // (G4int k - interval number, G4double en1 - energy point)
293 
295  G4double energy1 )
296 {
297  G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
298 
299  a1 = (*(*fMatSandiaMatrix)[k])[1];
300  a2 = (*(*fMatSandiaMatrix)[k])[2];
301  a3 = (*(*fMatSandiaMatrix)[k])[3];
302  a4 = (*(*fMatSandiaMatrix)[k])[4];
303 
304  energy2 = energy1*energy1;
305  energy3 = energy2*energy1;
306  energy4 = energy3*energy1;
307 
308  result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
309  result *= hbarc/energy1 ;
310 
311  return result ;
312 
313 } // end of ImPartDielectricConst
314 
316 //
317 // Modulus squared of dielectric constant
318 // (G4int k - interval number, G4double omega - energy point)
319 
321  G4double omega )
322 {
323  G4double eIm2, eRe2, result;
324 
325  result = ImPartDielectricConst(k,omega);
326  eIm2 = result*result;
327 
328  result = RePartDielectricConst(omega);
329  eRe2 = result*result;
330 
331  result = eIm2 + eRe2;
332 
333  return result ;
334 }
335 
336 
338 //
339 // Real part of dielectric constant minus unit: epsilon_1 - 1
340 // (G4double enb - energy point)
341 //
342 
344 {
345  G4int i;
346  G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
347  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
348 
349  x0 = enb ;
350  result = 0 ;
351 
352  for( i = 0; i < fIntervalNumber-1; i++)
353  {
354  x1 = (*(*fMatSandiaMatrix)[i])[0];
355  x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
356 
357  a1 = (*(*fMatSandiaMatrix)[i])[1];
358  a2 = (*(*fMatSandiaMatrix)[i])[2];
359  a3 = (*(*fMatSandiaMatrix)[i])[3];
360  a4 = (*(*fMatSandiaMatrix)[i])[4];
361 
362  if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta )
363  {
364  if(x0 >= x1) x0 = x1*(1+fDelta);
365  else x0 = x1*(1-fDelta);
366  }
367  if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta )
368  {
369  if(x0 >= x2) x0 = x2*(1+fDelta);
370  else x0 = x2*(1-fDelta);
371  }
372  xx1 = x1 - x0 ;
373  xx2 = x2 - x0 ;
374  xx12 = xx2/xx1 ;
375 
376  if( xx12 < 0 ) xx12 = -xx12;
377 
378  xln1 = log(x2/x1) ;
379  xln2 = log(xx12) ;
380  xln3 = log((x2 + x0)/(x1 + x0)) ;
381 
382  x02 = x0*x0 ;
383  x03 = x02*x0 ;
384  x04 = x03*x0 ;
385  x05 = x04*x0;
386 
387  c1 = (x2 - x1)/x1/x2 ;
388  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
389  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
390 
391  result -= (a1/x02 + a3/x04)*xln1 ;
392  result -= (a2/x02 + a4/x04)*c1 ;
393  result -= a3*c2/2/x02 ;
394  result -= a4*c3/3/x02 ;
395 
396  cof1 = a1/x02 + a3/x04 ;
397  cof2 = a2/x03 + a4/x05 ;
398 
399  result += 0.5*(cof1 +cof2)*xln2 ;
400  result += 0.5*(cof1 - cof2)*xln3 ;
401  }
402  result *= 2*hbarc/pi ;
403 
404  return result ;
405 
406 } // end of RePartDielectricConst
407 
409 //
410 // PAI differential cross-section in terms of
411 // simplified Allison's equation
412 //
413 
415 {
417  G4double betaGammaSq = fBetaGammaSq;
418  G4double integralTerm = IntegralTerm(omega);
419  G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
420  G4double epsilonRe = RePartDielectricConst(omega);
421  G4double epsilonIm = ImPartDielectricConst(i,omega);
422  G4double be4 ;
423  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
424  static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
425  be2 = betaGammaSq/(1 + betaGammaSq) ;
426  be4 = be2*be2 ;
427 
428  cof = 1 ;
429  x1 = log(2*electron_mass_c2/omega) ;
430 
431  if( betaGammaSq < 0.01 ) x2 = log(be2) ;
432  else
433  {
434  x2 = -log( (1/betaGammaSq - epsilonRe)*
435  (1/betaGammaSq - epsilonRe) +
436  epsilonIm*epsilonIm )/2 ;
437  }
438  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
439  {
440  x6=0 ;
441  }
442  else
443  {
444  x3 = -epsilonRe + 1/betaGammaSq ;
445  x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
446  epsilonIm*epsilonIm) ;
447 
448  x7 = atan2(epsilonIm,x3) ;
449  x6 = x5 * x7 ;
450  }
451  // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
452 
453  x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
454  // if( x4 < 0.0 ) x4 = 0.0 ;
455  x8 = (1 + epsilonRe)*(1 + epsilonRe) +
456  epsilonIm*epsilonIm;
457 
458  result = (x4 + cof*integralTerm/omega/omega) ;
459  if(result < 1.0e-8) result = 1.0e-8 ;
460  result *= fine_structure_const/be2/pi ;
461  // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
462  // result *= (1-exp(-be2/betaBohr2)) ;
463  result *= (1-exp(-be4/betaBohr4)) ;
464  if(fDensity >= fSolidDensity)
465  {
466  result /= x8 ;
467  }
468  return result ;
469 
470 } // end of DifPAIxSection
471 
473 //
474 // Differential PAI dEdx(omega)=omega*dNdx(omega)
475 //
476 
478 {
479  G4double dEdx = omega*DifPAIxSection(omega);
480  return dEdx;
481 }
482 
484 //
485 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
486 
488 {
490  G4double betaGammaSq = fBetaGammaSq;
491  G4double epsilonRe = RePartDielectricConst(omega);
492  G4double epsilonIm = ImPartDielectricConst(i,omega);
493 
494  G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ;
495  G4double be2, be4;
496 
497  //cof = 1.0 ;
498  static const G4double cofBetaBohr = 4.0 ;
499  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
500  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
501 
502  be2 = betaGammaSq/(1 + betaGammaSq) ;
503  be4 = be2*be2 ;
504 
505  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
506  else
507  {
508  logarithm = -log( (1/betaGammaSq - epsilonRe)*
509  (1/betaGammaSq - epsilonRe) +
510  epsilonIm*epsilonIm )*0.5 ;
511  logarithm += log(1+1.0/betaGammaSq) ;
512  }
513 
514  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
515  {
516  argument = 0.0 ;
517  }
518  else
519  {
520  x3 = -epsilonRe + 1.0/betaGammaSq ;
521  x5 = -1.0 - epsilonRe +
522  be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
523  epsilonIm*epsilonIm) ;
524  if( x3 == 0.0 ) argument = 0.5*pi;
525  else argument = atan2(epsilonIm,x3) ;
526  argument *= x5 ;
527  }
528  dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
529 
530  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
531 
532  dNdxC *= fine_structure_const/be2/pi ;
533 
534  dNdxC *= (1-exp(-be4/betaBohr4)) ;
535 
536  if(fDensity >= fSolidDensity)
537  {
538  modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
539  epsilonIm*epsilonIm;
540  dNdxC /= modul2 ;
541  }
542  return dNdxC ;
543 
544 } // end of PAIdNdxCerenkov
545 
547 //
548 // Calculation od dN/dx of collisions with creation of longitudinal EM
549 // excitations (plasmons, delta-electrons)
550 
552 {
554  G4double betaGammaSq = fBetaGammaSq;
555  G4double integralTerm = IntegralTerm(omega);
556  G4double epsilonRe = RePartDielectricConst(omega);
557  G4double epsilonIm = ImPartDielectricConst(i,omega);
558 
559  G4double cof, resonance, modul2, dNdxP ;
560  G4double be2, be4;
561 
562  cof = 1 ;
563  static const G4double cofBetaBohr = 4.0 ;
564  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
565  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
566 
567  be2 = betaGammaSq/(1 + betaGammaSq) ;
568  be4 = be2*be2 ;
569 
570  resonance = log(2*electron_mass_c2*be2/omega) ;
571  resonance *= epsilonIm/hbarc ;
572 
573 
574  dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
575 
576  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
577 
578  dNdxP *= fine_structure_const/be2/pi ;
579  dNdxP *= (1-exp(-be4/betaBohr4)) ;
580 
581  if( fDensity >= fSolidDensity )
582  {
583  modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
584  epsilonIm*epsilonIm;
585  dNdxP /= modul2 ;
586  }
587  return dNdxP ;
588 
589 } // end of PAIdNdxPlasmon
590 
592 //
593 // Calculation of the PAI integral cross-section
594 // = specific primary ionisation, 1/cm
595 //
596 
598 {
599  G4int i,k,i1,i2;
600  G4double energy1, energy2, result = 0.;
601 
602  fBetaGammaSq = bg2;
603  fTmax = Tmax;
604 
605  if(fPAIxscVector) delete fPAIxscVector;
606 
608  fPAIxscVector->PutValue(fPAIbin-1,result);
609 
610  for( i = fIntervalNumber - 1; i >= 0; i-- )
611  {
612  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
613  }
614  if (i < 0) i = 0; // Tmax should be more than
615  // first ionisation potential
616  fIntervalTmax = i;
617 
619 
620  for( k = fPAIbin - 2; k >= 0; k-- )
621  {
622  energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
623  energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
624 
625  for( i = fIntervalTmax; i >= 0; i-- )
626  {
627  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
628  }
629  if(i < 0) i = 0;
630  i2 = i;
631 
632  for( i = fIntervalTmax; i >= 0; i-- )
633  {
634  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
635  }
636  if(i < 0) i = 0;
637  i1 = i;
638 
639  if( i1 == i2 )
640  {
641  fCurrentInterval = i1;
642  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
643  energy1,energy2);
644  fPAIxscVector->PutValue(k,result);
645  }
646  else
647  {
648  for( i = i2; i >= i1; i-- )
649  {
650  fCurrentInterval = i;
651 
652  if( i==i2 ) result += integral.Legendre10(this,
654  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
655 
656  else if( i == i1 ) result += integral.Legendre10(this,
658  (*(*fMatSandiaMatrix)[i+1])[0]);
659 
660  else result += integral.Legendre10(this,
662  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
663  }
664  fPAIxscVector->PutValue(k,result);
665  }
666  // G4cout<<k<<"\t"<<result<<G4endl;
667  }
668  return ;
669 }
670 
671 
673 //
674 // Calculation of the PAI integral dEdx
675 // = mean energy loss per unit length, keV/cm
676 //
677 
679 {
680  G4int i,k,i1,i2;
681  G4double energy1, energy2, result = 0.;
682 
683  fBetaGammaSq = bg2;
684  fTmax = Tmax;
685 
686  if(fPAIdEdxVector) delete fPAIdEdxVector;
687 
689  fPAIdEdxVector->PutValue(fPAIbin-1,result);
690 
691  for( i = fIntervalNumber - 1; i >= 0; i-- )
692  {
693  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
694  }
695  if (i < 0) i = 0; // Tmax should be more than
696  // first ionisation potential
697  fIntervalTmax = i;
698 
700 
701  for( k = fPAIbin - 2; k >= 0; k-- )
702  {
703  energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
704  energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
705 
706  for( i = fIntervalTmax; i >= 0; i-- )
707  {
708  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
709  }
710  if(i < 0) i = 0;
711  i2 = i;
712 
713  for( i = fIntervalTmax; i >= 0; i-- )
714  {
715  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
716  }
717  if(i < 0) i = 0;
718  i1 = i;
719 
720  if( i1 == i2 )
721  {
722  fCurrentInterval = i1;
723  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
724  energy1,energy2);
725  fPAIdEdxVector->PutValue(k,result);
726  }
727  else
728  {
729  for( i = i2; i >= i1; i-- )
730  {
731  fCurrentInterval = i;
732 
733  if( i==i2 ) result += integral.Legendre10(this,
735  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
736 
737  else if( i == i1 ) result += integral.Legendre10(this,
738  &G4InitXscPAI::DifPAIdEdx,energy1,
739  (*(*fMatSandiaMatrix)[i+1])[0]);
740 
741  else result += integral.Legendre10(this,
743  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
744  }
745  fPAIdEdxVector->PutValue(k,result);
746  }
747  // G4cout<<k<<"\t"<<result<<G4endl;
748  }
749  return ;
750 }
751 
753 //
754 // Calculation of the PAI Cerenkov integral cross-section
755 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
756 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
757 
759 {
760  G4int i,k,i1,i2;
761  G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
762 
763  fBetaGammaSq = bg2;
764  fTmax = Tmax;
765  beta2 = bg2/(1+bg2);
766 
768  if(fChCosSqVector) delete fChCosSqVector;
769  if(fChWidthVector) delete fChWidthVector;
770 
774 
775  fPAIphotonVector->PutValue(fPAIbin-1,result);
778 
779  for( i = fIntervalNumber - 1; i >= 0; i-- )
780  {
781  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
782  }
783  if (i < 0) i = 0; // Tmax should be more than
784  // first ionisation potential
785  fIntervalTmax = i;
786 
788 
789  for( k = fPAIbin - 2; k >= 0; k-- )
790  {
791  energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
792  energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
793 
794  for( i = fIntervalTmax; i >= 0; i-- )
795  {
796  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
797  }
798  if(i < 0) i = 0;
799  i2 = i;
800 
801  for( i = fIntervalTmax; i >= 0; i-- )
802  {
803  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
804  }
805  if(i < 0) i = 0;
806  i1 = i;
807 
808  module2 = ModuleSqDielectricConst(i1,energy1);
809  cos2 = RePartDielectricConst(energy1)/module2/beta2;
810  width = ImPartDielectricConst(i1,energy1)/module2/beta2;
811 
812  fChCosSqVector->PutValue(k,cos2);
813  fChWidthVector->PutValue(k,width);
814 
815  if( i1 == i2 )
816  {
817  fCurrentInterval = i1;
818  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
819  energy1,energy2);
820  fPAIphotonVector->PutValue(k,result);
821 
822  }
823  else
824  {
825  for( i = i2; i >= i1; i-- )
826  {
827  fCurrentInterval = i;
828 
829  if( i==i2 ) result += integral.Legendre10(this,
831  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
832 
833  else if( i == i1 ) result += integral.Legendre10(this,
835  (*(*fMatSandiaMatrix)[i+1])[0]);
836 
837  else result += integral.Legendre10(this,
839  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
840  }
841  fPAIphotonVector->PutValue(k,result);
842  }
843  // G4cout<<k<<"\t"<<result<<G4endl;
844  }
845  return;
846 } // end of IntegralCerenkov
847 
849 //
850 // Calculation of the PAI Plasmon integral cross-section
851 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
852 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
853 
855 {
856  G4int i,k,i1,i2;
857  G4double energy1, energy2, result = 0.;
858 
859  fBetaGammaSq = bg2;
860  fTmax = Tmax;
861 
863 
866 
867  for( i = fIntervalNumber - 1; i >= 0; i-- )
868  {
869  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
870  }
871  if (i < 0) i = 0; // Tmax should be more than
872  // first ionisation potential
873  fIntervalTmax = i;
874 
876 
877  for( k = fPAIbin - 2; k >= 0; k-- )
878  {
879  energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
880  energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
881 
882  for( i = fIntervalTmax; i >= 0; i-- )
883  {
884  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
885  }
886  if(i < 0) i = 0;
887  i2 = i;
888 
889  for( i = fIntervalTmax; i >= 0; i-- )
890  {
891  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
892  }
893  if(i < 0) i = 0;
894  i1 = i;
895 
896  if( i1 == i2 )
897  {
898  fCurrentInterval = i1;
899  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
900  energy1,energy2);
901  fPAIelectronVector->PutValue(k,result);
902  }
903  else
904  {
905  for( i = i2; i >= i1; i-- )
906  {
907  fCurrentInterval = i;
908 
909  if( i==i2 ) result += integral.Legendre10(this,
911  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
912 
913  else if( i == i1 ) result += integral.Legendre10(this,
915  (*(*fMatSandiaMatrix)[i+1])[0]);
916 
917  else result += integral.Legendre10(this,
919  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
920  }
921  fPAIelectronVector->PutValue(k,result);
922  }
923  // G4cout<<k<<"\t"<<result<<G4endl;
924  }
925  return;
926 } // end of IntegralPlasmon
927 
928 
930 //
931 //
932 
934 {
935  G4int i ;
936  G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
937 
938  omega2 = omega*omega ;
939  omega3 = omega2*omega ;
940  omega4 = omega2*omega2 ;
941 
942  for(i = 0; i < fIntervalNumber;i++)
943  {
944  if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
945  }
946  if( i == 0 )
947  {
948  G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
949  }
950  else i-- ;
951 
952  a1 = (*(*fMatSandiaMatrix)[i])[1];
953  a2 = (*(*fMatSandiaMatrix)[i])[2];
954  a3 = (*(*fMatSandiaMatrix)[i])[3];
955  a4 = (*(*fMatSandiaMatrix)[i])[4];
956 
957  lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
958 
959  return lambda ;
960 }
961 
963 //
964 //
965 
967 //
968 //
969 
971 {
972  G4double loss = 0.0 ;
973  loss *= step;
974 
975  return loss ;
976 }
977 
979 //
980 //
981 
983 {
984  G4double loss = 0.0 ;
985  loss *= step;
986 
987  return loss ;
988 }
989 
991 //
992 //
993 
995 {
996 
997 
998  G4double loss = 0.0 ;
999  loss *= step;
1000  return loss ;
1001 }
1002 
1003 
1004 //
1005 // end of G4InitXscPAI implementation file
1006 //
1008 
G4double fTmax
G4PhysicsLogVector * fPAIxscVector
G4int fCurrentInterval
G4PhysicsLogVector * fPAIelectronVector
G4int fIntervalTmax
G4double GetLowEdgeEnergy(size_t binNumber) const
G4OrderedTable * fMatSandiaMatrix
G4double ModuleSqDielectricConst(G4int intervalNumber, G4double energy)
static const G4double fSolidDensity
Float_t x1[n_points_granero]
Definition: compare.C:5
size_t GetIndex() const
Definition: G4Material.hh:262
G4double DifPAIdEdx(G4double omega)
#define G4endl
Definition: G4ios.hh:61
static constexpr double hbarc
G4double fDensity
G4PhysicsLogVector * fPAIphotonVector
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double GetStepPlasmonLoss(G4double step)
G4int GetMaxInterval() const
void IntegralPAIdEdx(G4double bg2, G4double Tmax)
void IntegralPlasmon(G4double bg2, G4double Tmax)
G4int fIntervalNumber
void IntegralPAIxSection(G4double bg2, G4double Tmax)
static const G4int fPAIbin
static constexpr double g
Definition: G4SIunits.hh:183
G4PhysicsLogVector * fChWidthVector
double G4double
Definition: G4Types.hh:76
#define width
static const G4double fDelta
G4PhysicsLogVector * fPAIdEdxVector
TCanvas * c2
Definition: plot_hist.C:75
static constexpr double electron_mass_c2
G4PhysicsLogVector * fChCosSqVector
G4double PAIdNdxPlasmon(G4double omega)
G4double fNormalizationCof
G4double GetStepEnergyLoss(G4double step)
G4double DifPAIxSection(G4double omega)
G4double RePartDielectricConst(G4double energy)
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
virtual ~G4InitXscPAI()
G4double GetStepCerenkovLoss(G4double step)
G4InitXscPAI(const G4MaterialCutsCouple *matCC)
Definition: G4InitXscPAI.cc:70
G4double fElectronDensity
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
void KillCloseIntervals()
G4GLOB_DLL std::ostream G4cout
void Normalisation()
G4double fBetaGammaSq
G4SandiaTable * fSandia
const G4Material * GetMaterial() const
static constexpr double pi
Definition: G4SIunits.hh:75
G4double GetSandiaMatTable(G4int, G4int) const
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4double PAIdNdxCherenkov(G4double omega)
static constexpr double cm3
Definition: G4SIunits.hh:121
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
void IntegralCherenkov(G4double bg2, G4double Tmax)
G4double GetElectronDensity() const
Definition: G4Material.hh:218
G4double IntegralTerm(G4double omega)
void PutValue(size_t index, G4double theValue)
G4double GetPhotonLambda(G4double omega)
G4double GetDensity() const
Definition: G4Material.hh:181