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