Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4InitXscPAI.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 62 //////////////////////////////////////////////////////////////////
 63 //
 64 // Constructor
 65 //
 66 
 67 using namespace std;
 68 
 69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC)
 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();
 80   fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
 81   matIndex         = (G4int)matCC->GetMaterial()->GetIndex();
 82 
 83   fSandia          = new G4SandiaTable(matIndex);
 84   fIntervalNumber  = fSandia->GetMaxInterval()-1;
 85 
 86   fMatSandiaMatrix = new G4OrderedTable();
 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   }
101   KillCloseIntervals();
102   Normalisation();
103   fBetaGammaSq = fTmax = 0.0;
104   fIntervalTmax = fCurrentInterval = 0;
105 }
106 
107 
108 
109 
110 ////////////////////////////////////////////////////////////////////////////
111 //
112 // Destructor
113 
114 G4InitXscPAI::~G4InitXscPAI()
115 {
116   delete fPAIxscVector;
117   delete fPAIdEdxVector;
118   delete fPAIphotonVector;
119   delete fPAIelectronVector;
120   delete fChCosSqVector;
121   delete fChWidthVector;
122   delete fSandia;
123   delete fMatSandiaMatrix;
124 }
125 
126 ////////////////////////////////////////////////////////////////////////
127 //
128 // Kill close intervals, recalculate fIntervalNumber
129 
130 void G4InitXscPAI::KillCloseIntervals()
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 
157 ////////////////////////////////////////////////////////////////////////
158 //
159 // Kill close intervals, recalculate fIntervalNumber
160 
161 void G4InitXscPAI::Normalisation()
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   }
180   fNormalizationCof  = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
181   fNormalizationCof *= fElectronDensity;
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 
226 ////////////////////////////////////////////////////////////////////
227 //
228 // Integration over electrons that could be considered
229 // quasi-free at energy transfer of interest
230 
231 G4double G4InitXscPAI::RutherfordIntegral( G4int k,
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 
251 ///////////////////////////////////////////////////////////////
252 //
253 //  Integrate photo-absorption cross-section from I1 up to omega
254 
255 G4double G4InitXscPAI::IntegralTerm(G4double omega)
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 
288 ////////////////////////////////////////////////////////////////
289 //
290 // Imaginary part of dielectric constant
291 // (G4int k - interval number, G4double en1 - energy point)
292 
293 G4double G4InitXscPAI::ImPartDielectricConst( G4int    k ,
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 
314 ////////////////////////////////////////////////////////////////
315 //
316 // Modulus squared of dielectric constant
317 // (G4int k - interval number, G4double omega - energy point)
318 
319 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int    k ,
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 
336 //////////////////////////////////////////////////////////////////////////////
337 //
338 // Real part of dielectric constant minus unit: epsilon_1 - 1
339 // (G4double enb - energy point)
340 //
341 
342 G4double G4InitXscPAI::RePartDielectricConst(G4double enb)
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 
407 //////////////////////////////////////////////////////////////////////
408 //
409 // PAI differential cross-section in terms of
410 // simplified Allison's equation
411 //
412 
413 G4double G4InitXscPAI::DifPAIxSection( G4double omega )
414 {
415   G4int i = fCurrentInterval;
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 
471 //////////////////////////////////////////////////////////////////////
472 //
473 // Differential PAI dEdx(omega)=omega*dNdx(omega)
474 //
475 
476 G4double G4InitXscPAI::DifPAIdEdx( G4double omega )
477 {
478   G4double dEdx = omega*DifPAIxSection(omega);
479   return dEdx;
480 }
481 
482 //////////////////////////////////////////////////////////////////////////
483 //
484 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
485 
486 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega  )
487 {        
488   G4int i = fCurrentInterval;
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 
545 //////////////////////////////////////////////////////////////////////////
546 //
547 // Calculation od dN/dx of collisions with creation of longitudinal EM
548 // excitations (plasmons, delta-electrons)
549 
550 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega )
551 {        
552   G4int i = fCurrentInterval;
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 
590 ////////////////////////////////////////////////////////////////////////
591 //
592 // Calculation of the PAI integral cross-section
593 // = specific primary ionisation, 1/cm
594 // 
595 
596 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax)
597 {
598   G4int i,k,i1,i2;
599   G4double energy1, energy2, result = 0.;
600  
601   fBetaGammaSq = bg2;
602   fTmax        = Tmax;
603 
604   delete fPAIxscVector;
605   
606   fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
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 
617   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
652                            &G4InitXscPAI::DifPAIxSection,
653                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
654 
655   else if( i == i1 ) result += integral.Legendre10(this,
656                            &G4InitXscPAI::DifPAIxSection,energy1,
657                            (*(*fMatSandiaMatrix)[i+1])[0]);
658 
659         else               result += integral.Legendre10(this,
660                            &G4InitXscPAI::DifPAIxSection,
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 
671 ////////////////////////////////////////////////////////////////////////
672 //
673 // Calculation of the PAI integral dEdx
674 // = mean energy loss per unit length, keV/cm
675 // 
676 
677 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax)
678 {
679   G4int i,k,i1,i2;
680   G4double energy1, energy2, result = 0.;
681  
682   fBetaGammaSq = bg2;
683   fTmax        = Tmax;
684 
685   delete fPAIdEdxVector;
686   
687   fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
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 
698   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
733                            &G4InitXscPAI::DifPAIdEdx,
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,
741                            &G4InitXscPAI::DifPAIdEdx,
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 
751 ////////////////////////////////////////////////////////////////////////
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 
757 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax)
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 
766   delete fPAIphotonVector;
767   delete fChCosSqVector;
768   delete fChWidthVector;
769   
770   fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
771   fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
772   fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
773 
774   fPAIphotonVector->PutValue(fPAIbin-1,result);
775   fChCosSqVector->PutValue(fPAIbin-1,1.);
776   fChWidthVector->PutValue(fPAIbin-1,1e-7);
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 
786   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
829                            &G4InitXscPAI::PAIdNdxCherenkov,
830                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
831 
832   else if( i == i1 ) result += integral.Legendre10(this,
833                            &G4InitXscPAI::PAIdNdxCherenkov,energy1,
834                            (*(*fMatSandiaMatrix)[i+1])[0]);
835 
836         else               result += integral.Legendre10(this,
837                            &G4InitXscPAI::PAIdNdxCherenkov,
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 
847 ////////////////////////////////////////////////////////////////////////
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 
853 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax)
854 {
855   G4int i,k,i1,i2;
856   G4double energy1, energy2, result = 0.;
857  
858   fBetaGammaSq = bg2;
859   fTmax        = Tmax;
860 
861   delete fPAIelectronVector;
862   
863   fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
864   fPAIelectronVector->PutValue(fPAIbin-1,result);
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 
874   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
909                            &G4InitXscPAI::PAIdNdxPlasmon,
910                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
911 
912   else if( i == i1 ) result += integral.Legendre10(this,
913                            &G4InitXscPAI::PAIdNdxPlasmon,energy1,
914                            (*(*fMatSandiaMatrix)[i+1])[0]);
915 
916         else               result += integral.Legendre10(this,
917                            &G4InitXscPAI::PAIdNdxPlasmon,
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 
928 /////////////////////////////////////////////////////////////////////////
929 //
930 //
931 
932 G4double G4InitXscPAI::GetPhotonLambda( G4double omega )
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 
961 /////////////////////////////////////////////////////////////////////////
962 //
963 //
964 
965 /////////////////////////////////////////////////////////////////////////
966 //
967 //
968 
969 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step )
970 {  
971   G4double loss = 0.0 ;
972   loss *= step;
973 
974   return loss ;
975 }
976 
977 /////////////////////////////////////////////////////////////////////////
978 //
979 //
980 
981 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step )
982 {  
983   G4double loss = 0.0 ;
984   loss *= step;
985 
986   return loss ;
987 }
988 
989 /////////////////////////////////////////////////////////////////////////
990 //
991 //
992 
993 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step )
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 //
1006 ////////////////////////////////////////////////////////////////////////////
1007 
1008