source: trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.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: 64.0 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: G4PAIxSection.cc,v 1.24 2008/05/30 16:04:40 grichine Exp $
28// GEANT4 tag $Name: geant4-09-03 $
29//
30//
31// G4PAIxSection.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// 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy
43// 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation
44// 17.05.01 V. Grichine, low energy extension down to 10*keV of proton
45// 20.11.98 adapted to a new Material/SandiaTable interface, mma
46// 11.06.97 V. Grichine, 1st version
47//
48
49
50
51#include "G4PAIxSection.hh"
52
53#include "globals.hh"
54#include "G4ios.hh"
55#include "G4Poisson.hh"
56#include "G4Material.hh"
57#include "G4MaterialCutsCouple.hh"
58#include "G4SandiaTable.hh"
59
60using namespace std;
61
62/* ******************************************************************
63
64// Init  array of Lorentz factors
65
66const G4double G4PAIxSection::fLorentzFactor[22] =
67{
68          0.0 ,     1.1 ,   1.2 ,   1.3 ,    1.5 ,    1.8 ,  2.0 ,
69          2.5 ,     3.0 ,   4.0 ,   7.0 ,   10.0 ,   20.0 , 40.0 ,
70         70.0 ,   100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 ,
71      10000.0 , 50000.0
72};
73
74const G4int G4PAIxSection::
75fRefGammaNumber = 29;         // The number of gamma for creation of
76                               // spline (9)
77
78***************************************************************** */ 
79
80// Local class constants
81
82const G4double G4PAIxSection::fDelta = 0.005; // energy shift from interval border
83const G4double G4PAIxSection::fError = 0.005; // error in lin-log approximation
84
85const G4int G4PAIxSection::fMaxSplineSize = 500;  // Max size of output spline
86                                                    // arrays
87
88//////////////////////////////////////////////////////////////////
89//
90// Constructor
91//
92
93G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
94{
95  fDensity       = matCC->GetMaterial()->GetDensity();
96  G4int matIndex = matCC->GetMaterial()->GetIndex();
97  fMaterialIndex = matIndex;   
98  fSandia        = new G4SandiaTable(matIndex);
99
100  G4int i, j; 
101  fMatSandiaMatrix = new G4OrderedTable();
102 
103  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
104  {
105     fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
106  }                     
107  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
108  {
109    (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
110
111    for(j = 1; j < 5; j++)
112    {
113      (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
114    }     
115  }                     
116}
117
118////////////////////////////////////////////////////////////////
119
120G4PAIxSection::G4PAIxSection(G4int materialIndex,
121                             G4double maxEnergyTransfer)
122{
123   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
124   G4int i, j;   
125
126      fMaterialIndex   = materialIndex;   
127      fDensity                = (*theMaterialTable)[materialIndex]->GetDensity();
128      fElectronDensity        = (*theMaterialTable)[materialIndex]->
129                             GetElectronDensity();
130      fIntervalNumber         = (*theMaterialTable)[materialIndex]->
131                             GetSandiaTable()->GetMatNbOfIntervals();
132      fIntervalNumber--;     
133      // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
134
135      fEnergyInterval = new G4double[fIntervalNumber+2];
136      fA1             = new G4double[fIntervalNumber+2];
137      fA2             = new G4double[fIntervalNumber+2];
138      fA3             = new G4double[fIntervalNumber+2];
139      fA4             = new G4double[fIntervalNumber+2];
140
141      for(i = 1; i <= fIntervalNumber; i++ )
142      {
143         if(((*theMaterialTable)[materialIndex]->
144    GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
145              i > fIntervalNumber               )
146         {
147            fEnergyInterval[i] = maxEnergyTransfer;
148            fIntervalNumber = i;
149            break;
150         }
151         fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
152                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
153         fA1[i]             = (*theMaterialTable)[materialIndex]->
154                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
155         fA2[i]             = (*theMaterialTable)[materialIndex]->
156                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
157         fA3[i]             = (*theMaterialTable)[materialIndex]->
158                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
159         fA4[i]             = (*theMaterialTable)[materialIndex]->
160                              GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
161         // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
162         //                               <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
163      }   
164      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
165      {
166         fIntervalNumber++;
167         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
168      }
169
170      // Now checking, if two borders are too close together
171
172      for(i=1;i<fIntervalNumber;i++)
173      {
174        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
175           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
176        {
177          continue;
178        }
179        else
180        {
181          for(j=i;j<fIntervalNumber;j++)
182          {
183            fEnergyInterval[j] = fEnergyInterval[j+1];
184                        fA1[j] = fA1[j+1];
185                        fA2[j] = fA2[j+1];
186                        fA3[j] = fA3[j+1];
187                        fA4[j] = fA4[j+1];
188          }
189          fIntervalNumber--;
190          i--;
191        }
192      }
193
194
195      /* *********************************
196
197      fSplineEnergy          = new G4double[fMaxSplineSize];   
198      fRePartDielectricConst = new G4double[fMaxSplineSize];   
199      fImPartDielectricConst = new G4double[fMaxSplineSize];   
200      fIntegralTerm          = new G4double[fMaxSplineSize];   
201      fDifPAIxSection        = new G4double[fMaxSplineSize];   
202      fIntegralPAIxSection   = new G4double[fMaxSplineSize];   
203     
204      for(i=0;i<fMaxSplineSize;i++)
205      {
206         fSplineEnergy[i]          = 0.0;   
207         fRePartDielectricConst[i] = 0.0;   
208         fImPartDielectricConst[i] = 0.0;   
209         fIntegralTerm[i]          = 0.0;   
210         fDifPAIxSection[i]        = 0.0;   
211         fIntegralPAIxSection[i]   = 0.0;   
212      }
213      **************************************************  */   
214
215      InitPAI();  // create arrays allocated above
216     
217      delete[] fEnergyInterval;
218      delete[] fA1;
219      delete[] fA2;
220      delete[] fA3;
221      delete[] fA4;   
222}
223
224////////////////////////////////////////////////////////////////////////
225//
226// Constructor with beta*gamma square value
227
228G4PAIxSection::G4PAIxSection( G4int materialIndex,
229                              G4double maxEnergyTransfer,
230                              G4double betaGammaSq,
231                              G4double** photoAbsCof, 
232                              G4int intNumber                   )
233{
234   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
235   G4int i, j; 
236 
237      fMaterialIndex   = materialIndex;     
238      fDensity                = (*theMaterialTable)[materialIndex]->GetDensity();
239      fElectronDensity        = (*theMaterialTable)[materialIndex]->
240                             GetElectronDensity();
241
242      fIntervalNumber         = intNumber;
243      fIntervalNumber--;
244      //   G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
245 
246      fEnergyInterval = new G4double[fIntervalNumber+2];
247      fA1             = new G4double[fIntervalNumber+2];
248      fA2             = new G4double[fIntervalNumber+2];
249      fA3             = new G4double[fIntervalNumber+2];
250      fA4             = new G4double[fIntervalNumber+2];
251
252      for( i = 1; i <= fIntervalNumber; i++ )
253      {
254         if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
255             i > fIntervalNumber )
256         {
257            fEnergyInterval[i] = maxEnergyTransfer;
258            fIntervalNumber = i;
259            break;
260         }
261         fEnergyInterval[i] = photoAbsCof[i-1][0];
262         fA1[i]             = photoAbsCof[i-1][1];
263         fA2[i]             = photoAbsCof[i-1][2];
264         fA3[i]             = photoAbsCof[i-1][3];
265         fA4[i]             = photoAbsCof[i-1][4];
266         // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
267         //        <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
268      }
269  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;   
270      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
271      {
272         fIntervalNumber++;
273         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
274      }
275      for(i=1;i<=fIntervalNumber;i++)
276      {
277        //  G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
278        //    <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
279      }
280      // Now checking, if two borders are too close together
281
282      for( i = 1; i < fIntervalNumber; i++ )
283      {
284        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
285           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
286        {
287          continue;
288        }
289        else
290        {
291          for(j=i;j<fIntervalNumber;j++)
292          {
293            fEnergyInterval[j] = fEnergyInterval[j+1];
294                        fA1[j] = fA1[j+1];
295                        fA2[j] = fA2[j+1];
296                        fA3[j] = fA3[j+1];
297                        fA4[j] = fA4[j+1];
298          }
299          fIntervalNumber--;
300          i--;
301        }
302      }
303
304      // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
305     
306      G4double   betaGammaSqRef = 
307      fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
308
309      NormShift(betaGammaSqRef);             
310      SplainPAI(betaGammaSqRef);
311     
312      // Preparation of integral PAI cross section for input betaGammaSq
313   
314      for(i = 1; i <= fSplineNumber; i++)
315      {
316         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
317         fdNdxMM[i]   = PAIdNdxMM(i,betaGammaSq);
318         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
319         fdNdxResonance[i]  = PAIdNdxResonance(i,betaGammaSq);
320         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
321
322         // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
323         //    <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
324      }
325      IntegralCerenkov();
326      IntegralMM();
327      IntegralPlasmon();
328      IntegralResonance();
329      IntegralPAIxSection();
330     
331      delete[] fEnergyInterval;
332      delete[] fA1;
333      delete[] fA2;
334      delete[] fA3;
335      delete[] fA4;   
336}
337
338////////////////////////////////////////////////////////////////////////
339//
340// Test Constructor with beta*gamma square value
341
342G4PAIxSection::G4PAIxSection( G4int materialIndex,
343                              G4double maxEnergyTransfer,
344                              G4double betaGammaSq          )
345{
346   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
347
348   G4int i, j, numberOfElements;   
349
350   fMaterialIndex   = materialIndex;   
351   fDensity         = (*theMaterialTable)[materialIndex]->GetDensity();
352   fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
353   numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
354
355   G4int* thisMaterialZ = new G4int[numberOfElements];
356   
357   for( i = 0; i < numberOfElements; i++ )
358   {
359         thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
360                                      GetElement(i)->GetZ();
361   }
362   // fSandia = new G4SandiaTable(materialIndex);
363   fSandia = (*theMaterialTable)[materialIndex]->
364     GetSandiaTable();
365   G4SandiaTable     thisMaterialSandiaTable(materialIndex);
366   fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
367                           (thisMaterialZ,numberOfElements);   
368   fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
369                           ( thisMaterialZ ,
370                      (*theMaterialTable)[materialIndex]->GetFractionVector() ,
371                             numberOfElements,fIntervalNumber);
372
373   fIntervalNumber--;
374
375      fEnergyInterval = new G4double[fIntervalNumber+2];
376      fA1             = new G4double[fIntervalNumber+2];
377      fA2             = new G4double[fIntervalNumber+2];
378      fA3             = new G4double[fIntervalNumber+2];
379      fA4             = new G4double[fIntervalNumber+2];
380
381      for( i = 1; i <= fIntervalNumber; i++ )
382      {
383  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
384          i > fIntervalNumber)
385         {
386            fEnergyInterval[i] = maxEnergyTransfer;
387            fIntervalNumber = i;
388            break;
389         }
390   fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
391   fA1[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
392   fA2[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
393   fA3[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
394   fA4[i]             = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
395
396      }   
397      if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
398      {
399         fIntervalNumber++;
400         fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
401         fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
402         fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
403         fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
404         fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
405      }
406      for(i=1;i<=fIntervalNumber;i++)
407      {
408        // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
409        //   <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
410      }
411      // Now checking, if two borders are too close together
412
413      for( i = 1; i < fIntervalNumber; i++ )
414      {
415        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
416           1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
417        {
418          continue;
419        }
420        else
421        {
422          for( j = i; j < fIntervalNumber; j++ )
423          {
424            fEnergyInterval[j] = fEnergyInterval[j+1];
425                        fA1[j] = fA1[j+1];
426                        fA2[j] = fA2[j+1];
427                        fA3[j] = fA3[j+1];
428                        fA4[j] = fA4[j+1];
429          }
430          fIntervalNumber--;
431          i--;
432        }
433      }
434
435      /* *********************************
436      fSplineEnergy          = new G4double[fMaxSplineSize];   
437      fRePartDielectricConst = new G4double[fMaxSplineSize];   
438      fImPartDielectricConst = new G4double[fMaxSplineSize];   
439      fIntegralTerm          = new G4double[fMaxSplineSize];   
440      fDifPAIxSection        = new G4double[fMaxSplineSize];   
441      fIntegralPAIxSection   = new G4double[fMaxSplineSize];   
442     
443      for(i=0;i<fMaxSplineSize;i++)
444      {
445         fSplineEnergy[i]          = 0.0;   
446         fRePartDielectricConst[i] = 0.0;   
447         fImPartDielectricConst[i] = 0.0;   
448         fIntegralTerm[i]          = 0.0;   
449         fDifPAIxSection[i]        = 0.0;   
450         fIntegralPAIxSection[i]   = 0.0;   
451      }
452      */ ////////////////////////
453
454      // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
455     
456      G4double   betaGammaSqRef = 
457      fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
458
459      NormShift(betaGammaSqRef);             
460      SplainPAI(betaGammaSqRef);
461     
462      // Preparation of integral PAI cross section for input betaGammaSq
463   
464      for(i = 1; i <= fSplineNumber; i++)
465      {
466         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
467         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
468         fdNdxMM[i]   = PAIdNdxMM(i,betaGammaSq);
469         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
470         fdNdxResonance[i]  = PAIdNdxResonance(i,betaGammaSq);
471      }
472      IntegralPAIxSection();
473      IntegralCerenkov();
474      IntegralMM();
475      IntegralPlasmon();
476      IntegralResonance();
477     
478      //   delete[] fEnergyInterval;
479      delete[] fA1;
480      delete[] fA2;
481      delete[] fA3;
482      delete[] fA4;   
483}
484
485
486////////////////////////////////////////////////////////////////////////////
487//
488// Destructor
489
490G4PAIxSection::~G4PAIxSection()
491{
492   /* ************************
493   delete[] fSplineEnergy         ;   
494   delete[] fRePartDielectricConst;   
495   delete[] fImPartDielectricConst;   
496   delete[] fIntegralTerm         ;   
497   delete[] fDifPAIxSection       ;   
498   delete[] fIntegralPAIxSection  ;
499   */ ////////////////////////
500}
501
502/////////////////////////////////////////////////////////////////////////
503//
504// General control function for class G4PAIxSection
505//
506
507void G4PAIxSection::InitPAI()
508{   
509   G4int i;
510   G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
511                          fLorentzFactor[fRefGammaNumber] - 1;
512
513   // Preparation of integral PAI cross section for reference gamma
514   
515   NormShift(betaGammaSq);             
516   SplainPAI(betaGammaSq);
517
518   IntegralPAIxSection();
519   IntegralCerenkov();
520   IntegralMM();
521   IntegralPlasmon();
522   IntegralResonance();
523
524   for(i = 0; i<= fSplineNumber; i++)
525   {
526      fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
527      if(i != 0) 
528      {
529         fPAItable[i][0] = fSplineEnergy[i];
530      }
531   }
532   fPAItable[0][0] = fSplineNumber;
533   
534   for(G4int j = 1; j < 112; j++)       // for other gammas
535   {
536      if( j == fRefGammaNumber ) continue;
537     
538      betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
539     
540      for(i = 1; i <= fSplineNumber; i++)
541      {
542         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
543         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
544         fdNdxMM[i]   = PAIdNdxMM(i,betaGammaSq);
545         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
546         fdNdxResonance[i]  = PAIdNdxResonance(i,betaGammaSq);
547      }
548      IntegralPAIxSection();
549      IntegralCerenkov();
550      IntegralMM();
551      IntegralPlasmon();
552      IntegralResonance();
553     
554      for(i = 0; i <= fSplineNumber; i++)
555      {
556         fPAItable[i][j] = fIntegralPAIxSection[i];
557      }
558   } 
559
560} 
561
562///////////////////////////////////////////////////////////////////////
563//
564// Shifting from borders to intervals Creation of first energy points
565//
566
567void G4PAIxSection::NormShift(G4double betaGammaSq)
568{
569  G4int i, j;
570
571  for( i = 1; i <= fIntervalNumber-1; i++ )
572  {
573    for( j = 1; j <= 2; j++ )
574    {
575      fSplineNumber = (i-1)*2 + j;
576
577      if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i  ]*(1+fDelta);
578      else         fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta); 
579      //    G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
580      //  <<fSplineEnergy[fSplineNumber]<<G4endl;
581    }
582  }
583  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
584
585  j = 1;
586
587  for( i = 2; i <= fSplineNumber; i++ )
588  {
589    if(fSplineEnergy[i]<fEnergyInterval[j+1])
590    {
591         fIntegralTerm[i] = fIntegralTerm[i-1] + 
592                            RutherfordIntegral(j,fSplineEnergy[i-1],
593                                                 fSplineEnergy[i]   );
594    }
595    else
596    {
597       G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
598                                           fEnergyInterval[j+1]   );
599         j++;
600         fIntegralTerm[i] = fIntegralTerm[i-1] + x + 
601                            RutherfordIntegral(j,fEnergyInterval[j],
602                                                 fSplineEnergy[i]    );
603    }
604    // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
605  } 
606  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
607  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
608
609  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
610
611          // Calculation of PAI differrential cross-section (1/(keV*cm))
612          // in the energy points near borders of energy intervals
613
614   for(G4int k = 1; k <= fIntervalNumber-1; k++ )
615   {
616      for( j = 1; j <= 2; j++ )
617      {
618         i = (k-1)*2 + j;
619         fImPartDielectricConst[i] = fNormalizationCof*
620                                     ImPartDielectricConst(k,fSplineEnergy[i]);
621         fRePartDielectricConst[i] = fNormalizationCof*
622                                     RePartDielectricConst(fSplineEnergy[i]);
623         fIntegralTerm[i] *= fNormalizationCof;
624
625         fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
626         fdNdxCerenkov[i]   = PAIdNdxCerenkov(i,betaGammaSq);
627         fdNdxMM[i]   = PAIdNdxMM(i,betaGammaSq);
628         fdNdxPlasmon[i]    = PAIdNdxPlasmon(i,betaGammaSq);
629         fdNdxResonance[i]    = PAIdNdxResonance(i,betaGammaSq);
630      }
631   }
632
633}  // end of NormShift
634
635/////////////////////////////////////////////////////////////////////////
636//
637// Creation of new energy points as geometrical mean of existing
638// one, calculation PAI_cs for them, while the error of logarithmic
639// linear approximation would be smaller than 'fError'
640
641void G4PAIxSection::SplainPAI(G4double betaGammaSq)
642{
643   G4int k = 1;
644   G4int i = 1;
645
646   while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
647   {
648      if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
649      {
650          k++;   // Here next energy point is in next energy interval
651          i++;
652          continue;
653      }
654                       // Shifting of arrayes for inserting the geometrical
655                       // average of 'i' and 'i+1' energy points to 'i+1' place
656      fSplineNumber++;
657
658      for(G4int j = fSplineNumber; j >= i+2; j-- )
659      {
660         fSplineEnergy[j]          = fSplineEnergy[j-1];
661         fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
662         fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
663         fIntegralTerm[j]          = fIntegralTerm[j-1];
664
665         fDifPAIxSection[j] = fDifPAIxSection[j-1];
666         fdNdxCerenkov[j]   = fdNdxCerenkov[j-1];
667         fdNdxMM[j]   = fdNdxMM[j-1];
668         fdNdxPlasmon[j]    = fdNdxPlasmon[j-1];
669         fdNdxResonance[j]  = fdNdxResonance[j-1];
670      }
671      G4double x1  = fSplineEnergy[i];
672      G4double x2  = fSplineEnergy[i+1];
673      G4double yy1 = fDifPAIxSection[i];
674      G4double y2  = fDifPAIxSection[i+1];
675
676      G4double en1 = sqrt(x1*x2);
677      fSplineEnergy[i+1] = en1;
678
679                 // Calculation of logarithmic linear approximation
680                 // in this (enr) energy point, which number is 'i+1' now
681
682      G4double a = log10(y2/yy1)/log10(x2/x1);
683      G4double b = log10(yy1) - a*log10(x1);
684      G4double y = a*log10(en1) + b;
685      y = pow(10.,y);
686
687                 // Calculation of the PAI dif. cross-section at this point
688
689      fImPartDielectricConst[i+1] = fNormalizationCof*
690                                    ImPartDielectricConst(k,fSplineEnergy[i+1]);
691      fRePartDielectricConst[i+1] = fNormalizationCof*
692                                    RePartDielectricConst(fSplineEnergy[i+1]);
693      fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
694                           RutherfordIntegral(k,fSplineEnergy[i],
695                                                fSplineEnergy[i+1]);
696
697      fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
698      fdNdxCerenkov[i+1]   = PAIdNdxCerenkov(i+1,betaGammaSq);
699      fdNdxMM[i+1]   = PAIdNdxMM(i+1,betaGammaSq);
700      fdNdxPlasmon[i+1]    = PAIdNdxPlasmon(i+1,betaGammaSq);
701      fdNdxResonance[i+1]  = PAIdNdxResonance(i+1,betaGammaSq);
702
703                  // Condition for next division of this segment or to pass
704                  // to higher energies
705
706      G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
707
708      if( x < 0 ) 
709      {
710         x = -x;
711      }
712      if( x > fError && fSplineNumber < fMaxSplineSize-1 )
713      {
714         continue;  // next division
715      }
716      i += 2;  // pass to next segment
717
718   }   // close 'while'
719
720}  // end of SplainPAI
721
722
723////////////////////////////////////////////////////////////////////
724//
725// Integration over electrons that could be considered
726// quasi-free at energy transfer of interest
727
728G4double G4PAIxSection::RutherfordIntegral( G4int k,
729                                            G4double x1,
730                                            G4double x2   )
731{
732   G4double  c1, c2, c3;
733   // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;   
734   c1 = (x2 - x1)/x1/x2;
735   c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
736   c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
737   // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
738   
739   return  fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
740
741}   // end of RutherfordIntegral
742
743
744/////////////////////////////////////////////////////////////////
745//
746// Imaginary part of dielectric constant
747// (G4int k - interval number, G4double en1 - energy point)
748
749G4double G4PAIxSection::ImPartDielectricConst( G4int    k ,
750                                               G4double energy1 )
751{
752   G4double energy2,energy3,energy4,result;
753
754   energy2 = energy1*energy1;
755   energy3 = energy2*energy1;
756   energy4 = energy3*energy1;
757   
758   result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4; 
759   result *=hbarc/energy1;
760   
761   return result;
762
763}  // end of ImPartDielectricConst
764
765/////////////////////////////////////////////////////////////////
766//
767// Returns lambda of photon with energy1 in current material
768
769G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
770{
771  G4int i;
772  G4double energy2, energy3, energy4, result, lambda;
773
774  energy2 = energy1*energy1;
775  energy3 = energy2*energy1;
776  energy4 = energy3*energy1;
777
778  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
779  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
780  // result *= fDensity;
781
782  for( i = 1; i <= fIntervalNumber; i++ )
783  {
784     if( energy1 < fEnergyInterval[i]) break;
785  }
786  i--;
787  if(i == 0) i = 1;
788
789  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4; 
790
791  if( result > DBL_MIN ) lambda = 1./result;
792  else                   lambda = DBL_MAX;
793   
794  return lambda;
795} 
796
797/////////////////////////////////////////////////////////////////
798//
799// Return lambda of electron with energy1 in current material
800// parametrisation from NIM A554(2005)474-493
801
802G4double G4PAIxSection::GetElectronRange( G4double energy )
803{
804  G4double range;
805  /*
806  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
807
808  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
809  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
810
811  energy /= keV; // energy in keV in parametrised formula
812
813  if (energy < 10.)
814  {
815    range = 3.872e-3*A/Z;
816    range *= pow(energy,1.492);
817  }
818  else
819  {
820    range = 6.97e-3*pow(energy,1.6);
821  }
822  */
823  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
824
825  G4double cofA = 5.37e-4*g/cm2/keV;
826  G4double cofB = 0.9815;
827  G4double cofC = 3.123e-3/keV;
828  // energy /= keV;
829
830  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) ); 
831
832  // range *= g/cm2;
833  range /= fDensity;
834
835  return range;
836}
837
838//////////////////////////////////////////////////////////////////////////////
839//
840// Real part of dielectric constant minus unit: epsilon_1 - 1
841// (G4double enb - energy point)
842//
843
844G4double G4PAIxSection::RePartDielectricConst(G4double enb)
845{       
846   G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
847            c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
848
849   x0 = enb;
850   result = 0;
851   
852   for(G4int i=1;i<=fIntervalNumber-1;i++)
853   {
854      x1 = fEnergyInterval[i];
855      x2 = fEnergyInterval[i+1];
856      xx1 = x1 - x0;
857      xx2 = x2 - x0;
858      xx12 = xx2/xx1;
859     
860      if(xx12<0)
861      {
862         xx12 = -xx12;
863      }
864      xln1 = log(x2/x1);
865      xln2 = log(xx12);
866      xln3 = log((x2 + x0)/(x1 + x0));
867      x02 = x0*x0;
868      x03 = x02*x0;
869      x04 = x03*x0;
870      x05 = x04*x0;
871      c1  = (x2 - x1)/x1/x2;
872      c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
873      c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
874
875      result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
876      result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
877      result -= fA3[i]*c2/2/x02;
878      result -= fA4[i]*c3/3/x02;
879
880      cof1 = fA1[i]/x02 + fA3[i]/x04;
881      cof2 = fA2[i]/x03 + fA4[i]/x05;
882
883      result += 0.5*(cof1 +cof2)*xln2;
884      result += 0.5*(cof1 - cof2)*xln3;
885   } 
886   result *= 2*hbarc/pi;
887   
888   return result;
889
890}   // end of RePartDielectricConst
891
892//////////////////////////////////////////////////////////////////////
893//
894// PAI differential cross-section in terms of
895// simplified Allison's equation
896//
897
898G4double G4PAIxSection::DifPAIxSection( G4int              i ,
899                                        G4double betaGammaSq  )
900{       
901   G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
902   //G4double beta, be4;
903   G4double be4;
904   G4double betaBohr2 = fine_structure_const*fine_structure_const;
905   G4double betaBohr4 = betaBohr2*betaBohr2*4.0;
906   be2 = betaGammaSq/(1 + betaGammaSq);
907   be4 = be2*be2;
908   //  beta = sqrt(be2);
909   cof = 1;
910   x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
911
912   if( betaGammaSq < 0.01 ) x2 = log(be2);
913   else
914   {
915     x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
916                (1/betaGammaSq - fRePartDielectricConst[i]) + 
917                fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
918   }
919   if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
920   {
921     x6=0;
922   }
923   else
924   {
925     x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
926     x5 = -1 - fRePartDielectricConst[i] +
927          be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
928          fImPartDielectricConst[i]*fImPartDielectricConst[i]);
929
930     x7 = atan2(fImPartDielectricConst[i],x3);
931     x6 = x5 * x7;
932   }
933    // if(fImPartDielectricConst[i] == 0) x6 = 0;
934   
935   x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
936   //   if( x4 < 0.0 ) x4 = 0.0;
937   x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
938        fImPartDielectricConst[i]*fImPartDielectricConst[i];
939
940   result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
941   if(result < 1.0e-8) result = 1.0e-8;
942   result *= fine_structure_const/be2/pi;
943   //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr));
944   //  result *= (1-exp(-be2/betaBohr2));
945   result *= (1-exp(-be4/betaBohr4));
946   if(fDensity >= 0.1)
947   { 
948      result /= x8;
949   }
950   return result;
951
952} // end of DifPAIxSection
953
954//////////////////////////////////////////////////////////////////////////
955//
956// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
957
958G4double G4PAIxSection::PAIdNdxCerenkov( G4int    i ,
959                                         G4double betaGammaSq  )
960{       
961   G4double logarithm, x3, x5, argument, modul2, dNdxC; 
962   G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
963
964   cofBetaBohr = 4.0;
965   betaBohr2   = fine_structure_const*fine_structure_const;
966   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr;
967
968   be2 = betaGammaSq/(1 + betaGammaSq);
969   be4 = be2*be2;
970
971   if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
972   else
973   {
974     logarithm  = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
975                        (1/betaGammaSq - fRePartDielectricConst[i]) + 
976                        fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
977     logarithm += log(1+1.0/betaGammaSq);
978   }
979
980   if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
981   {
982     argument = 0.0;
983   }
984   else
985   {
986     x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
987     x5 = -1.0 - fRePartDielectricConst[i] +
988          be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
989          fImPartDielectricConst[i]*fImPartDielectricConst[i]);
990     if( x3 == 0.0 ) argument = 0.5*pi;
991     else            argument = atan2(fImPartDielectricConst[i],x3);
992     argument *= x5 ;
993   }   
994   dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
995 
996   if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
997
998   dNdxC *= fine_structure_const/be2/pi;
999
1000   dNdxC *= (1-exp(-be4/betaBohr4));
1001
1002   if(fDensity >= 0.1)
1003   { 
1004      modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 
1005                    fImPartDielectricConst[i]*fImPartDielectricConst[i];
1006      dNdxC /= modul2;
1007   }
1008   return dNdxC;
1009
1010} // end of PAIdNdxCerenkov
1011
1012//////////////////////////////////////////////////////////////////////////
1013//
1014// Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons
1015
1016G4double G4PAIxSection::PAIdNdxMM( G4int    i ,
1017                                         G4double betaGammaSq  )
1018{       
1019   G4double logarithm, x3, x5, argument, dNdxC; 
1020   G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
1021
1022   cofBetaBohr = 4.0;
1023   betaBohr2   = fine_structure_const*fine_structure_const;
1024   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr;
1025
1026   be2 = betaGammaSq/(1 + betaGammaSq);
1027   be4 = be2*be2;
1028
1029   if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1030   else
1031   {
1032     logarithm  = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1033                        (1/betaGammaSq - fRePartDielectricConst[i]) + 
1034                        fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1035     logarithm += log(1+1.0/betaGammaSq);
1036   }
1037
1038   if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1039   {
1040     argument = 0.0;
1041   }
1042   else
1043   {
1044     x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1045     x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1046     if( x3 == 0.0 ) argument = 0.5*pi;
1047     else            argument = atan2(fImPartDielectricConst[i],x3);
1048     argument *= x5 ;
1049   }   
1050   dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1051 
1052   if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1053
1054   dNdxC *= fine_structure_const/be2/pi;
1055
1056   dNdxC *= (1-exp(-be4/betaBohr4));
1057   return dNdxC;
1058
1059} // end of PAIdNdxMM
1060
1061//////////////////////////////////////////////////////////////////////////
1062//
1063// Calculation od dN/dx of collisions with creation of longitudinal EM
1064// excitations (plasmons, delta-electrons)
1065
1066G4double G4PAIxSection::PAIdNdxPlasmon( G4int    i ,
1067                                        G4double betaGammaSq  )
1068{       
1069   G4double resonance, modul2, dNdxP, cof = 1.;
1070   G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1071
1072 
1073   cofBetaBohr = 4.0;
1074   betaBohr2   = fine_structure_const*fine_structure_const;
1075   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr;
1076
1077   be2 = betaGammaSq/(1 + betaGammaSq);
1078   be4 = be2*be2;
1079 
1080   resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]); 
1081   resonance *= fImPartDielectricConst[i]/hbarc;
1082
1083
1084   dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1085
1086   if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1087
1088   dNdxP *= fine_structure_const/be2/pi;
1089   dNdxP *= (1-exp(-be4/betaBohr4));
1090
1091   if( fDensity >= 0.1 )
1092   { 
1093     modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
1094        fImPartDielectricConst[i]*fImPartDielectricConst[i];
1095     dNdxP /= modul2;
1096   }
1097   return dNdxP;
1098
1099} // end of PAIdNdxPlasmon
1100
1101//////////////////////////////////////////////////////////////////////////
1102//
1103// Calculation od dN/dx of collisions with creation of longitudinal EM
1104// resonance excitations (plasmons, delta-electrons)
1105
1106G4double G4PAIxSection::PAIdNdxResonance( G4int    i ,
1107                                        G4double betaGammaSq  )
1108{       
1109   G4double resonance, modul2, dNdxP;
1110   G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1111
1112   cofBetaBohr = 4.0;
1113   betaBohr2   = fine_structure_const*fine_structure_const;
1114   betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr;
1115
1116   be2 = betaGammaSq/(1 + betaGammaSq);
1117   be4 = be2*be2;
1118 
1119   resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]); 
1120   resonance *= fImPartDielectricConst[i]/hbarc;
1121
1122
1123   dNdxP = resonance;
1124
1125   if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1126
1127   dNdxP *= fine_structure_const/be2/pi;
1128   dNdxP *= (1-exp(-be4/betaBohr4));
1129
1130   if( fDensity >= 0.1 )
1131   { 
1132     modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
1133        fImPartDielectricConst[i]*fImPartDielectricConst[i];
1134     dNdxP /= modul2;
1135   }
1136   return dNdxP;
1137
1138} // end of PAIdNdxResonance
1139
1140////////////////////////////////////////////////////////////////////////
1141//
1142// Calculation of the PAI integral cross-section
1143// fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
1144// and fIntegralPAIxSection[0] = mean energy loss per cm  in keV/cm
1145
1146void G4PAIxSection::IntegralPAIxSection()
1147{
1148  fIntegralPAIxSection[fSplineNumber] = 0;
1149  fIntegralPAIdEdx[fSplineNumber]     = 0;
1150  fIntegralPAIxSection[0]             = 0;
1151  G4int k = fIntervalNumber -1;
1152
1153  for(G4int i = fSplineNumber-1; i >= 1; i--)
1154  {
1155    if(fSplineEnergy[i] >= fEnergyInterval[k])
1156    {
1157      fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1158      fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1159    }
1160    else
1161    {
1162      fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + 
1163                                   SumOverBorder(i+1,fEnergyInterval[k]);
1164      fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + 
1165                                   SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1166      k--;
1167    }
1168  }
1169}   // end of IntegralPAIxSection
1170
1171////////////////////////////////////////////////////////////////////////
1172//
1173// Calculation of the PAI Cerenkov integral cross-section
1174// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
1175// and fIntegralCerenkov[0] = mean Cerenkov loss per cm  in keV/cm
1176
1177void G4PAIxSection::IntegralCerenkov()
1178{
1179  G4int i, k;
1180   fIntegralCerenkov[fSplineNumber] = 0;
1181   fIntegralCerenkov[0] = 0;
1182   k = fIntervalNumber -1;
1183
1184   for( i = fSplineNumber-1; i >= 1; i-- )
1185   {
1186      if(fSplineEnergy[i] >= fEnergyInterval[k])
1187      {
1188        fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1189        // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1190      }
1191      else
1192      {
1193        fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + 
1194                                   SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1195        k--;
1196        // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1197      }
1198   }
1199
1200}   // end of IntegralCerenkov
1201
1202////////////////////////////////////////////////////////////////////////
1203//
1204// Calculation of the PAI MM-Cerenkov integral cross-section
1205// fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm
1206// and fIntegralMM[0] = mean MM-Cerenkov loss per cm  in keV/cm
1207
1208void G4PAIxSection::IntegralMM()
1209{
1210  G4int i, k;
1211   fIntegralMM[fSplineNumber] = 0;
1212   fIntegralMM[0] = 0;
1213   k = fIntervalNumber -1;
1214
1215   for( i = fSplineNumber-1; i >= 1; i-- )
1216   {
1217      if(fSplineEnergy[i] >= fEnergyInterval[k])
1218      {
1219        fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1220        // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1221      }
1222      else
1223      {
1224        fIntegralMM[i] = fIntegralMM[i+1] + 
1225                                   SumOverBordMM(i+1,fEnergyInterval[k]);
1226        k--;
1227        // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1228      }
1229   }
1230
1231}   // end of IntegralMM
1232
1233////////////////////////////////////////////////////////////////////////
1234//
1235// Calculation of the PAI Plasmon integral cross-section
1236// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1237// and fIntegralPlasmon[0] = mean plasmon loss per cm  in keV/cm
1238
1239void G4PAIxSection::IntegralPlasmon()
1240{
1241   fIntegralPlasmon[fSplineNumber] = 0;
1242   fIntegralPlasmon[0] = 0;
1243   G4int k = fIntervalNumber -1;
1244   for(G4int i=fSplineNumber-1;i>=1;i--)
1245   {
1246      if(fSplineEnergy[i] >= fEnergyInterval[k])
1247      {
1248        fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1249      }
1250      else
1251      {
1252        fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + 
1253                                   SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1254        k--;
1255      }
1256   }
1257
1258}   // end of IntegralPlasmon
1259
1260////////////////////////////////////////////////////////////////////////
1261//
1262// Calculation of the PAI resonance integral cross-section
1263// fIntegralResonance[1] = resonance primary ionisation, 1/cm
1264// and fIntegralResonance[0] = mean resonance loss per cm  in keV/cm
1265
1266void G4PAIxSection::IntegralResonance()
1267{
1268   fIntegralResonance[fSplineNumber] = 0;
1269   fIntegralResonance[0] = 0;
1270   G4int k = fIntervalNumber -1;
1271   for(G4int i=fSplineNumber-1;i>=1;i--)
1272   {
1273      if(fSplineEnergy[i] >= fEnergyInterval[k])
1274      {
1275        fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1276      }
1277      else
1278      {
1279        fIntegralResonance[i] = fIntegralResonance[i+1] + 
1280                                   SumOverBordResonance(i+1,fEnergyInterval[k]);
1281        k--;
1282      }
1283   }
1284
1285}   // end of IntegralResonance
1286
1287//////////////////////////////////////////////////////////////////////
1288//
1289// Calculation the PAI integral cross-section inside
1290// of interval of continuous values of photo-ionisation
1291// cross-section. Parameter  'i' is the number of interval.
1292
1293G4double G4PAIxSection::SumOverInterval( G4int i )
1294{         
1295   G4double x0,x1,y0,yy1,a,b,c,result;
1296
1297   x0 = fSplineEnergy[i];
1298   x1 = fSplineEnergy[i+1];
1299   y0 = fDifPAIxSection[i];
1300   yy1 = fDifPAIxSection[i+1];
1301   c = x1/x0;
1302   a = log10(yy1/y0)/log10(c);
1303   // b = log10(y0) - a*log10(x0);
1304   b = y0/pow(x0,a);
1305   a += 1;
1306   if(a == 0) 
1307   {
1308      result = b*log(x1/x0);
1309   }
1310   else
1311   {
1312      result = y0*(x1*pow(c,a-1) - x0)/a;
1313   }
1314   a++;
1315   if(a == 0) 
1316   {
1317      fIntegralPAIxSection[0] += b*log(x1/x0);
1318   }
1319   else
1320   {
1321      fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1322   }
1323   return result;
1324
1325} //  end of SumOverInterval
1326
1327/////////////////////////////////
1328
1329G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1330{         
1331   G4double x0,x1,y0,yy1,a,b,c,result;
1332
1333   x0 = fSplineEnergy[i];
1334   x1 = fSplineEnergy[i+1];
1335   y0 = fDifPAIxSection[i];
1336   yy1 = fDifPAIxSection[i+1];
1337   c = x1/x0;
1338   a = log10(yy1/y0)/log10(c);
1339   // b = log10(y0) - a*log10(x0);
1340   b = y0/pow(x0,a);
1341   a += 2;
1342   if(a == 0) 
1343   {
1344     result = b*log(x1/x0);
1345   }
1346   else
1347   {
1348     result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1349   }
1350   return result;
1351
1352} //  end of SumOverInterval
1353
1354//////////////////////////////////////////////////////////////////////
1355//
1356// Calculation the PAI Cerenkov integral cross-section inside
1357// of interval of continuous values of photo-ionisation Cerenkov
1358// cross-section. Parameter  'i' is the number of interval.
1359
1360G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1361{         
1362   G4double x0,x1,y0,yy1,a,b,c,result;
1363
1364   x0  = fSplineEnergy[i];
1365   x1  = fSplineEnergy[i+1];
1366   y0  = fdNdxCerenkov[i];
1367   yy1 = fdNdxCerenkov[i+1];
1368   // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1369   //   <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1370
1371   c = x1/x0;
1372   a = log10(yy1/y0)/log10(c);
1373   b = y0/pow(x0,a);
1374
1375   a += 1.0;
1376   if(a == 0) result = b*log(c);
1377   else       result = y0*(x1*pow(c,a-1) - x0)/a;   
1378   a += 1.0;
1379
1380   if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1381   else         fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1382   //  G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;   
1383   return result;
1384
1385} //  end of SumOverInterCerenkov
1386
1387//////////////////////////////////////////////////////////////////////
1388//
1389// Calculation the PAI MM-Cerenkov integral cross-section inside
1390// of interval of continuous values of photo-ionisation Cerenkov
1391// cross-section. Parameter  'i' is the number of interval.
1392
1393G4double G4PAIxSection::SumOverInterMM( G4int i )
1394{         
1395   G4double x0,x1,y0,yy1,a,b,c,result;
1396
1397   x0  = fSplineEnergy[i];
1398   x1  = fSplineEnergy[i+1];
1399   y0  = fdNdxMM[i];
1400   yy1 = fdNdxMM[i+1];
1401   // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1402   //   <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1403
1404   c = x1/x0;
1405   a = log10(yy1/y0)/log10(c);
1406   b = y0/pow(x0,a);
1407
1408   a += 1.0;
1409   if(a == 0) result = b*log(c);
1410   else       result = y0*(x1*pow(c,a-1) - x0)/a;   
1411   a += 1.0;
1412
1413   if( a == 0 ) fIntegralMM[0] += b*log(x1/x0);
1414   else         fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1415   //  G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;   
1416   return result;
1417
1418} //  end of SumOverInterMM
1419
1420//////////////////////////////////////////////////////////////////////
1421//
1422// Calculation the PAI Plasmon integral cross-section inside
1423// of interval of continuous values of photo-ionisation Plasmon
1424// cross-section. Parameter  'i' is the number of interval.
1425
1426G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1427{         
1428   G4double x0,x1,y0,yy1,a,b,c,result;
1429
1430   x0  = fSplineEnergy[i];
1431   x1  = fSplineEnergy[i+1];
1432   y0  = fdNdxPlasmon[i];
1433   yy1 = fdNdxPlasmon[i+1];
1434   c =x1/x0;
1435   a = log10(yy1/y0)/log10(c);
1436   // b = log10(y0) - a*log10(x0);
1437   b = y0/pow(x0,a);
1438
1439   a += 1.0;
1440   if(a == 0) result = b*log(x1/x0);
1441   else       result = y0*(x1*pow(c,a-1) - x0)/a;   
1442   a += 1.0;
1443
1444   if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1445   else         fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1446   
1447   return result;
1448
1449} //  end of SumOverInterPlasmon
1450
1451//////////////////////////////////////////////////////////////////////
1452//
1453// Calculation the PAI resonance integral cross-section inside
1454// of interval of continuous values of photo-ionisation resonance
1455// cross-section. Parameter  'i' is the number of interval.
1456
1457G4double G4PAIxSection::SumOverInterResonance( G4int i )
1458{         
1459   G4double x0,x1,y0,yy1,a,b,c,result;
1460
1461   x0  = fSplineEnergy[i];
1462   x1  = fSplineEnergy[i+1];
1463   y0  = fdNdxResonance[i];
1464   yy1 = fdNdxResonance[i+1];
1465   c =x1/x0;
1466   a = log10(yy1/y0)/log10(c);
1467   // b = log10(y0) - a*log10(x0);
1468   b = y0/pow(x0,a);
1469
1470   a += 1.0;
1471   if(a == 0) result = b*log(x1/x0);
1472   else       result = y0*(x1*pow(c,a-1) - x0)/a;   
1473   a += 1.0;
1474
1475   if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1476   else         fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1477   
1478   return result;
1479
1480} //  end of SumOverInterResonance
1481
1482///////////////////////////////////////////////////////////////////////////////
1483//
1484// Integration of PAI cross-section for the case of
1485// passing across border between intervals
1486
1487G4double G4PAIxSection::SumOverBorder( G4int      i , 
1488                                       G4double en0    )
1489{               
1490   G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1491
1492   e0 = en0;
1493   x0 = fSplineEnergy[i];
1494   x1 = fSplineEnergy[i+1];
1495   y0 = fDifPAIxSection[i];
1496   yy1 = fDifPAIxSection[i+1];
1497
1498   c = x1/x0;
1499   d = e0/x0;   
1500   a = log10(yy1/y0)/log10(x1/x0);
1501   // b0 = log10(y0) - a*log10(x0);
1502   b = y0/pow(x0,a);  // pow(10.,b);
1503   
1504   a += 1;
1505   if(a == 0)
1506   {
1507      result = b*log(x0/e0);
1508   }
1509   else
1510   {
1511      result = y0*(x0 - e0*pow(d,a-1))/a;
1512   }
1513   a++;
1514   if(a == 0)
1515   {
1516      fIntegralPAIxSection[0] += b*log(x0/e0);
1517   }
1518   else 
1519   {
1520      fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1521   }
1522   x0 = fSplineEnergy[i - 1];
1523   x1 = fSplineEnergy[i - 2];
1524   y0 = fDifPAIxSection[i - 1];
1525   yy1 = fDifPAIxSection[i - 2];
1526
1527   c = x1/x0;
1528   d = e0/x0;   
1529   a = log10(yy1/y0)/log10(x1/x0);
1530   //  b0 = log10(y0) - a*log10(x0);
1531   b = y0/pow(x0,a);
1532   a += 1;
1533   if(a == 0)
1534   {
1535      result += b*log(e0/x0);
1536   }
1537   else
1538   {
1539      result += y0*(e0*pow(d,a-1) - x0)/a;
1540   }
1541   a++;
1542   if(a == 0) 
1543   {
1544      fIntegralPAIxSection[0] += b*log(e0/x0);
1545   }
1546   else
1547   {
1548      fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1549   }
1550   return result;
1551
1552} 
1553
1554///////////////////////////////////////////////////////////////////////
1555
1556G4double G4PAIxSection::SumOverBorderdEdx( G4int      i , 
1557                                       G4double en0    )
1558{               
1559   G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1560
1561   e0 = en0;
1562   x0 = fSplineEnergy[i];
1563   x1 = fSplineEnergy[i+1];
1564   y0 = fDifPAIxSection[i];
1565   yy1 = fDifPAIxSection[i+1];
1566
1567   c = x1/x0;
1568   d = e0/x0;   
1569   a = log10(yy1/y0)/log10(x1/x0);
1570   // b0 = log10(y0) - a*log10(x0);
1571   b = y0/pow(x0,a);  // pow(10.,b);
1572   
1573   a += 2;
1574   if(a == 0)
1575   {
1576      result = b*log(x0/e0);
1577   }
1578   else 
1579   {
1580      result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1581   }
1582   x0 = fSplineEnergy[i - 1];
1583   x1 = fSplineEnergy[i - 2];
1584   y0 = fDifPAIxSection[i - 1];
1585   yy1 = fDifPAIxSection[i - 2];
1586
1587   c = x1/x0;
1588   d = e0/x0;   
1589   a = log10(yy1/y0)/log10(x1/x0);
1590   //  b0 = log10(y0) - a*log10(x0);
1591   b = y0/pow(x0,a);
1592   a += 2;
1593   if(a == 0) 
1594   {
1595      result += b*log(e0/x0);
1596   }
1597   else
1598   {
1599      result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1600   }
1601   return result;
1602
1603} 
1604
1605///////////////////////////////////////////////////////////////////////////////
1606//
1607// Integration of Cerenkov cross-section for the case of
1608// passing across border between intervals
1609
1610G4double G4PAIxSection::SumOverBordCerenkov( G4int      i , 
1611                                             G4double en0    )
1612{               
1613   G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1614
1615   e0 = en0;
1616   x0 = fSplineEnergy[i];
1617   x1 = fSplineEnergy[i+1];
1618   y0 = fdNdxCerenkov[i];
1619   yy1 = fdNdxCerenkov[i+1];
1620
1621   //  G4cout<<G4endl;
1622   //  G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1623   //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1624   c = x1/x0;
1625   d = e0/x0;
1626   a = log10(yy1/y0)/log10(c);
1627   // b0 = log10(y0) - a*log10(x0);
1628   b = y0/pow(x0,a); // pow(10.,b0);   
1629   
1630   a += 1.0;
1631   if( a == 0 ) result = b*log(x0/e0);
1632   else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1633   a += 1.0;
1634
1635   if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
1636   else         fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1637
1638// G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1639   
1640   x0  = fSplineEnergy[i - 1];
1641   x1  = fSplineEnergy[i - 2];
1642   y0  = fdNdxCerenkov[i - 1];
1643   yy1 = fdNdxCerenkov[i - 2];
1644
1645   // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1646   //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1647
1648   c = x1/x0;
1649   d = e0/x0;
1650   a  = log10(yy1/y0)/log10(x1/x0);
1651   // b0 = log10(y0) - a*log10(x0);
1652   b  =  y0/pow(x0,a);  // pow(10.,b0);
1653
1654   a += 1.0;
1655   if( a == 0 ) result += b*log(e0/x0);
1656   else         result += y0*(e0*pow(d,a-1) - x0 )/a;
1657   a += 1.0;
1658
1659   if( a == 0 )   fIntegralCerenkov[0] += b*log(e0/x0);
1660   else           fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1661
1662   // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1663   // <<b<<"; result = "<<result<<G4endl;   
1664
1665   return result;
1666
1667} 
1668
1669///////////////////////////////////////////////////////////////////////////////
1670//
1671// Integration of MM-Cerenkov cross-section for the case of
1672// passing across border between intervals
1673
1674G4double G4PAIxSection::SumOverBordMM( G4int      i , 
1675                                             G4double en0    )
1676{               
1677   G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1678
1679   e0 = en0;
1680   x0 = fSplineEnergy[i];
1681   x1 = fSplineEnergy[i+1];
1682   y0 = fdNdxMM[i];
1683   yy1 = fdNdxMM[i+1];
1684
1685   //  G4cout<<G4endl;
1686   //  G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1687   //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1688   c = x1/x0;
1689   d = e0/x0;
1690   a = log10(yy1/y0)/log10(c);
1691   // b0 = log10(y0) - a*log10(x0);
1692   b = y0/pow(x0,a); // pow(10.,b0);   
1693   
1694   a += 1.0;
1695   if( a == 0 ) result = b*log(x0/e0);
1696   else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1697   a += 1.0;
1698
1699   if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
1700   else         fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1701
1702// G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1703   
1704   x0  = fSplineEnergy[i - 1];
1705   x1  = fSplineEnergy[i - 2];
1706   y0  = fdNdxMM[i - 1];
1707   yy1 = fdNdxMM[i - 2];
1708
1709   // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1710   //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1711
1712   c = x1/x0;
1713   d = e0/x0;
1714   a  = log10(yy1/y0)/log10(x1/x0);
1715   // b0 = log10(y0) - a*log10(x0);
1716   b  =  y0/pow(x0,a);  // pow(10.,b0);
1717
1718   a += 1.0;
1719   if( a == 0 ) result += b*log(e0/x0);
1720   else         result += y0*(e0*pow(d,a-1) - x0 )/a;
1721   a += 1.0;
1722
1723   if( a == 0 )   fIntegralMM[0] += b*log(e0/x0);
1724   else           fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1725
1726   // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1727   // <<b<<"; result = "<<result<<G4endl;   
1728
1729   return result;
1730
1731} 
1732
1733///////////////////////////////////////////////////////////////////////////////
1734//
1735// Integration of Plasmon cross-section for the case of
1736// passing across border between intervals
1737
1738G4double G4PAIxSection::SumOverBordPlasmon( G4int      i , 
1739                                             G4double en0    )
1740{               
1741   G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1742
1743   e0 = en0;
1744   x0 = fSplineEnergy[i];
1745   x1 = fSplineEnergy[i+1];
1746   y0 = fdNdxPlasmon[i];
1747   yy1 = fdNdxPlasmon[i+1];
1748
1749   c = x1/x0;
1750   d = e0/x0;   
1751   a = log10(yy1/y0)/log10(c);
1752   //  b0 = log10(y0) - a*log10(x0);
1753   b = y0/pow(x0,a); //pow(10.,b);
1754   
1755   a += 1.0;
1756   if( a == 0 ) result = b*log(x0/e0);
1757   else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1758   a += 1.0;
1759
1760   if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
1761   else         fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1762   
1763   x0 = fSplineEnergy[i - 1];
1764   x1 = fSplineEnergy[i - 2];
1765   y0 = fdNdxPlasmon[i - 1];
1766   yy1 = fdNdxPlasmon[i - 2];
1767
1768   c = x1/x0;
1769   d = e0/x0;
1770   a = log10(yy1/y0)/log10(c);
1771   // b0 = log10(y0) - a*log10(x0);
1772   b = y0/pow(x0,a);// pow(10.,b0);
1773
1774   a += 1.0;
1775   if( a == 0 ) result += b*log(e0/x0);
1776   else         result += y0*(e0*pow(d,a-1) - x0)/a;
1777   a += 1.0;
1778
1779   if( a == 0 )   fIntegralPlasmon[0] += b*log(e0/x0);
1780   else           fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1781   
1782   return result;
1783
1784} 
1785
1786///////////////////////////////////////////////////////////////////////////////
1787//
1788// Integration of resonance cross-section for the case of
1789// passing across border between intervals
1790
1791G4double G4PAIxSection::SumOverBordResonance( G4int      i , 
1792                                             G4double en0    )
1793{               
1794   G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
1795
1796   e0 = en0;
1797   x0 = fSplineEnergy[i];
1798   x1 = fSplineEnergy[i+1];
1799   y0 = fdNdxResonance[i];
1800   yy1 = fdNdxResonance[i+1];
1801
1802   c = x1/x0;
1803   d = e0/x0;   
1804   a = log10(yy1/y0)/log10(c);
1805   //  b0 = log10(y0) - a*log10(x0);
1806   b = y0/pow(x0,a); //pow(10.,b);
1807   
1808   a += 1.0;
1809   if( a == 0 ) result = b*log(x0/e0);
1810   else         result = y0*(x0 - e0*pow(d,a-1))/a;   
1811   a += 1.0;
1812
1813   if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
1814   else         fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1815   
1816   x0 = fSplineEnergy[i - 1];
1817   x1 = fSplineEnergy[i - 2];
1818   y0 = fdNdxResonance[i - 1];
1819   yy1 = fdNdxResonance[i - 2];
1820
1821   c = x1/x0;
1822   d = e0/x0;
1823   a = log10(yy1/y0)/log10(c);
1824   // b0 = log10(y0) - a*log10(x0);
1825   b = y0/pow(x0,a);// pow(10.,b0);
1826
1827   a += 1.0;
1828   if( a == 0 ) result += b*log(e0/x0);
1829   else         result += y0*(e0*pow(d,a-1) - x0)/a;
1830   a += 1.0;
1831
1832   if( a == 0 )   fIntegralResonance[0] += b*log(e0/x0);
1833   else           fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1834   
1835   return result;
1836
1837} 
1838
1839/////////////////////////////////////////////////////////////////////////
1840//
1841// Returns random PAI-total energy loss over step
1842
1843G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
1844{ 
1845  G4long numOfCollisions;
1846  G4double meanNumber, loss = 0.0;
1847
1848  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
1849
1850  meanNumber = fIntegralPAIxSection[1]*step;
1851  numOfCollisions = G4Poisson(meanNumber);
1852
1853  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1854
1855  while(numOfCollisions)
1856  {
1857   loss += GetEnergyTransfer();
1858   numOfCollisions--;
1859  }
1860  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
1861
1862  return loss;
1863}
1864
1865/////////////////////////////////////////////////////////////////////////
1866//
1867// Returns random PAI-total energy transfer in one collision
1868
1869G4double G4PAIxSection::GetEnergyTransfer()
1870{ 
1871  G4int iTransfer ;
1872
1873  G4double energyTransfer, position;
1874
1875  position = fIntegralPAIxSection[1]*G4UniformRand();
1876
1877  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
1878  {
1879        if( position >= fIntegralPAIxSection[iTransfer] ) break;
1880  }
1881  if(iTransfer > fSplineNumber) iTransfer--;
1882 
1883  energyTransfer = fSplineEnergy[iTransfer];
1884
1885  if(iTransfer > 1)
1886  {
1887    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
1888  }
1889  return energyTransfer;
1890}
1891
1892/////////////////////////////////////////////////////////////////////////
1893//
1894// Returns random Cerenkov energy loss over step
1895
1896G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
1897{ 
1898  G4long numOfCollisions;
1899  G4double meanNumber, loss = 0.0;
1900
1901  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
1902
1903  meanNumber = fIntegralCerenkov[1]*step;
1904  numOfCollisions = G4Poisson(meanNumber);
1905
1906  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1907
1908  while(numOfCollisions)
1909  {
1910    loss += GetCerenkovEnergyTransfer();
1911    numOfCollisions--;
1912  }
1913  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1914
1915  return loss;
1916}
1917
1918/////////////////////////////////////////////////////////////////////////
1919//
1920// Returns random MM-Cerenkov energy loss over step
1921
1922G4double G4PAIxSection::GetStepMMLoss( G4double step )
1923{ 
1924  G4long numOfCollisions;
1925  G4double meanNumber, loss = 0.0;
1926
1927  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
1928
1929  meanNumber = fIntegralMM[1]*step;
1930  numOfCollisions = G4Poisson(meanNumber);
1931
1932  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1933
1934  while(numOfCollisions)
1935  {
1936    loss += GetMMEnergyTransfer();
1937    numOfCollisions--;
1938  }
1939  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
1940
1941  return loss;
1942}
1943
1944/////////////////////////////////////////////////////////////////////////
1945//
1946// Returns Cerenkov energy transfer in one collision
1947
1948G4double G4PAIxSection::GetCerenkovEnergyTransfer()
1949{ 
1950  G4int iTransfer ;
1951
1952  G4double energyTransfer, position;
1953
1954  position = fIntegralCerenkov[1]*G4UniformRand();
1955
1956  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
1957  {
1958        if( position >= fIntegralCerenkov[iTransfer] ) break;
1959  }
1960  if(iTransfer > fSplineNumber) iTransfer--;
1961 
1962  energyTransfer = fSplineEnergy[iTransfer];
1963
1964  if(iTransfer > 1)
1965  {
1966    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
1967  }
1968  return energyTransfer;
1969}
1970
1971/////////////////////////////////////////////////////////////////////////
1972//
1973// Returns MM-Cerenkov energy transfer in one collision
1974
1975G4double G4PAIxSection::GetMMEnergyTransfer()
1976{ 
1977  G4int iTransfer ;
1978
1979  G4double energyTransfer, position;
1980
1981  position = fIntegralMM[1]*G4UniformRand();
1982
1983  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
1984  {
1985        if( position >= fIntegralMM[iTransfer] ) break;
1986  }
1987  if(iTransfer > fSplineNumber) iTransfer--;
1988 
1989  energyTransfer = fSplineEnergy[iTransfer];
1990
1991  if(iTransfer > 1)
1992  {
1993    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
1994  }
1995  return energyTransfer;
1996}
1997
1998/////////////////////////////////////////////////////////////////////////
1999//
2000// Returns random plasmon energy loss over step
2001
2002G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
2003{ 
2004  G4long numOfCollisions;
2005  G4double  meanNumber, loss = 0.0;
2006
2007  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2008
2009  meanNumber = fIntegralPlasmon[1]*step;
2010  numOfCollisions = G4Poisson(meanNumber);
2011
2012  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2013
2014  while(numOfCollisions)
2015  {
2016    loss += GetPlasmonEnergyTransfer();
2017    numOfCollisions--;
2018  }
2019  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2020
2021  return loss;
2022}
2023
2024/////////////////////////////////////////////////////////////////////////
2025//
2026// Returns plasmon energy transfer in one collision
2027
2028G4double G4PAIxSection::GetPlasmonEnergyTransfer()
2029{ 
2030  G4int iTransfer ;
2031
2032  G4double energyTransfer, position;
2033
2034  position = fIntegralPlasmon[1]*G4UniformRand();
2035
2036  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2037  {
2038        if( position >= fIntegralPlasmon[iTransfer] ) break;
2039  }
2040  if(iTransfer > fSplineNumber) iTransfer--;
2041 
2042  energyTransfer = fSplineEnergy[iTransfer];
2043
2044  if(iTransfer > 1)
2045  {
2046    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2047  }
2048  return energyTransfer;
2049}
2050
2051/////////////////////////////////////////////////////////////////////////
2052//
2053// Returns random resonance energy loss over step
2054
2055G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
2056{ 
2057  G4long numOfCollisions;
2058  G4double meanNumber, loss = 0.0;
2059
2060  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2061
2062  meanNumber = fIntegralResonance[1]*step;
2063  numOfCollisions = G4Poisson(meanNumber);
2064
2065  //   G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2066
2067  while(numOfCollisions)
2068  {
2069    loss += GetResonanceEnergyTransfer();
2070    numOfCollisions--;
2071  }
2072  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2073
2074  return loss;
2075}
2076
2077
2078/////////////////////////////////////////////////////////////////////////
2079//
2080// Returns resonance energy transfer in one collision
2081
2082G4double G4PAIxSection::GetResonanceEnergyTransfer()
2083{ 
2084  G4int iTransfer ;
2085
2086  G4double energyTransfer, position;
2087
2088  position = fIntegralResonance[1]*G4UniformRand();
2089
2090  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2091  {
2092        if( position >= fIntegralResonance[iTransfer] ) break;
2093  }
2094  if(iTransfer > fSplineNumber) iTransfer--;
2095 
2096  energyTransfer = fSplineEnergy[iTransfer];
2097
2098  if(iTransfer > 1)
2099  {
2100    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2101  }
2102  return energyTransfer;
2103}
2104
2105
2106/////////////////////////////////////////////////////////////////////////
2107//
2108// Returns Rutherford energy transfer in one collision
2109
2110G4double G4PAIxSection::GetRutherfordEnergyTransfer()
2111{ 
2112  G4int iTransfer ;
2113
2114  G4double energyTransfer, position;
2115
2116  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2117
2118  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2119  {
2120        if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2121  }
2122  if(iTransfer > fSplineNumber) iTransfer--;
2123 
2124  energyTransfer = fSplineEnergy[iTransfer];
2125
2126  if(iTransfer > 1)
2127  {
2128    energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2129  }
2130  return energyTransfer;
2131}
2132
2133
2134/////////////////////////////////////////////////////////////////////////////
2135//
2136// Init  array of Lorentz factors
2137//
2138
2139G4int G4PAIxSection::fNumberOfGammas = 111;
2140
2141const G4double G4PAIxSection::fLorentzFactor[112] =     // fNumberOfGammas+1
2142{
21430.0,
21441.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
21451.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
21461.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
21471.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
21482.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
21493.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
21505.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
21518.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
21521.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
21532.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
21545.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
21551.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
21561.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
21573.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
21586.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
21591.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
21602.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
21614.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
21628.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
21631.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
21643.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
21655.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
21661.065799e+05
2167};
2168
2169///////////////////////////////////////////////////////////////////////
2170//
2171// The number of gamma for creation of  spline (near ion-min , G ~ 4 )
2172//
2173
2174const
2175G4int G4PAIxSection::fRefGammaNumber = 29; 
2176
2177   
2178//   
2179// end of G4PAIxSection implementation file
2180//
2181////////////////////////////////////////////////////////////////////////////
2182
Note: See TracBrowser for help on using the repository browser.