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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
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-cand-01 $
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.