source: trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc@ 899

Last change on this file since 899 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 47.6 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.21 2006/06/29 19:53:20 gunter Exp $
28// GEANT4 tag $Name: $
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 fSandia = new G4SandiaTable(matIndex);
98
99 G4int i, j;
100 fMatSandiaMatrix = new G4OrderedTable();
101
102 for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
103 {
104 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
105 }
106 for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
107 {
108 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
109
110 for(j = 1; j < 5 ; j++)
111 {
112 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
113 }
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 fDensity = (*theMaterialTable)[materialIndex]->GetDensity() ;
127 fElectronDensity = (*theMaterialTable)[materialIndex]->
128 GetElectronDensity() ;
129 fIntervalNumber = (*theMaterialTable)[materialIndex]->
130 GetSandiaTable()->GetMatNbOfIntervals() ;
131 fIntervalNumber--;
132 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl ;
133
134 fEnergyInterval = new G4double[fIntervalNumber+2] ;
135 fA1 = new G4double[fIntervalNumber+2] ;
136 fA2 = new G4double[fIntervalNumber+2] ;
137 fA3 = new G4double[fIntervalNumber+2] ;
138 fA4 = new G4double[fIntervalNumber+2] ;
139
140 for(i = 1 ; i <= fIntervalNumber ; i++ )
141 {
142 if(((*theMaterialTable)[materialIndex]->
143 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
144 i > fIntervalNumber )
145 {
146 fEnergyInterval[i] = maxEnergyTransfer ;
147 fIntervalNumber = i ;
148 break;
149 }
150 fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
151 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
152 fA1[i] = (*theMaterialTable)[materialIndex]->
153 GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
154 fA2[i] = (*theMaterialTable)[materialIndex]->
155 GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
156 fA3[i] = (*theMaterialTable)[materialIndex]->
157 GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
158 fA4[i] = (*theMaterialTable)[materialIndex]->
159 GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
160 // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
161 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
162 }
163 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
164 {
165 fIntervalNumber++;
166 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
167 }
168
169 // Now checking, if two borders are too close together
170
171 for(i=1;i<fIntervalNumber;i++)
172 {
173 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
174 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
175 {
176 continue ;
177 }
178 else
179 {
180 for(j=i;j<fIntervalNumber;j++)
181 {
182 fEnergyInterval[j] = fEnergyInterval[j+1] ;
183 fA1[j] = fA1[j+1] ;
184 fA2[j] = fA2[j+1] ;
185 fA3[j] = fA3[j+1] ;
186 fA4[j] = fA4[j+1] ;
187 }
188 fIntervalNumber-- ;
189 i-- ;
190 }
191 }
192
193
194 /* *********************************
195
196 fSplineEnergy = new G4double[fMaxSplineSize] ;
197 fRePartDielectricConst = new G4double[fMaxSplineSize] ;
198 fImPartDielectricConst = new G4double[fMaxSplineSize] ;
199 fIntegralTerm = new G4double[fMaxSplineSize] ;
200 fDifPAIxSection = new G4double[fMaxSplineSize] ;
201 fIntegralPAIxSection = new G4double[fMaxSplineSize] ;
202
203 for(i=0;i<fMaxSplineSize;i++)
204 {
205 fSplineEnergy[i] = 0.0 ;
206 fRePartDielectricConst[i] = 0.0 ;
207 fImPartDielectricConst[i] = 0.0 ;
208 fIntegralTerm[i] = 0.0 ;
209 fDifPAIxSection[i] = 0.0 ;
210 fIntegralPAIxSection[i] = 0.0 ;
211 }
212 ************************************************** */
213
214 InitPAI() ; // create arrays allocated above
215
216 delete[] fEnergyInterval ;
217 delete[] fA1 ;
218 delete[] fA2 ;
219 delete[] fA3 ;
220 delete[] fA4 ;
221}
222
223////////////////////////////////////////////////////////////////////////
224//
225// Constructor with beta*gamma square value
226
227G4PAIxSection::G4PAIxSection( G4int materialIndex,
228 G4double maxEnergyTransfer,
229 G4double betaGammaSq,
230 G4double** photoAbsCof,
231 G4int intNumber )
232{
233 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
234 G4int i, j ;
235
236 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
237 fElectronDensity = (*theMaterialTable)[materialIndex]->
238 GetElectronDensity() ;
239
240 fIntervalNumber = intNumber ;
241 fIntervalNumber--;
242 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl ;
243
244 fEnergyInterval = new G4double[fIntervalNumber+2] ;
245 fA1 = new G4double[fIntervalNumber+2] ;
246 fA2 = new G4double[fIntervalNumber+2] ;
247 fA3 = new G4double[fIntervalNumber+2] ;
248 fA4 = new G4double[fIntervalNumber+2] ;
249
250 for( i = 1 ; i <= fIntervalNumber ; i++ )
251 {
252 if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
253 i > fIntervalNumber )
254 {
255 fEnergyInterval[i] = maxEnergyTransfer ;
256 fIntervalNumber = i ;
257 break;
258 }
259 fEnergyInterval[i] = photoAbsCof[i-1][0] ;
260 fA1[i] = photoAbsCof[i-1][1] ;
261 fA2[i] = photoAbsCof[i-1][2] ;
262 fA3[i] = photoAbsCof[i-1][3] ;
263 fA4[i] = photoAbsCof[i-1][4] ;
264 // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
265 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
266 }
267 // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
268 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
269 {
270 fIntervalNumber++;
271 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
272 }
273 for(i=1;i<=fIntervalNumber;i++)
274 {
275 // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
276 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
277 }
278 // Now checking, if two borders are too close together
279
280 for( i = 1 ; i < fIntervalNumber ; i++ )
281 {
282 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
283 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
284 {
285 continue ;
286 }
287 else
288 {
289 for(j=i;j<fIntervalNumber;j++)
290 {
291 fEnergyInterval[j] = fEnergyInterval[j+1] ;
292 fA1[j] = fA1[j+1] ;
293 fA2[j] = fA2[j+1] ;
294 fA3[j] = fA3[j+1] ;
295 fA4[j] = fA4[j+1] ;
296 }
297 fIntervalNumber-- ;
298 i-- ;
299 }
300 }
301
302 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
303
304 G4double betaGammaSqRef =
305 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
306
307 NormShift(betaGammaSqRef) ;
308 SplainPAI(betaGammaSqRef) ;
309
310 // Preparation of integral PAI cross section for input betaGammaSq
311
312 for(i = 1 ; i <= fSplineNumber ; i++)
313 {
314 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
315 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
316 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
317 // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
318 // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
319 }
320 IntegralCerenkov() ;
321 IntegralPlasmon() ;
322 IntegralPAIxSection() ;
323
324 delete[] fEnergyInterval ;
325 delete[] fA1 ;
326 delete[] fA2 ;
327 delete[] fA3 ;
328 delete[] fA4 ;
329}
330
331////////////////////////////////////////////////////////////////////////
332//
333// Test Constructor with beta*gamma square value
334
335G4PAIxSection::G4PAIxSection( G4int materialIndex,
336 G4double maxEnergyTransfer,
337 G4double betaGammaSq )
338{
339 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
340
341 G4int i, j, numberOfElements ;
342
343 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
344 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity() ;
345 numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements() ;
346
347 G4int* thisMaterialZ = new G4int[numberOfElements] ;
348
349 for( i = 0 ; i < numberOfElements ; i++ )
350 {
351 thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
352 GetElement(i)->GetZ() ;
353 }
354 G4SandiaTable thisMaterialSandiaTable(materialIndex) ;
355 fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
356 (thisMaterialZ,numberOfElements) ;
357 fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
358 ( thisMaterialZ ,
359 (*theMaterialTable)[materialIndex]->GetFractionVector() ,
360 numberOfElements,fIntervalNumber) ;
361
362 fIntervalNumber--;
363
364 fEnergyInterval = new G4double[fIntervalNumber+2] ;
365 fA1 = new G4double[fIntervalNumber+2] ;
366 fA2 = new G4double[fIntervalNumber+2] ;
367 fA3 = new G4double[fIntervalNumber+2] ;
368 fA4 = new G4double[fIntervalNumber+2] ;
369
370 for(i=1;i<=fIntervalNumber;i++)
371 {
372 if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
373 i > fIntervalNumber)
374 {
375 fEnergyInterval[i] = maxEnergyTransfer ;
376 fIntervalNumber = i ;
377 break;
378 }
379 fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) ;
380 fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity ;
381 fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity ;
382 fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity ;
383 fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity ;
384
385 }
386 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
387 {
388 fIntervalNumber++;
389 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ;
390 fA1[fIntervalNumber] = fA1[fIntervalNumber-1] ;
391 fA2[fIntervalNumber] = fA2[fIntervalNumber-1] ;
392 fA3[fIntervalNumber] = fA3[fIntervalNumber-1] ;
393 fA4[fIntervalNumber] = fA4[fIntervalNumber-1] ;
394 }
395 for(i=1;i<=fIntervalNumber;i++)
396 {
397 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
398 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
399 }
400 // Now checking, if two borders are too close together
401
402 for(i=1;i<fIntervalNumber;i++)
403 {
404 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
405 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
406 {
407 continue ;
408 }
409 else
410 {
411 for(j=i;j<fIntervalNumber;j++)
412 {
413 fEnergyInterval[j] = fEnergyInterval[j+1] ;
414 fA1[j] = fA1[j+1] ;
415 fA2[j] = fA2[j+1] ;
416 fA3[j] = fA3[j+1] ;
417 fA4[j] = fA4[j+1] ;
418 }
419 fIntervalNumber-- ;
420 i-- ;
421 }
422 }
423
424 /* *********************************
425 fSplineEnergy = new G4double[fMaxSplineSize] ;
426 fRePartDielectricConst = new G4double[fMaxSplineSize] ;
427 fImPartDielectricConst = new G4double[fMaxSplineSize] ;
428 fIntegralTerm = new G4double[fMaxSplineSize] ;
429 fDifPAIxSection = new G4double[fMaxSplineSize] ;
430 fIntegralPAIxSection = new G4double[fMaxSplineSize] ;
431
432 for(i=0;i<fMaxSplineSize;i++)
433 {
434 fSplineEnergy[i] = 0.0 ;
435 fRePartDielectricConst[i] = 0.0 ;
436 fImPartDielectricConst[i] = 0.0 ;
437 fIntegralTerm[i] = 0.0 ;
438 fDifPAIxSection[i] = 0.0 ;
439 fIntegralPAIxSection[i] = 0.0 ;
440 }
441 */ ////////////////////////
442
443 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
444
445 G4double betaGammaSqRef =
446 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
447
448 NormShift(betaGammaSqRef) ;
449 SplainPAI(betaGammaSqRef) ;
450
451 // Preparation of integral PAI cross section for input betaGammaSq
452
453 for(i = 1 ; i <= fSplineNumber ; i++)
454 {
455 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
456 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
457 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
458 }
459 IntegralPAIxSection() ;
460 IntegralCerenkov() ;
461 IntegralPlasmon() ;
462
463 // delete[] fEnergyInterval ;
464 delete[] fA1 ;
465 delete[] fA2 ;
466 delete[] fA3 ;
467 delete[] fA4 ;
468}
469
470
471////////////////////////////////////////////////////////////////////////////
472//
473// Destructor
474
475G4PAIxSection::~G4PAIxSection()
476{
477 /* ************************
478 delete[] fSplineEnergy ;
479 delete[] fRePartDielectricConst ;
480 delete[] fImPartDielectricConst ;
481 delete[] fIntegralTerm ;
482 delete[] fDifPAIxSection ;
483 delete[] fIntegralPAIxSection ;
484 */ ////////////////////////
485}
486
487/////////////////////////////////////////////////////////////////////////
488//
489// General control function for class G4PAIxSection
490//
491
492void G4PAIxSection::InitPAI()
493{
494 G4int i ;
495 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
496 fLorentzFactor[fRefGammaNumber] - 1;
497
498 // Preparation of integral PAI cross section for reference gamma
499
500 NormShift(betaGammaSq) ;
501 SplainPAI(betaGammaSq) ;
502
503 IntegralPAIxSection() ;
504 IntegralCerenkov() ;
505 IntegralPlasmon() ;
506
507 for(i = 0 ; i<=fSplineNumber ; i++)
508 {
509 fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i] ;
510 if(i != 0)
511 {
512 fPAItable[i][0] = fSplineEnergy[i] ;
513 }
514 }
515 fPAItable[0][0] = fSplineNumber ;
516
517 for(G4int j = 1 ; j < 112 ; j++) // for other gammas
518 {
519 if( j == fRefGammaNumber ) continue ;
520
521 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1 ;
522
523 for(i = 1 ; i <= fSplineNumber ; i++)
524 {
525 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
526 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
527 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
528 }
529 IntegralPAIxSection() ;
530 IntegralCerenkov() ;
531 IntegralPlasmon() ;
532
533 for(i = 0 ; i <= fSplineNumber ; i++)
534 {
535 fPAItable[i][j] = fIntegralPAIxSection[i] ;
536 }
537 }
538
539}
540
541///////////////////////////////////////////////////////////////////////
542//
543// Shifting from borders to intervals Creation of first energy points
544//
545
546void G4PAIxSection::NormShift(G4double betaGammaSq)
547{
548 G4int i, j ;
549
550 for( i = 1 ; i <= fIntervalNumber-1 ; i++ )
551 {
552 for( j = 1 ; j <= 2 ; j++ )
553 {
554 fSplineNumber = (i-1)*2 + j ;
555
556 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
557 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
558 // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = "
559 // <<fSplineEnergy[fSplineNumber]<<G4endl;
560 }
561 }
562 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
563
564 j = 1 ;
565
566 for(i=2;i<=fSplineNumber;i++)
567 {
568 if(fSplineEnergy[i]<fEnergyInterval[j+1])
569 {
570 fIntegralTerm[i] = fIntegralTerm[i-1] +
571 RutherfordIntegral(j,fSplineEnergy[i-1],
572 fSplineEnergy[i] ) ;
573 }
574 else
575 {
576 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
577 fEnergyInterval[j+1] ) ;
578 j++;
579 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
580 RutherfordIntegral(j,fEnergyInterval[j],
581 fSplineEnergy[i] ) ;
582 }
583 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl;
584 }
585 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
586 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber] ;
587
588 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl ;
589
590 // Calculation of PAI differrential cross-section (1/(keV*cm))
591 // in the energy points near borders of energy intervals
592
593 for(G4int k=1;k<=fIntervalNumber-1;k++)
594 {
595 for(j=1;j<=2;j++)
596 {
597 i = (k-1)*2 + j ;
598 fImPartDielectricConst[i] = fNormalizationCof*
599 ImPartDielectricConst(k,fSplineEnergy[i]);
600 fRePartDielectricConst[i] = fNormalizationCof*
601 RePartDielectricConst(fSplineEnergy[i]);
602 fIntegralTerm[i] *= fNormalizationCof;
603
604 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
605 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
606 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
607 }
608 }
609
610} // end of NormShift
611
612/////////////////////////////////////////////////////////////////////////
613//
614// Creation of new energy points as geometrical mean of existing
615// one, calculation PAI_cs for them, while the error of logarithmic
616// linear approximation would be smaller than 'fError'
617
618void
619 G4PAIxSection::SplainPAI(G4double betaGammaSq)
620{
621 G4int k = 1 ;
622 G4int i = 1 ;
623
624 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
625 {
626 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
627 {
628 k++ ; // Here next energy point is in next energy interval
629 i++;
630 continue;
631 }
632 // Shifting of arrayes for inserting the geometrical
633 // average of 'i' and 'i+1' energy points to 'i+1' place
634 fSplineNumber++;
635
636 for(G4int j = fSplineNumber; j >= i+2 ; j-- )
637 {
638 fSplineEnergy[j] = fSplineEnergy[j-1];
639 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
640 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
641 fIntegralTerm[j] = fIntegralTerm[j-1];
642
643 fDifPAIxSection[j] = fDifPAIxSection[j-1];
644 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
645 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
646 }
647 G4double x1 = fSplineEnergy[i];
648 G4double x2 = fSplineEnergy[i+1];
649 G4double yy1 = fDifPAIxSection[i];
650 G4double y2 = fDifPAIxSection[i+1];
651
652 G4double en1 = sqrt(x1*x2);
653 fSplineEnergy[i+1] = en1;
654
655 // Calculation of logarithmic linear approximation
656 // in this (enr) energy point, which number is 'i+1' now
657
658 G4double a = log10(y2/yy1)/log10(x2/x1);
659 G4double b = log10(yy1) - a*log10(x1);
660 G4double y = a*log10(en1) + b ;
661 y = pow(10.,y);
662
663 // Calculation of the PAI dif. cross-section at this point
664
665 fImPartDielectricConst[i+1] = fNormalizationCof*
666 ImPartDielectricConst(k,fSplineEnergy[i+1]);
667 fRePartDielectricConst[i+1] = fNormalizationCof*
668 RePartDielectricConst(fSplineEnergy[i+1]);
669 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
670 RutherfordIntegral(k,fSplineEnergy[i],
671 fSplineEnergy[i+1]);
672
673 fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
674 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
675 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
676
677 // Condition for next division of this segment or to pass
678 // to higher energies
679
680 G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
681
682 if( x < 0 )
683 {
684 x = -x ;
685 }
686 if( x > fError && fSplineNumber < fMaxSplineSize-1 )
687 {
688 continue; // next division
689 }
690 i += 2; // pass to next segment
691
692 } // close 'while'
693
694} // end of SplainPAI
695
696
697////////////////////////////////////////////////////////////////////
698//
699// Integration over electrons that could be considered
700// quasi-free at energy transfer of interest
701
702G4double G4PAIxSection::RutherfordIntegral( G4int k,
703 G4double x1,
704 G4double x2 )
705{
706 G4double c1, c2, c3 ;
707 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
708 c1 = (x2 - x1)/x1/x2 ;
709 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
710 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
711 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
712
713 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3 ;
714
715} // end of RutherfordIntegral
716
717
718/////////////////////////////////////////////////////////////////
719//
720// Imaginary part of dielectric constant
721// (G4int k - interval number, G4double en1 - energy point)
722
723G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
724 G4double energy1 )
725{
726 G4double energy2,energy3,energy4,result;
727
728 energy2 = energy1*energy1;
729 energy3 = energy2*energy1;
730 energy4 = energy3*energy1;
731
732 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4 ;
733 result *=hbarc/energy1 ;
734
735 return result ;
736
737} // end of ImPartDielectricConst
738
739
740//////////////////////////////////////////////////////////////////////////////
741//
742// Real part of dielectric constant minus unit: epsilon_1 - 1
743// (G4double enb - energy point)
744//
745
746G4double G4PAIxSection::RePartDielectricConst(G4double enb)
747{
748 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
749 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
750
751 x0 = enb ;
752 result = 0 ;
753
754 for(G4int i=1;i<=fIntervalNumber-1;i++)
755 {
756 x1 = fEnergyInterval[i] ;
757 x2 = fEnergyInterval[i+1] ;
758 xx1 = x1 - x0 ;
759 xx2 = x2 - x0 ;
760 xx12 = xx2/xx1 ;
761
762 if(xx12<0)
763 {
764 xx12 = -xx12;
765 }
766 xln1 = log(x2/x1) ;
767 xln2 = log(xx12) ;
768 xln3 = log((x2 + x0)/(x1 + x0)) ;
769 x02 = x0*x0 ;
770 x03 = x02*x0 ;
771 x04 = x03*x0 ;
772 x05 = x04*x0;
773 c1 = (x2 - x1)/x1/x2 ;
774 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
775 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
776
777 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 ;
778 result -= (fA2[i]/x02 + fA4[i]/x04)*c1 ;
779 result -= fA3[i]*c2/2/x02 ;
780 result -= fA4[i]*c3/3/x02 ;
781
782 cof1 = fA1[i]/x02 + fA3[i]/x04 ;
783 cof2 = fA2[i]/x03 + fA4[i]/x05 ;
784
785 result += 0.5*(cof1 +cof2)*xln2 ;
786 result += 0.5*(cof1 - cof2)*xln3 ;
787 }
788 result *= 2*hbarc/pi ;
789
790 return result ;
791
792} // end of RePartDielectricConst
793
794//////////////////////////////////////////////////////////////////////
795//
796// PAI differential cross-section in terms of
797// simplified Allison's equation
798//
799
800G4double G4PAIxSection::DifPAIxSection( G4int i ,
801 G4double betaGammaSq )
802{
803 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
804 //G4double beta, be4 ;
805 G4double be4 ;
806 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
807 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
808 be2 = betaGammaSq/(1 + betaGammaSq) ;
809 be4 = be2*be2 ;
810 // beta = sqrt(be2) ;
811 cof = 1 ;
812 x1 = log(2*electron_mass_c2/fSplineEnergy[i]) ;
813
814 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
815 else
816 {
817 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
818 (1/betaGammaSq - fRePartDielectricConst[i]) +
819 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2 ;
820 }
821 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
822 {
823 x6=0 ;
824 }
825 else
826 {
827 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq ;
828 x5 = -1 - fRePartDielectricConst[i] +
829 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
830 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
831
832 x7 = atan2(fImPartDielectricConst[i],x3) ;
833 x6 = x5 * x7 ;
834 }
835 // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
836
837 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc ;
838 // if( x4 < 0.0 ) x4 = 0.0 ;
839 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
840 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
841
842 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]) ;
843 if(result < 1.0e-8) result = 1.0e-8 ;
844 result *= fine_structure_const/be2/pi ;
845 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
846 // result *= (1-exp(-be2/betaBohr2)) ;
847 result *= (1-exp(-be4/betaBohr4)) ;
848 if(fDensity >= 0.1)
849 {
850 result /= x8 ;
851 }
852 return result ;
853
854} // end of DifPAIxSection
855
856//////////////////////////////////////////////////////////////////////////
857//
858// Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons
859
860G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
861 G4double betaGammaSq )
862{
863 G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ;
864 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
865
866 cof = 1.0 ;
867 cofBetaBohr = 4.0 ;
868 betaBohr2 = fine_structure_const*fine_structure_const ;
869 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
870
871 be2 = betaGammaSq/(1 + betaGammaSq) ;
872 be4 = be2*be2 ;
873
874 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
875 else
876 {
877 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
878 (1/betaGammaSq - fRePartDielectricConst[i]) +
879 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5 ;
880 logarithm += log(1+1.0/betaGammaSq) ;
881 }
882
883 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
884 {
885 argument = 0.0 ;
886 }
887 else
888 {
889 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq ;
890 x5 = -1.0 - fRePartDielectricConst[i] +
891 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
892 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ;
893 if( x3 == 0.0 ) argument = 0.5*pi;
894 else argument = atan2(fImPartDielectricConst[i],x3) ;
895 argument *= x5 ;
896 }
897 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc ;
898
899 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
900
901 dNdxC *= fine_structure_const/be2/pi ;
902
903 dNdxC *= (1-exp(-be4/betaBohr4)) ;
904
905 if(fDensity >= 0.1)
906 {
907 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
908 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
909 dNdxC /= modul2 ;
910 }
911 return dNdxC ;
912
913} // end of PAIdNdxCerenkov
914
915//////////////////////////////////////////////////////////////////////////
916//
917// Calculation od dN/dx of collisions with creation of longitudinal EM
918// excitations (plasmons, delta-electrons)
919
920G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
921 G4double betaGammaSq )
922{
923 G4double cof, resonance, modul2, dNdxP ;
924 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
925
926 cof = 1 ;
927 cofBetaBohr = 4.0 ;
928 betaBohr2 = fine_structure_const*fine_structure_const ;
929 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
930
931 be2 = betaGammaSq/(1 + betaGammaSq) ;
932 be4 = be2*be2 ;
933
934 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]) ;
935 resonance *= fImPartDielectricConst[i]/hbarc ;
936
937
938 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ) ;
939
940 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
941
942 dNdxP *= fine_structure_const/be2/pi ;
943 dNdxP *= (1-exp(-be4/betaBohr4)) ;
944
945 if( fDensity >= 0.1 )
946 {
947 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
948 fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
949 dNdxP /= modul2 ;
950 }
951 return dNdxP ;
952
953} // end of PAIdNdxPlasmon
954
955////////////////////////////////////////////////////////////////////////
956//
957// Calculation of the PAI integral cross-section
958// fIntegralPAIxSection[1] = specific primary ionisation, 1/cm
959// and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm
960
961void G4PAIxSection::IntegralPAIxSection()
962{
963 fIntegralPAIxSection[fSplineNumber] = 0 ;
964 fIntegralPAIdEdx[fSplineNumber] = 0 ;
965 fIntegralPAIxSection[0] = 0 ;
966 G4int k = fIntervalNumber -1 ;
967
968 for(G4int i = fSplineNumber-1 ; i >= 1 ; i--)
969 {
970 if(fSplineEnergy[i] >= fEnergyInterval[k])
971 {
972 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i) ;
973 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i) ;
974 }
975 else
976 {
977 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
978 SumOverBorder(i+1,fEnergyInterval[k]) ;
979 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
980 SumOverBorderdEdx(i+1,fEnergyInterval[k]) ;
981 k-- ;
982 }
983 }
984} // end of IntegralPAIxSection
985
986////////////////////////////////////////////////////////////////////////
987//
988// Calculation of the PAI Cerenkov integral cross-section
989// fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm
990// and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm
991
992void G4PAIxSection::IntegralCerenkov()
993{
994 G4int i, k ;
995 fIntegralCerenkov[fSplineNumber] = 0 ;
996 fIntegralCerenkov[0] = 0 ;
997 k = fIntervalNumber -1 ;
998
999 for( i = fSplineNumber-1 ; i >= 1 ; i-- )
1000 {
1001 if(fSplineEnergy[i] >= fEnergyInterval[k])
1002 {
1003 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i) ;
1004 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1005 }
1006 else
1007 {
1008 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1009 SumOverBordCerenkov(i+1,fEnergyInterval[k]) ;
1010 k-- ;
1011 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1012 }
1013 }
1014
1015} // end of IntegralCerenkov
1016
1017////////////////////////////////////////////////////////////////////////
1018//
1019// Calculation of the PAI Plasmon integral cross-section
1020// fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm
1021// and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm
1022
1023void G4PAIxSection::IntegralPlasmon()
1024{
1025 fIntegralPlasmon[fSplineNumber] = 0 ;
1026 fIntegralPlasmon[0] = 0 ;
1027 G4int k = fIntervalNumber -1 ;
1028 for(G4int i=fSplineNumber-1;i>=1;i--)
1029 {
1030 if(fSplineEnergy[i] >= fEnergyInterval[k])
1031 {
1032 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i) ;
1033 }
1034 else
1035 {
1036 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1037 SumOverBordPlasmon(i+1,fEnergyInterval[k]) ;
1038 k-- ;
1039 }
1040 }
1041
1042} // end of IntegralPlasmon
1043
1044//////////////////////////////////////////////////////////////////////
1045//
1046// Calculation the PAI integral cross-section inside
1047// of interval of continuous values of photo-ionisation
1048// cross-section. Parameter 'i' is the number of interval.
1049
1050G4double G4PAIxSection::SumOverInterval( G4int i )
1051{
1052 G4double x0,x1,y0,yy1,a,b,c,result ;
1053
1054 x0 = fSplineEnergy[i] ;
1055 x1 = fSplineEnergy[i+1] ;
1056 y0 = fDifPAIxSection[i] ;
1057 yy1 = fDifPAIxSection[i+1];
1058 c = x1/x0;
1059 a = log10(yy1/y0)/log10(c) ;
1060 // b = log10(y0) - a*log10(x0) ;
1061 b = y0/pow(x0,a) ;
1062 a += 1 ;
1063 if(a == 0)
1064 {
1065 result = b*log(x1/x0) ;
1066 }
1067 else
1068 {
1069 result = y0*(x1*pow(c,a-1) - x0)/a ;
1070 }
1071 a++;
1072 if(a == 0)
1073 {
1074 fIntegralPAIxSection[0] += b*log(x1/x0) ;
1075 }
1076 else
1077 {
1078 fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1079 }
1080 return result ;
1081
1082} // end of SumOverInterval
1083
1084/////////////////////////////////
1085
1086G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
1087{
1088 G4double x0,x1,y0,yy1,a,b,c,result ;
1089
1090 x0 = fSplineEnergy[i] ;
1091 x1 = fSplineEnergy[i+1] ;
1092 y0 = fDifPAIxSection[i] ;
1093 yy1 = fDifPAIxSection[i+1];
1094 c = x1/x0;
1095 a = log10(yy1/y0)/log10(c) ;
1096 // b = log10(y0) - a*log10(x0) ;
1097 b = y0/pow(x0,a) ;
1098 a += 2 ;
1099 if(a == 0)
1100 {
1101 result = b*log(x1/x0) ;
1102 }
1103 else
1104 {
1105 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1106 }
1107 return result ;
1108
1109} // end of SumOverInterval
1110
1111//////////////////////////////////////////////////////////////////////
1112//
1113// Calculation the PAI Cerenkov integral cross-section inside
1114// of interval of continuous values of photo-ionisation Cerenkov
1115// cross-section. Parameter 'i' is the number of interval.
1116
1117G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
1118{
1119 G4double x0,x1,y0,yy1,a,b,c,result ;
1120
1121 x0 = fSplineEnergy[i] ;
1122 x1 = fSplineEnergy[i+1] ;
1123 y0 = fdNdxCerenkov[i] ;
1124 yy1 = fdNdxCerenkov[i+1];
1125 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1126 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1127
1128 c = x1/x0;
1129 a = log10(yy1/y0)/log10(c) ;
1130 b = y0/pow(x0,a) ;
1131
1132 a += 1.0 ;
1133 if(a == 0) result = b*log(c) ;
1134 else result = y0*(x1*pow(c,a-1) - x0)/a ;
1135 a += 1.0 ;
1136
1137 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0) ;
1138 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1139 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1140 return result ;
1141
1142} // end of SumOverInterCerenkov
1143
1144//////////////////////////////////////////////////////////////////////
1145//
1146// Calculation the PAI Plasmon integral cross-section inside
1147// of interval of continuous values of photo-ionisation Plasmon
1148// cross-section. Parameter 'i' is the number of interval.
1149
1150G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
1151{
1152 G4double x0,x1,y0,yy1,a,b,c,result ;
1153
1154 x0 = fSplineEnergy[i] ;
1155 x1 = fSplineEnergy[i+1] ;
1156 y0 = fdNdxPlasmon[i] ;
1157 yy1 = fdNdxPlasmon[i+1];
1158 c =x1/x0;
1159 a = log10(yy1/y0)/log10(c) ;
1160 // b = log10(y0) - a*log10(x0) ;
1161 b = y0/pow(x0,a) ;
1162
1163 a += 1.0 ;
1164 if(a == 0) result = b*log(x1/x0) ;
1165 else result = y0*(x1*pow(c,a-1) - x0)/a ;
1166 a += 1.0 ;
1167
1168 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0) ;
1169 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ;
1170
1171 return result ;
1172
1173} // end of SumOverInterPlasmon
1174
1175///////////////////////////////////////////////////////////////////////////////
1176//
1177// Integration of PAI cross-section for the case of
1178// passing across border between intervals
1179
1180G4double G4PAIxSection::SumOverBorder( G4int i ,
1181 G4double en0 )
1182{
1183 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1184
1185 e0 = en0 ;
1186 x0 = fSplineEnergy[i] ;
1187 x1 = fSplineEnergy[i+1] ;
1188 y0 = fDifPAIxSection[i] ;
1189 yy1 = fDifPAIxSection[i+1] ;
1190
1191 c = x1/x0;
1192 d = e0/x0;
1193 a = log10(yy1/y0)/log10(x1/x0) ;
1194 // b0 = log10(y0) - a*log10(x0) ;
1195 b = y0/pow(x0,a); // pow(10.,b) ;
1196
1197 a += 1 ;
1198 if(a == 0)
1199 {
1200 result = b*log(x0/e0) ;
1201 }
1202 else
1203 {
1204 result = y0*(x0 - e0*pow(d,a-1))/a ;
1205 }
1206 a++ ;
1207 if(a == 0)
1208 {
1209 fIntegralPAIxSection[0] += b*log(x0/e0) ;
1210 }
1211 else
1212 {
1213 fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1214 }
1215 x0 = fSplineEnergy[i - 1] ;
1216 x1 = fSplineEnergy[i - 2] ;
1217 y0 = fDifPAIxSection[i - 1] ;
1218 yy1 = fDifPAIxSection[i - 2] ;
1219
1220 c = x1/x0;
1221 d = e0/x0;
1222 a = log10(yy1/y0)/log10(x1/x0) ;
1223 // b0 = log10(y0) - a*log10(x0) ;
1224 b = y0/pow(x0,a) ;
1225 a += 1 ;
1226 if(a == 0)
1227 {
1228 result += b*log(e0/x0) ;
1229 }
1230 else
1231 {
1232 result += y0*(e0*pow(d,a-1) - x0)/a ;
1233 }
1234 a++ ;
1235 if(a == 0)
1236 {
1237 fIntegralPAIxSection[0] += b*log(e0/x0) ;
1238 }
1239 else
1240 {
1241 fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1242 }
1243 return result ;
1244
1245}
1246
1247///////////////////////////////////////////////////////////////////////
1248
1249G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
1250 G4double en0 )
1251{
1252 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1253
1254 e0 = en0 ;
1255 x0 = fSplineEnergy[i] ;
1256 x1 = fSplineEnergy[i+1] ;
1257 y0 = fDifPAIxSection[i] ;
1258 yy1 = fDifPAIxSection[i+1] ;
1259
1260 c = x1/x0;
1261 d = e0/x0;
1262 a = log10(yy1/y0)/log10(x1/x0) ;
1263 // b0 = log10(y0) - a*log10(x0) ;
1264 b = y0/pow(x0,a); // pow(10.,b) ;
1265
1266 a += 2 ;
1267 if(a == 0)
1268 {
1269 result = b*log(x0/e0) ;
1270 }
1271 else
1272 {
1273 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1274 }
1275 x0 = fSplineEnergy[i - 1] ;
1276 x1 = fSplineEnergy[i - 2] ;
1277 y0 = fDifPAIxSection[i - 1] ;
1278 yy1 = fDifPAIxSection[i - 2] ;
1279
1280 c = x1/x0;
1281 d = e0/x0;
1282 a = log10(yy1/y0)/log10(x1/x0) ;
1283 // b0 = log10(y0) - a*log10(x0) ;
1284 b = y0/pow(x0,a) ;
1285 a += 2 ;
1286 if(a == 0)
1287 {
1288 result += b*log(e0/x0) ;
1289 }
1290 else
1291 {
1292 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1293 }
1294 return result ;
1295
1296}
1297
1298///////////////////////////////////////////////////////////////////////////////
1299//
1300// Integration of Cerenkov cross-section for the case of
1301// passing across border between intervals
1302
1303G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
1304 G4double en0 )
1305{
1306 G4double x0,x1,y0,yy1,a,b,e0,c,d,result ;
1307
1308 e0 = en0 ;
1309 x0 = fSplineEnergy[i] ;
1310 x1 = fSplineEnergy[i+1] ;
1311 y0 = fdNdxCerenkov[i] ;
1312 yy1 = fdNdxCerenkov[i+1] ;
1313
1314 // G4cout<<G4endl;
1315 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1316 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1317 c = x1/x0 ;
1318 d = e0/x0 ;
1319 a = log10(yy1/y0)/log10(c) ;
1320 // b0 = log10(y0) - a*log10(x0) ;
1321 b = y0/pow(x0,a); // pow(10.,b0) ;
1322
1323 a += 1.0 ;
1324 if( a == 0 ) result = b*log(x0/e0) ;
1325 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1326 a += 1.0 ;
1327
1328 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0) ;
1329 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1330
1331// G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
1332
1333 x0 = fSplineEnergy[i - 1] ;
1334 x1 = fSplineEnergy[i - 2] ;
1335 y0 = fdNdxCerenkov[i - 1] ;
1336 yy1 = fdNdxCerenkov[i - 2] ;
1337
1338 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
1339 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1340
1341 c = x1/x0 ;
1342 d = e0/x0 ;
1343 a = log10(yy1/y0)/log10(x1/x0) ;
1344 // b0 = log10(y0) - a*log10(x0) ;
1345 b = y0/pow(x0,a); // pow(10.,b0) ;
1346
1347 a += 1.0 ;
1348 if( a == 0 ) result += b*log(e0/x0) ;
1349 else result += y0*(e0*pow(d,a-1) - x0 )/a ;
1350 a += 1.0 ;
1351
1352 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0) ;
1353 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1354
1355 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
1356 // <<b<<"; result = "<<result<<G4endl;
1357
1358 return result ;
1359
1360}
1361
1362///////////////////////////////////////////////////////////////////////////////
1363//
1364// Integration of Plasmon cross-section for the case of
1365// passing across border between intervals
1366
1367G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
1368 G4double en0 )
1369{
1370 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
1371
1372 e0 = en0 ;
1373 x0 = fSplineEnergy[i] ;
1374 x1 = fSplineEnergy[i+1] ;
1375 y0 = fdNdxPlasmon[i] ;
1376 yy1 = fdNdxPlasmon[i+1] ;
1377
1378 c = x1/x0 ;
1379 d = e0/x0 ;
1380 a = log10(yy1/y0)/log10(c) ;
1381 // b0 = log10(y0) - a*log10(x0) ;
1382 b = y0/pow(x0,a); //pow(10.,b) ;
1383
1384 a += 1.0 ;
1385 if( a == 0 ) result = b*log(x0/e0) ;
1386 else result = y0*(x0 - e0*pow(d,a-1))/a ;
1387 a += 1.0 ;
1388
1389 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0) ;
1390 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
1391
1392 x0 = fSplineEnergy[i - 1] ;
1393 x1 = fSplineEnergy[i - 2] ;
1394 y0 = fdNdxPlasmon[i - 1] ;
1395 yy1 = fdNdxPlasmon[i - 2] ;
1396
1397 c = x1/x0 ;
1398 d = e0/x0 ;
1399 a = log10(yy1/y0)/log10(c) ;
1400 // b0 = log10(y0) - a*log10(x0) ;
1401 b = y0/pow(x0,a);// pow(10.,b0) ;
1402
1403 a += 1.0 ;
1404 if( a == 0 ) result += b*log(e0/x0) ;
1405 else result += y0*(e0*pow(d,a-1) - x0)/a ;
1406 a += 1.0 ;
1407
1408 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0) ;
1409 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
1410
1411 return result ;
1412
1413}
1414
1415/////////////////////////////////////////////////////////////////////////
1416//
1417//
1418
1419G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
1420{
1421 G4int iTransfer ;
1422 G4long numOfCollisions ;
1423 G4double loss = 0.0 ;
1424 G4double meanNumber, position ;
1425
1426 // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl ;
1427
1428
1429
1430 meanNumber = fIntegralPAIxSection[1]*step ;
1431 numOfCollisions = G4Poisson(meanNumber) ;
1432
1433 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1434
1435 while(numOfCollisions)
1436 {
1437 position = fIntegralPAIxSection[1]*G4UniformRand() ;
1438
1439 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1440 {
1441 if( position >= fIntegralPAIxSection[iTransfer] ) break ;
1442 }
1443 loss += fSplineEnergy[iTransfer] ;
1444 numOfCollisions-- ;
1445 }
1446 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl ;
1447
1448 return loss ;
1449}
1450
1451/////////////////////////////////////////////////////////////////////////
1452//
1453//
1454
1455G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
1456{
1457 G4int iTransfer ;
1458 G4long numOfCollisions ;
1459 G4double loss = 0.0 ;
1460 G4double meanNumber, position ;
1461
1462 // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl ;
1463
1464
1465
1466 meanNumber = fIntegralCerenkov[1]*step ;
1467 numOfCollisions = G4Poisson(meanNumber) ;
1468
1469 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1470
1471 while(numOfCollisions)
1472 {
1473 position = fIntegralCerenkov[1]*G4UniformRand() ;
1474
1475 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1476 {
1477 if( position >= fIntegralCerenkov[iTransfer] ) break ;
1478 }
1479 loss += fSplineEnergy[iTransfer] ;
1480 numOfCollisions-- ;
1481 }
1482 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl ;
1483
1484 return loss ;
1485}
1486
1487/////////////////////////////////////////////////////////////////////////
1488//
1489//
1490
1491G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
1492{
1493 G4int iTransfer ;
1494 G4long numOfCollisions ;
1495 G4double loss = 0.0 ;
1496 G4double meanNumber, position ;
1497
1498 // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl ;
1499
1500
1501
1502 meanNumber = fIntegralPlasmon[1]*step ;
1503 numOfCollisions = G4Poisson(meanNumber) ;
1504
1505 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1506
1507 while(numOfCollisions)
1508 {
1509 position = fIntegralPlasmon[1]*G4UniformRand() ;
1510
1511 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ )
1512 {
1513 if( position >= fIntegralPlasmon[iTransfer] ) break ;
1514 }
1515 loss += fSplineEnergy[iTransfer] ;
1516 numOfCollisions-- ;
1517 }
1518 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl ;
1519
1520 return loss ;
1521}
1522
1523
1524
1525/////////////////////////////////////////////////////////////////////////////
1526//
1527// Init array of Lorentz factors
1528//
1529
1530G4int G4PAIxSection::fNumberOfGammas = 111 ;
1531
1532const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1
1533{
15340.0,
15351.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
15361.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10
15371.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
15381.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20
15392.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
15403.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30
15415.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
15428.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40
15431.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
15442.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50
15455.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
15461.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60
15471.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
15483.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70
15496.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
15501.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80
15512.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
15524.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90
15538.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
15541.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100
15553.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
15565.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110
15571.065799e+05
1558} ;
1559
1560///////////////////////////////////////////////////////////////////////
1561//
1562// The number of gamma for creation of spline (near ion-min , G ~ 4 )
1563//
1564
1565const
1566G4int G4PAIxSection::fRefGammaNumber = 29 ;
1567
1568
1569//
1570// end of G4PAIxSection implementation file
1571//
1572////////////////////////////////////////////////////////////////////////////
1573
Note: See TracBrowser for help on using the repository browser.