Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InitXscPAI.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 // $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 
70 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC)
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();
81  fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
82  matIndex = matCC->GetMaterial()->GetIndex();
83 
84  fSandia = new G4SandiaTable(matIndex);
85  fIntervalNumber = fSandia->GetMaxInterval()-1;
86 
87  fMatSandiaMatrix = new G4OrderedTable();
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  }
102  KillCloseIntervals();
103  Normalisation();
104  fBetaGammaSq = fTmax = 0.0;
105  fIntervalTmax = fCurrentInterval = 0;
106 }
107 
108 
109 
110 
112 //
113 // Destructor
114 
115 G4InitXscPAI::~G4InitXscPAI()
116 {
117  if(fPAIxscVector) delete fPAIxscVector;
118  if(fPAIdEdxVector) delete fPAIdEdxVector;
119  if(fPAIphotonVector) delete fPAIphotonVector;
120  if(fPAIelectronVector) delete fPAIelectronVector;
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 
131 void G4InitXscPAI::KillCloseIntervals()
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 
162 void G4InitXscPAI::Normalisation()
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  }
181  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
182  fNormalizationCof *= fElectronDensity;
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 
232 G4double G4InitXscPAI::RutherfordIntegral( G4int k,
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 
256 G4double G4InitXscPAI::IntegralTerm(G4double omega)
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 
294 G4double G4InitXscPAI::ImPartDielectricConst( G4int k ,
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 
320 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int k ,
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 
343 G4double G4InitXscPAI::RePartDielectricConst(G4double enb)
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 
414 G4double G4InitXscPAI::DifPAIxSection( G4double omega )
415 {
416  G4int i = fCurrentInterval;
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 
477 G4double G4InitXscPAI::DifPAIdEdx( G4double omega )
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 
487 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega )
488 {
489  G4int i = fCurrentInterval;
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 
551 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega )
552 {
553  G4int i = fCurrentInterval;
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 
597 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax)
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 
607  fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
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 
618  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
653  &G4InitXscPAI::DifPAIxSection,
654  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
655 
656  else if( i == i1 ) result += integral.Legendre10(this,
657  &G4InitXscPAI::DifPAIxSection,energy1,
658  (*(*fMatSandiaMatrix)[i+1])[0]);
659 
660  else result += integral.Legendre10(this,
661  &G4InitXscPAI::DifPAIxSection,
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 
678 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax)
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 
688  fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
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 
699  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
734  &G4InitXscPAI::DifPAIdEdx,
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,
742  &G4InitXscPAI::DifPAIdEdx,
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 
758 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax)
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 
767  if(fPAIphotonVector) delete fPAIphotonVector;
768  if(fChCosSqVector) delete fChCosSqVector;
769  if(fChWidthVector) delete fChWidthVector;
770 
771  fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
772  fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
773  fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
774 
775  fPAIphotonVector->PutValue(fPAIbin-1,result);
776  fChCosSqVector->PutValue(fPAIbin-1,1.);
777  fChWidthVector->PutValue(fPAIbin-1,1e-7);
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 
787  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
830  &G4InitXscPAI::PAIdNdxCherenkov,
831  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
832 
833  else if( i == i1 ) result += integral.Legendre10(this,
834  &G4InitXscPAI::PAIdNdxCherenkov,energy1,
835  (*(*fMatSandiaMatrix)[i+1])[0]);
836 
837  else result += integral.Legendre10(this,
838  &G4InitXscPAI::PAIdNdxCherenkov,
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 
854 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax)
855 {
856  G4int i,k,i1,i2;
857  G4double energy1, energy2, result = 0.;
858 
859  fBetaGammaSq = bg2;
860  fTmax = Tmax;
861 
862  if(fPAIelectronVector) delete fPAIelectronVector;
863 
864  fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
865  fPAIelectronVector->PutValue(fPAIbin-1,result);
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 
875  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
910  &G4InitXscPAI::PAIdNdxPlasmon,
911  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
912 
913  else if( i == i1 ) result += integral.Legendre10(this,
914  &G4InitXscPAI::PAIdNdxPlasmon,energy1,
915  (*(*fMatSandiaMatrix)[i+1])[0]);
916 
917  else result += integral.Legendre10(this,
918  &G4InitXscPAI::PAIdNdxPlasmon,
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 
933 G4double G4InitXscPAI::GetPhotonLambda( G4double omega )
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 
970 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step )
971 {
972  G4double loss = 0.0 ;
973  loss *= step;
974 
975  return loss ;
976 }
977 
979 //
980 //
981 
982 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step )
983 {
984  G4double loss = 0.0 ;
985  loss *= step;
986 
987  return loss ;
988 }
989 
991 //
992 //
993 
994 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step )
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 
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4Integrator::Legendre10
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
c2
static c2_factory< G4double > c2
Definition: G4ScreenedNuclearRecoil.cc:126
G4MaterialCutsCouple.hh
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4InitXscPAI::PAIdNdxCherenkov
G4double PAIdNdxCherenkov(G4double omega)
Definition: G4InitXscPAI.cc:487
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:180
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.cc:716
G4InitXscPAI::ModuleSqDielectricConst
G4double ModuleSqDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4InitXscPAI.cc:320
width
#define width
G4InitXscPAI::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4InitXscPAI.cc:294
G4InitXscPAI::KillCloseIntervals
void KillCloseIntervals()
Definition: G4InitXscPAI.cc:131
G4InitXscPAI::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition: G4InitXscPAI.cc:343
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:151
G4int
int G4int
Definition: G4Types.hh:78
G4InitXscPAI::GetStepCerenkovLoss
G4double GetStepCerenkovLoss(G4double step)
Definition: G4InitXscPAI.cc:982
G4InitXscPAI.hh
G4InitXscPAI::DifPAIxSection
G4double DifPAIxSection(G4double omega)
Definition: G4InitXscPAI.cc:414
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4InitXscPAI::IntegralCherenkov
void IntegralCherenkov(G4double bg2, G4double Tmax)
Definition: G4InitXscPAI.cc:758
G4cout
G4GLOB_DLL std::ostream G4cout
G4InitXscPAI::Normalisation
void Normalisation()
Definition: G4InitXscPAI.cc:162
G4InitXscPAI::GetStepEnergyLoss
G4double GetStepEnergyLoss(G4double step)
Definition: G4InitXscPAI.cc:970
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4Poisson.hh
G4Material.hh
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4InitXscPAI::G4InitXscPAI
G4InitXscPAI(const G4MaterialCutsCouple *matCC)
Definition: G4InitXscPAI.cc:70
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4InitXscPAI::IntegralPAIxSection
void IntegralPAIxSection(G4double bg2, G4double Tmax)
Definition: G4InitXscPAI.cc:597
cm3
static constexpr double cm3
Definition: G4SIunits.hh:121
G4PhysicalConstants.hh
G4InitXscPAI::DifPAIdEdx
G4double DifPAIdEdx(G4double omega)
Definition: G4InitXscPAI.cc:477
G4InitXscPAI::IntegralPlasmon
void IntegralPlasmon(G4double bg2, G4double Tmax)
Definition: G4InitXscPAI.cc:854
G4ios.hh
G4InitXscPAI::~G4InitXscPAI
virtual ~G4InitXscPAI()
Definition: G4InitXscPAI.cc:115
G4Integrator.hh
G4SandiaTable::GetSandiaMatTable
G4double GetSandiaMatTable(G4int, G4int) const
Definition: G4SandiaTable.cc:1016
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4InitXscPAI::IntegralPAIdEdx
void IntegralPAIdEdx(G4double bg2, G4double Tmax)
Definition: G4InitXscPAI.cc:678
G4SandiaTable.hh
G4InitXscPAI::RutherfordIntegral
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
Definition: G4InitXscPAI.cc:232
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4InitXscPAI::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4double omega)
Definition: G4InitXscPAI.cc:551
plottest35.c1
tuple c1
Definition: plottest35.py:14
G4InitXscPAI::IntegralTerm
G4double IntegralTerm(G4double omega)
Definition: G4InitXscPAI.cc:256
G4InitXscPAI::GetPhotonLambda
G4double GetPhotonLambda(G4double omega)
Definition: G4InitXscPAI.cc:933
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4InitXscPAI::GetStepPlasmonLoss
G4double GetStepPlasmonLoss(G4double step)
Definition: G4InitXscPAI.cc:994