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