source: trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc @ 1316

Last change on this file since 1316 was 1228, checked in by garnier, 14 years ago

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