ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InitXscPAI.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 //
28 // G4InitXscPAI.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: Vladimir.Grichine@cern.ch
36 //
37 // History:
38 //
39 
40 
41 
42 #include "G4InitXscPAI.hh"
43 
44 #include "globals.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4ios.hh"
48 #include "G4Poisson.hh"
49 #include "G4Integrator.hh"
50 #include "G4Material.hh"
51 #include "G4MaterialCutsCouple.hh"
52 #include "G4SandiaTable.hh"
53 
54 
55 
56 // Local class constants
57 
58 const G4double G4InitXscPAI::fDelta = 0.005 ; // energy shift from interval border
59 const G4int G4InitXscPAI::fPAIbin = 100 ; // size of energy transfer vectors
60 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border
61 
63 //
64 // Constructor
65 //
66 
67 using namespace std;
68 
70  : fPAIxscVector(nullptr),
71  fPAIdEdxVector(nullptr),
72  fPAIphotonVector(nullptr),
73  fPAIelectronVector(nullptr),
74  fChCosSqVector(nullptr),
75  fChWidthVector(nullptr)
76 {
77  G4int i, j, matIndex;
78 
79  fDensity = matCC->GetMaterial()->GetDensity();
81  matIndex = matCC->GetMaterial()->GetIndex();
82 
83  fSandia = new G4SandiaTable(matIndex);
85 
87 
88  for (i = 0; i < fIntervalNumber; i++)
89  {
90  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
91  }
92  for (i = 0; i < fIntervalNumber; i++)
93  {
94  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
95 
96  for(j = 1; j < 5 ; j++)
97  {
98  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
99  }
100  }
102  Normalisation();
103  fBetaGammaSq = fTmax = 0.0;
105 }
106 
107 
108 
109 
111 //
112 // Destructor
113 
115 {
116  if(fPAIxscVector) delete fPAIxscVector;
117  if(fPAIdEdxVector) delete fPAIdEdxVector;
120  if(fChCosSqVector) delete fChCosSqVector;
121  if(fChWidthVector) delete fChWidthVector;
122  delete fSandia;
123  delete fMatSandiaMatrix;
124 }
125 
127 //
128 // Kill close intervals, recalculate fIntervalNumber
129 
131 {
132  G4int i, j, k;
133  G4double energy1, energy2;
134 
135  for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
136  {
137  energy1 = (*(*fMatSandiaMatrix)[i])[0];
138  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
139 
140  if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ;
141  else
142  {
143  for(j = i; j < fIntervalNumber-1; j++)
144  {
145  for( k = 0; k < 5; k++ )
146  {
147  (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
148  }
149  }
150  fIntervalNumber-- ;
151  i-- ;
152  }
153  }
154 
155 }
156 
158 //
159 // Kill close intervals, recalculate fIntervalNumber
160 
162 {
163  G4int i, j;
164  G4double energy1, energy2, /*delta,*/ cof; // , shift;
165 
166  energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
167  energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
168 
169 
170  cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
171 
172  for( i = fIntervalNumber-2; i >= 0; i-- )
173  {
174  energy1 = (*(*fMatSandiaMatrix)[i])[0];
175  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
176 
177  cof += RutherfordIntegral(i,energy1,energy2);
178  // G4cout<<"norm. cof = "<<cof<<G4endl;
179  }
182  //delta = fNormalizationCof - cof;
183  fNormalizationCof /= cof;
184  // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
185  // <<"; at delta ="<<delta<<G4endl ;
186 
187  for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
188  {
189  for(j = 1; j < 5 ; j++)
190  {
191  (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
192  }
193  }
194  /*
195  if(delta > 0) // shift the first energy interval
196  {
197  for(i=1;i<100;i++)
198  {
199  energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
200  energy2 = (*(*fMatSandiaMatrix)[0])[0];
201  shift = RutherfordIntegral(0,energy1,energy2);
202  G4cout<<shift<<"\t";
203  if(shift >= delta) break;
204  }
205  (*(*fMatSandiaMatrix)[0])[0] = energy1;
206  cof += shift;
207  }
208  else if(delta < 0)
209  {
210  for(i=1;i<100;i++)
211  {
212  energy1 = (*(*fMatSandiaMatrix)[0])[0];
213  energy2 = (*(*fMatSandiaMatrix)[0])[0] +
214  ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
215  shift = RutherfordIntegral(0,energy1,energy2);
216  if( shift >= std::abs(delta) ) break;
217  }
218  (*(*fMatSandiaMatrix)[0])[0] = energy2;
219  cof -= shift;
220  }
221  G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
222  <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ;
223  */
224 }
225 
227 //
228 // Integration over electrons that could be considered
229 // quasi-free at energy transfer of interest
230 
232  G4double x1,
233  G4double x2 )
234 {
235  G4double c1, c2, c3, a1, a2, a3, a4 ;
236 
237  a1 = (*(*fMatSandiaMatrix)[k])[1];
238  a2 = (*(*fMatSandiaMatrix)[k])[2];
239  a3 = (*(*fMatSandiaMatrix)[k])[3];
240  a4 = (*(*fMatSandiaMatrix)[k])[4];
241  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
242  c1 = (x2 - x1)/x1/x2 ;
243  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
244  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
245  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
246 
247  return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
248 
249 } // end of RutherfordIntegral
250 
252 //
253 // Integrate photo-absorption cross-section from I1 up to omega
254 
256 {
257  G4int i;
258  G4double energy1, energy2, result = 0.;
259 
260  for( i = 0; i <= fIntervalTmax; i++ )
261  {
262  if(i == fIntervalTmax)
263  {
264  energy1 = (*(*fMatSandiaMatrix)[i])[0];
265  result += RutherfordIntegral(i,energy1,omega);
266  }
267  else
268  {
269  if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
270  {
271  energy1 = (*(*fMatSandiaMatrix)[i])[0];
272  result += RutherfordIntegral(i,energy1,omega);
273  break;
274  }
275  else
276  {
277  energy1 = (*(*fMatSandiaMatrix)[i])[0];
278  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
279  result += RutherfordIntegral(i,energy1,energy2);
280  }
281  }
282  // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
283  }
284  return result;
285 }
286 
287 
289 //
290 // Imaginary part of dielectric constant
291 // (G4int k - interval number, G4double en1 - energy point)
292 
294  G4double energy1 )
295 {
296  G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
297 
298  a1 = (*(*fMatSandiaMatrix)[k])[1];
299  a2 = (*(*fMatSandiaMatrix)[k])[2];
300  a3 = (*(*fMatSandiaMatrix)[k])[3];
301  a4 = (*(*fMatSandiaMatrix)[k])[4];
302 
303  energy2 = energy1*energy1;
304  energy3 = energy2*energy1;
305  energy4 = energy3*energy1;
306 
307  result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
308  result *= hbarc/energy1 ;
309 
310  return result ;
311 
312 } // end of ImPartDielectricConst
313 
315 //
316 // Modulus squared of dielectric constant
317 // (G4int k - interval number, G4double omega - energy point)
318 
320  G4double omega )
321 {
322  G4double eIm2, eRe2, result;
323 
324  result = ImPartDielectricConst(k,omega);
325  eIm2 = result*result;
326 
327  result = RePartDielectricConst(omega);
328  eRe2 = result*result;
329 
330  result = eIm2 + eRe2;
331 
332  return result ;
333 }
334 
335 
337 //
338 // Real part of dielectric constant minus unit: epsilon_1 - 1
339 // (G4double enb - energy point)
340 //
341 
343 {
344  G4int i;
345  G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
346  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
347 
348  x0 = enb ;
349  result = 0 ;
350 
351  for( i = 0; i < fIntervalNumber-1; i++)
352  {
353  x1 = (*(*fMatSandiaMatrix)[i])[0];
354  x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
355 
356  a1 = (*(*fMatSandiaMatrix)[i])[1];
357  a2 = (*(*fMatSandiaMatrix)[i])[2];
358  a3 = (*(*fMatSandiaMatrix)[i])[3];
359  a4 = (*(*fMatSandiaMatrix)[i])[4];
360 
361  if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta )
362  {
363  if(x0 >= x1) x0 = x1*(1+fDelta);
364  else x0 = x1*(1-fDelta);
365  }
366  if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta )
367  {
368  if(x0 >= x2) x0 = x2*(1+fDelta);
369  else x0 = x2*(1-fDelta);
370  }
371  xx1 = x1 - x0 ;
372  xx2 = x2 - x0 ;
373  xx12 = xx2/xx1 ;
374 
375  if( xx12 < 0 ) xx12 = -xx12;
376 
377  xln1 = log(x2/x1) ;
378  xln2 = log(xx12) ;
379  xln3 = log((x2 + x0)/(x1 + x0)) ;
380 
381  x02 = x0*x0 ;
382  x03 = x02*x0 ;
383  x04 = x03*x0 ;
384  x05 = x04*x0;
385 
386  c1 = (x2 - x1)/x1/x2 ;
387  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
388  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
389 
390  result -= (a1/x02 + a3/x04)*xln1 ;
391  result -= (a2/x02 + a4/x04)*c1 ;
392  result -= a3*c2/2/x02 ;
393  result -= a4*c3/3/x02 ;
394 
395  cof1 = a1/x02 + a3/x04 ;
396  cof2 = a2/x03 + a4/x05 ;
397 
398  result += 0.5*(cof1 +cof2)*xln2 ;
399  result += 0.5*(cof1 - cof2)*xln3 ;
400  }
401  result *= 2*hbarc/pi ;
402 
403  return result ;
404 
405 } // end of RePartDielectricConst
406 
408 //
409 // PAI differential cross-section in terms of
410 // simplified Allison's equation
411 //
412 
414 {
416  G4double betaGammaSq = fBetaGammaSq;
417  G4double integralTerm = IntegralTerm(omega);
418  G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
419  G4double epsilonRe = RePartDielectricConst(omega);
420  G4double epsilonIm = ImPartDielectricConst(i,omega);
421  G4double be4 ;
422  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
423  static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
424  be2 = betaGammaSq/(1 + betaGammaSq) ;
425  be4 = be2*be2 ;
426 
427  cof = 1 ;
428  x1 = log(2*electron_mass_c2/omega) ;
429 
430  if( betaGammaSq < 0.01 ) x2 = log(be2) ;
431  else
432  {
433  x2 = -log( (1/betaGammaSq - epsilonRe)*
434  (1/betaGammaSq - epsilonRe) +
435  epsilonIm*epsilonIm )/2 ;
436  }
437  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
438  {
439  x6=0 ;
440  }
441  else
442  {
443  x3 = -epsilonRe + 1/betaGammaSq ;
444  x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
445  epsilonIm*epsilonIm) ;
446 
447  x7 = atan2(epsilonIm,x3) ;
448  x6 = x5 * x7 ;
449  }
450  // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
451 
452  x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
453  // if( x4 < 0.0 ) x4 = 0.0 ;
454  x8 = (1 + epsilonRe)*(1 + epsilonRe) +
455  epsilonIm*epsilonIm;
456 
457  result = (x4 + cof*integralTerm/omega/omega) ;
458  if(result < 1.0e-8) result = 1.0e-8 ;
459  result *= fine_structure_const/be2/pi ;
460  // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
461  // result *= (1-exp(-be2/betaBohr2)) ;
462  result *= (1-exp(-be4/betaBohr4)) ;
463  if(fDensity >= fSolidDensity)
464  {
465  result /= x8 ;
466  }
467  return result ;
468 
469 } // end of DifPAIxSection
470 
472 //
473 // Differential PAI dEdx(omega)=omega*dNdx(omega)
474 //
475 
477 {
478  G4double dEdx = omega*DifPAIxSection(omega);
479  return dEdx;
480 }
481 
483 //
484 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
485 
487 {
489  G4double betaGammaSq = fBetaGammaSq;
490  G4double epsilonRe = RePartDielectricConst(omega);
491  G4double epsilonIm = ImPartDielectricConst(i,omega);
492 
493  G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ;
494  G4double be2, be4;
495 
496  //cof = 1.0 ;
497  static const G4double cofBetaBohr = 4.0 ;
498  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
499  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
500 
501  be2 = betaGammaSq/(1 + betaGammaSq) ;
502  be4 = be2*be2 ;
503 
504  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
505  else
506  {
507  logarithm = -log( (1/betaGammaSq - epsilonRe)*
508  (1/betaGammaSq - epsilonRe) +
509  epsilonIm*epsilonIm )*0.5 ;
510  logarithm += log(1+1.0/betaGammaSq) ;
511  }
512 
513  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
514  {
515  argument = 0.0 ;
516  }
517  else
518  {
519  x3 = -epsilonRe + 1.0/betaGammaSq ;
520  x5 = -1.0 - epsilonRe +
521  be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
522  epsilonIm*epsilonIm) ;
523  if( x3 == 0.0 ) argument = 0.5*pi;
524  else argument = atan2(epsilonIm,x3) ;
525  argument *= x5 ;
526  }
527  dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
528 
529  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
530 
531  dNdxC *= fine_structure_const/be2/pi ;
532 
533  dNdxC *= (1-exp(-be4/betaBohr4)) ;
534 
535  if(fDensity >= fSolidDensity)
536  {
537  modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
538  epsilonIm*epsilonIm;
539  dNdxC /= modul2 ;
540  }
541  return dNdxC ;
542 
543 } // end of PAIdNdxCerenkov
544 
546 //
547 // Calculation od dN/dx of collisions with creation of longitudinal EM
548 // excitations (plasmons, delta-electrons)
549 
551 {
553  G4double betaGammaSq = fBetaGammaSq;
554  G4double integralTerm = IntegralTerm(omega);
555  G4double epsilonRe = RePartDielectricConst(omega);
556  G4double epsilonIm = ImPartDielectricConst(i,omega);
557 
558  G4double cof, resonance, modul2, dNdxP ;
559  G4double be2, be4;
560 
561  cof = 1 ;
562  static const G4double cofBetaBohr = 4.0 ;
563  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
564  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
565 
566  be2 = betaGammaSq/(1 + betaGammaSq) ;
567  be4 = be2*be2 ;
568 
569  resonance = log(2*electron_mass_c2*be2/omega) ;
570  resonance *= epsilonIm/hbarc ;
571 
572 
573  dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
574 
575  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
576 
577  dNdxP *= fine_structure_const/be2/pi ;
578  dNdxP *= (1-exp(-be4/betaBohr4)) ;
579 
580  if( fDensity >= fSolidDensity )
581  {
582  modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
583  epsilonIm*epsilonIm;
584  dNdxP /= modul2 ;
585  }
586  return dNdxP ;
587 
588 } // end of PAIdNdxPlasmon
589 
591 //
592 // Calculation of the PAI integral cross-section
593 // = specific primary ionisation, 1/cm
594 //
595 
597 {
598  G4int i,k,i1,i2;
599  G4double energy1, energy2, result = 0.;
600 
601  fBetaGammaSq = bg2;
602  fTmax = Tmax;
603 
604  if(fPAIxscVector) delete fPAIxscVector;
605 
607  fPAIxscVector->PutValue(fPAIbin-1,result);
608 
609  for( i = fIntervalNumber - 1; i >= 0; i-- )
610  {
611  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
612  }
613  if (i < 0) i = 0; // Tmax should be more than
614  // first ionisation potential
615  fIntervalTmax = i;
616 
618 
619  for( k = fPAIbin - 2; k >= 0; k-- )
620  {
621  energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
622  energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
623 
624  for( i = fIntervalTmax; i >= 0; i-- )
625  {
626  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
627  }
628  if(i < 0) i = 0;
629  i2 = i;
630 
631  for( i = fIntervalTmax; i >= 0; i-- )
632  {
633  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
634  }
635  if(i < 0) i = 0;
636  i1 = i;
637 
638  if( i1 == i2 )
639  {
640  fCurrentInterval = i1;
641  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
642  energy1,energy2);
643  fPAIxscVector->PutValue(k,result);
644  }
645  else
646  {
647  for( i = i2; i >= i1; i-- )
648  {
649  fCurrentInterval = i;
650 
651  if( i==i2 ) result += integral.Legendre10(this,
653  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
654 
655  else if( i == i1 ) result += integral.Legendre10(this,
657  (*(*fMatSandiaMatrix)[i+1])[0]);
658 
659  else result += integral.Legendre10(this,
661  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
662  }
663  fPAIxscVector->PutValue(k,result);
664  }
665  // G4cout<<k<<"\t"<<result<<G4endl;
666  }
667  return ;
668 }
669 
670 
672 //
673 // Calculation of the PAI integral dEdx
674 // = mean energy loss per unit length, keV/cm
675 //
676 
678 {
679  G4int i,k,i1,i2;
680  G4double energy1, energy2, result = 0.;
681 
682  fBetaGammaSq = bg2;
683  fTmax = Tmax;
684 
685  if(fPAIdEdxVector) delete fPAIdEdxVector;
686 
688  fPAIdEdxVector->PutValue(fPAIbin-1,result);
689 
690  for( i = fIntervalNumber - 1; i >= 0; i-- )
691  {
692  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
693  }
694  if (i < 0) i = 0; // Tmax should be more than
695  // first ionisation potential
696  fIntervalTmax = i;
697 
699 
700  for( k = fPAIbin - 2; k >= 0; k-- )
701  {
702  energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
703  energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
704 
705  for( i = fIntervalTmax; i >= 0; i-- )
706  {
707  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
708  }
709  if(i < 0) i = 0;
710  i2 = i;
711 
712  for( i = fIntervalTmax; i >= 0; i-- )
713  {
714  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
715  }
716  if(i < 0) i = 0;
717  i1 = i;
718 
719  if( i1 == i2 )
720  {
721  fCurrentInterval = i1;
722  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
723  energy1,energy2);
724  fPAIdEdxVector->PutValue(k,result);
725  }
726  else
727  {
728  for( i = i2; i >= i1; i-- )
729  {
730  fCurrentInterval = i;
731 
732  if( i==i2 ) result += integral.Legendre10(this,
734  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
735 
736  else if( i == i1 ) result += integral.Legendre10(this,
737  &G4InitXscPAI::DifPAIdEdx,energy1,
738  (*(*fMatSandiaMatrix)[i+1])[0]);
739 
740  else result += integral.Legendre10(this,
742  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
743  }
744  fPAIdEdxVector->PutValue(k,result);
745  }
746  // G4cout<<k<<"\t"<<result<<G4endl;
747  }
748  return ;
749 }
750 
752 //
753 // Calculation of the PAI Cerenkov integral cross-section
754 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
755 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
756 
758 {
759  G4int i,k,i1,i2;
760  G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
761 
762  fBetaGammaSq = bg2;
763  fTmax = Tmax;
764  beta2 = bg2/(1+bg2);
765 
767  if(fChCosSqVector) delete fChCosSqVector;
768  if(fChWidthVector) delete fChWidthVector;
769 
773 
774  fPAIphotonVector->PutValue(fPAIbin-1,result);
777 
778  for( i = fIntervalNumber - 1; i >= 0; i-- )
779  {
780  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
781  }
782  if (i < 0) i = 0; // Tmax should be more than
783  // first ionisation potential
784  fIntervalTmax = i;
785 
787 
788  for( k = fPAIbin - 2; k >= 0; k-- )
789  {
790  energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
791  energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
792 
793  for( i = fIntervalTmax; i >= 0; i-- )
794  {
795  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
796  }
797  if(i < 0) i = 0;
798  i2 = i;
799 
800  for( i = fIntervalTmax; i >= 0; i-- )
801  {
802  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
803  }
804  if(i < 0) i = 0;
805  i1 = i;
806 
807  module2 = ModuleSqDielectricConst(i1,energy1);
808  cos2 = RePartDielectricConst(energy1)/module2/beta2;
809  width = ImPartDielectricConst(i1,energy1)/module2/beta2;
810 
811  fChCosSqVector->PutValue(k,cos2);
812  fChWidthVector->PutValue(k,width);
813 
814  if( i1 == i2 )
815  {
816  fCurrentInterval = i1;
817  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
818  energy1,energy2);
819  fPAIphotonVector->PutValue(k,result);
820 
821  }
822  else
823  {
824  for( i = i2; i >= i1; i-- )
825  {
826  fCurrentInterval = i;
827 
828  if( i==i2 ) result += integral.Legendre10(this,
830  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
831 
832  else if( i == i1 ) result += integral.Legendre10(this,
834  (*(*fMatSandiaMatrix)[i+1])[0]);
835 
836  else result += integral.Legendre10(this,
838  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
839  }
840  fPAIphotonVector->PutValue(k,result);
841  }
842  // G4cout<<k<<"\t"<<result<<G4endl;
843  }
844  return;
845 } // end of IntegralCerenkov
846 
848 //
849 // Calculation of the PAI Plasmon integral cross-section
850 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
851 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
852 
854 {
855  G4int i,k,i1,i2;
856  G4double energy1, energy2, result = 0.;
857 
858  fBetaGammaSq = bg2;
859  fTmax = Tmax;
860 
862 
865 
866  for( i = fIntervalNumber - 1; i >= 0; i-- )
867  {
868  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
869  }
870  if (i < 0) i = 0; // Tmax should be more than
871  // first ionisation potential
872  fIntervalTmax = i;
873 
875 
876  for( k = fPAIbin - 2; k >= 0; k-- )
877  {
878  energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
879  energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
880 
881  for( i = fIntervalTmax; i >= 0; i-- )
882  {
883  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
884  }
885  if(i < 0) i = 0;
886  i2 = i;
887 
888  for( i = fIntervalTmax; i >= 0; i-- )
889  {
890  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
891  }
892  if(i < 0) i = 0;
893  i1 = i;
894 
895  if( i1 == i2 )
896  {
897  fCurrentInterval = i1;
898  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
899  energy1,energy2);
900  fPAIelectronVector->PutValue(k,result);
901  }
902  else
903  {
904  for( i = i2; i >= i1; i-- )
905  {
906  fCurrentInterval = i;
907 
908  if( i==i2 ) result += integral.Legendre10(this,
910  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
911 
912  else if( i == i1 ) result += integral.Legendre10(this,
914  (*(*fMatSandiaMatrix)[i+1])[0]);
915 
916  else result += integral.Legendre10(this,
918  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
919  }
920  fPAIelectronVector->PutValue(k,result);
921  }
922  // G4cout<<k<<"\t"<<result<<G4endl;
923  }
924  return;
925 } // end of IntegralPlasmon
926 
927 
929 //
930 //
931 
933 {
934  G4int i ;
935  G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
936 
937  omega2 = omega*omega ;
938  omega3 = omega2*omega ;
939  omega4 = omega2*omega2 ;
940 
941  for(i = 0; i < fIntervalNumber;i++)
942  {
943  if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
944  }
945  if( i == 0 )
946  {
947  G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
948  }
949  else i-- ;
950 
951  a1 = (*(*fMatSandiaMatrix)[i])[1];
952  a2 = (*(*fMatSandiaMatrix)[i])[2];
953  a3 = (*(*fMatSandiaMatrix)[i])[3];
954  a4 = (*(*fMatSandiaMatrix)[i])[4];
955 
956  lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
957 
958  return lambda ;
959 }
960 
962 //
963 //
964 
966 //
967 //
968 
970 {
971  G4double loss = 0.0 ;
972  loss *= step;
973 
974  return loss ;
975 }
976 
978 //
979 //
980 
982 {
983  G4double loss = 0.0 ;
984  loss *= step;
985 
986  return loss ;
987 }
988 
990 //
991 //
992 
994 {
995 
996 
997  G4double loss = 0.0 ;
998  loss *= step;
999  return loss ;
1000 }
1001 
1002 
1003 //
1004 // end of G4InitXscPAI implementation file
1005 //
1007