00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 #include "G4PAIxSection.hh"
00052
00053 #include "globals.hh"
00054 #include "G4PhysicalConstants.hh"
00055 #include "G4SystemOfUnits.hh"
00056 #include "G4ios.hh"
00057 #include "G4Poisson.hh"
00058 #include "G4Material.hh"
00059 #include "G4MaterialCutsCouple.hh"
00060 #include "G4SandiaTable.hh"
00061
00062 using namespace std;
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 const G4double G4PAIxSection::fDelta = 0.005;
00085 const G4double G4PAIxSection::fError = 0.005;
00086
00087 const G4int G4PAIxSection::fMaxSplineSize = 500;
00088
00089
00091
00092
00093
00094
00095 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC)
00096 {
00097 fDensity = matCC->GetMaterial()->GetDensity();
00098 G4int matIndex = matCC->GetMaterial()->GetIndex();
00099 fMaterialIndex = matIndex;
00100 fSandia = new G4SandiaTable(matIndex);
00101
00102 G4int i, j;
00103 fMatSandiaMatrix = new G4OrderedTable();
00104
00105 for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
00106 {
00107 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
00108 }
00109 for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
00110 {
00111 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
00112
00113 for(j = 1; j < 5; j++)
00114 {
00115 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
00116 }
00117 }
00118 fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
00119 }
00120
00122
00123 G4PAIxSection::G4PAIxSection(G4int materialIndex,
00124 G4double maxEnergyTransfer)
00125 {
00126 fSandia = 0;
00127 fMatSandiaMatrix = 0;
00128 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00129 G4int i, j;
00130
00131 fMaterialIndex = materialIndex;
00132 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
00133 fElectronDensity = (*theMaterialTable)[materialIndex]->
00134 GetElectronDensity();
00135 fIntervalNumber = (*theMaterialTable)[materialIndex]->
00136 GetSandiaTable()->GetMatNbOfIntervals();
00137 fIntervalNumber--;
00138
00139
00140 fEnergyInterval = new G4double[fIntervalNumber+2];
00141 fA1 = new G4double[fIntervalNumber+2];
00142 fA2 = new G4double[fIntervalNumber+2];
00143 fA3 = new G4double[fIntervalNumber+2];
00144 fA4 = new G4double[fIntervalNumber+2];
00145
00146 for(i = 1; i <= fIntervalNumber; i++ )
00147 {
00148 if(((*theMaterialTable)[materialIndex]->
00149 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
00150 i > fIntervalNumber )
00151 {
00152 fEnergyInterval[i] = maxEnergyTransfer;
00153 fIntervalNumber = i;
00154 break;
00155 }
00156 fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
00157 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
00158 fA1[i] = (*theMaterialTable)[materialIndex]->
00159 GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
00160 fA2[i] = (*theMaterialTable)[materialIndex]->
00161 GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
00162 fA3[i] = (*theMaterialTable)[materialIndex]->
00163 GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
00164 fA4[i] = (*theMaterialTable)[materialIndex]->
00165 GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
00166
00167
00168 }
00169 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
00170 {
00171 fIntervalNumber++;
00172 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
00173 }
00174
00175
00176
00177 for(i=1;i<fIntervalNumber;i++)
00178 {
00179 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
00180 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
00181 {
00182 continue;
00183 }
00184 else
00185 {
00186 for(j=i;j<fIntervalNumber;j++)
00187 {
00188 fEnergyInterval[j] = fEnergyInterval[j+1];
00189 fA1[j] = fA1[j+1];
00190 fA2[j] = fA2[j+1];
00191 fA3[j] = fA3[j+1];
00192 fA4[j] = fA4[j+1];
00193 }
00194 fIntervalNumber--;
00195 i--;
00196 }
00197 }
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220 InitPAI();
00221
00222 delete[] fEnergyInterval;
00223 delete[] fA1;
00224 delete[] fA2;
00225 delete[] fA3;
00226 delete[] fA4;
00227 }
00228
00230
00231
00232
00233 G4PAIxSection::G4PAIxSection( G4int materialIndex,
00234 G4double maxEnergyTransfer,
00235 G4double betaGammaSq,
00236 G4double** photoAbsCof,
00237 G4int intNumber )
00238 {
00239 fSandia = 0;
00240 fMatSandiaMatrix = 0;
00241 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00242 G4int i, j;
00243
00244 fMaterialIndex = materialIndex;
00245 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
00246 fElectronDensity = (*theMaterialTable)[materialIndex]->
00247 GetElectronDensity();
00248
00249 fIntervalNumber = intNumber;
00250 fIntervalNumber--;
00251
00252
00253 fEnergyInterval = new G4double[fIntervalNumber+2];
00254 fA1 = new G4double[fIntervalNumber+2];
00255 fA2 = new G4double[fIntervalNumber+2];
00256 fA3 = new G4double[fIntervalNumber+2];
00257 fA4 = new G4double[fIntervalNumber+2];
00258
00259 for( i = 1; i <= fIntervalNumber; i++ )
00260 {
00261 if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
00262 i > fIntervalNumber )
00263 {
00264 fEnergyInterval[i] = maxEnergyTransfer;
00265 fIntervalNumber = i;
00266 break;
00267 }
00268 fEnergyInterval[i] = photoAbsCof[i-1][0];
00269 fA1[i] = photoAbsCof[i-1][1];
00270 fA2[i] = photoAbsCof[i-1][2];
00271 fA3[i] = photoAbsCof[i-1][3];
00272 fA4[i] = photoAbsCof[i-1][4];
00273
00274
00275 }
00276
00277 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
00278 {
00279 fIntervalNumber++;
00280 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
00281 }
00282 for(i=1;i<=fIntervalNumber;i++)
00283 {
00284
00285
00286 }
00287
00288
00289 for( i = 1; i < fIntervalNumber; i++ )
00290 {
00291 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
00292 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
00293 {
00294 continue;
00295 }
00296 else
00297 {
00298 for(j=i;j<fIntervalNumber;j++)
00299 {
00300 fEnergyInterval[j] = fEnergyInterval[j+1];
00301 fA1[j] = fA1[j+1];
00302 fA2[j] = fA2[j+1];
00303 fA3[j] = fA3[j+1];
00304 fA4[j] = fA4[j+1];
00305 }
00306 fIntervalNumber--;
00307 i--;
00308 }
00309 }
00310
00311
00312
00313 G4double betaGammaSqRef =
00314 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
00315
00316 NormShift(betaGammaSqRef);
00317 SplainPAI(betaGammaSqRef);
00318
00319
00320
00321 for(i = 1; i <= fSplineNumber; i++)
00322 {
00323 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00324 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
00325 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00326 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
00327 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
00328
00329
00330
00331 }
00332 IntegralCerenkov();
00333 IntegralMM();
00334 IntegralPlasmon();
00335 IntegralResonance();
00336 IntegralPAIxSection();
00337
00338 delete[] fEnergyInterval;
00339 delete[] fA1;
00340 delete[] fA2;
00341 delete[] fA3;
00342 delete[] fA4;
00343 }
00344
00346
00347
00348
00349 G4PAIxSection::G4PAIxSection( G4int materialIndex,
00350 G4double maxEnergyTransfer,
00351 G4double betaGammaSq )
00352 {
00353 fSandia = 0;
00354 fMatSandiaMatrix = 0;
00355 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00356
00357 G4int i, j, numberOfElements;
00358
00359 fMaterialIndex = materialIndex;
00360 fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
00361 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
00362 numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
00363
00364 G4int* thisMaterialZ = new G4int[numberOfElements];
00365
00366 for( i = 0; i < numberOfElements; i++ )
00367 {
00368 thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
00369 GetElement(i)->GetZ();
00370 }
00371
00372 fSandia = (*theMaterialTable)[materialIndex]->
00373 GetSandiaTable();
00374 G4SandiaTable thisMaterialSandiaTable(materialIndex);
00375 fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
00376 (thisMaterialZ,numberOfElements);
00377 fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
00378 ( thisMaterialZ ,
00379 (*theMaterialTable)[materialIndex]->GetFractionVector() ,
00380 numberOfElements,fIntervalNumber);
00381
00382 fIntervalNumber--;
00383
00384 fEnergyInterval = new G4double[fIntervalNumber+2];
00385 fA1 = new G4double[fIntervalNumber+2];
00386 fA2 = new G4double[fIntervalNumber+2];
00387 fA3 = new G4double[fIntervalNumber+2];
00388 fA4 = new G4double[fIntervalNumber+2];
00389
00390 for( i = 1; i <= fIntervalNumber; i++ )
00391 {
00392 if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
00393 i > fIntervalNumber)
00394 {
00395 fEnergyInterval[i] = maxEnergyTransfer;
00396 fIntervalNumber = i;
00397 break;
00398 }
00399 fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
00400 fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
00401 fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
00402 fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
00403 fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
00404
00405 }
00406 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
00407 {
00408 fIntervalNumber++;
00409 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
00410 fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
00411 fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
00412 fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
00413 fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
00414 }
00415 for(i=1;i<=fIntervalNumber;i++)
00416 {
00417
00418
00419 }
00420
00421
00422 for( i = 1; i < fIntervalNumber; i++ )
00423 {
00424 if(fEnergyInterval[i+1]-fEnergyInterval[i] >
00425 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
00426 {
00427 continue;
00428 }
00429 else
00430 {
00431 for( j = i; j < fIntervalNumber; j++ )
00432 {
00433 fEnergyInterval[j] = fEnergyInterval[j+1];
00434 fA1[j] = fA1[j+1];
00435 fA2[j] = fA2[j+1];
00436 fA3[j] = fA3[j+1];
00437 fA4[j] = fA4[j+1];
00438 }
00439 fIntervalNumber--;
00440 i--;
00441 }
00442 }
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 G4double betaGammaSqRef =
00466 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
00467
00468 NormShift(betaGammaSqRef);
00469 SplainPAI(betaGammaSqRef);
00470
00471
00472
00473 for(i = 1; i <= fSplineNumber; i++)
00474 {
00475 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
00476 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00477 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
00478 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00479 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
00480 }
00481 IntegralPAIxSection();
00482 IntegralCerenkov();
00483 IntegralMM();
00484 IntegralPlasmon();
00485 IntegralResonance();
00486
00487 delete[] fEnergyInterval;
00488 delete[] fA1;
00489 delete[] fA2;
00490 delete[] fA3;
00491 delete[] fA4;
00492 delete[] thisMaterialZ;
00493 }
00494
00495
00497
00498
00499
00500 G4PAIxSection::~G4PAIxSection()
00501 {
00502
00503
00504
00505
00506
00507
00508
00509
00510 delete fSandia;
00511 delete fMatSandiaMatrix;
00512 }
00513
00515
00516
00517
00518
00519 void G4PAIxSection::InitPAI()
00520 {
00521 G4int i;
00522 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
00523 fLorentzFactor[fRefGammaNumber] - 1;
00524
00525
00526
00527 NormShift(betaGammaSq);
00528 SplainPAI(betaGammaSq);
00529
00530 IntegralPAIxSection();
00531 IntegralCerenkov();
00532 IntegralMM();
00533 IntegralPlasmon();
00534 IntegralResonance();
00535
00536 for(i = 0; i<= fSplineNumber; i++)
00537 {
00538 fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
00539 if(i != 0)
00540 {
00541 fPAItable[i][0] = fSplineEnergy[i];
00542 }
00543 }
00544 fPAItable[0][0] = fSplineNumber;
00545
00546 for(G4int j = 1; j < 112; j++)
00547 {
00548 if( j == fRefGammaNumber ) continue;
00549
00550 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
00551
00552 for(i = 1; i <= fSplineNumber; i++)
00553 {
00554 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
00555 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00556 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
00557 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00558 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
00559 }
00560 IntegralPAIxSection();
00561 IntegralCerenkov();
00562 IntegralMM();
00563 IntegralPlasmon();
00564 IntegralResonance();
00565
00566 for(i = 0; i <= fSplineNumber; i++)
00567 {
00568 fPAItable[i][j] = fIntegralPAIxSection[i];
00569 }
00570 }
00571
00572 }
00573
00575
00576
00577
00578
00579 void G4PAIxSection::NormShift(G4double betaGammaSq)
00580 {
00581 G4int i, j;
00582
00583 for( i = 1; i <= fIntervalNumber-1; i++ )
00584 {
00585 for( j = 1; j <= 2; j++ )
00586 {
00587 fSplineNumber = (i-1)*2 + j;
00588
00589 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
00590 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
00591
00592
00593 }
00594 }
00595 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
00596
00597 j = 1;
00598
00599 for( i = 2; i <= fSplineNumber; i++ )
00600 {
00601 if(fSplineEnergy[i]<fEnergyInterval[j+1])
00602 {
00603 fIntegralTerm[i] = fIntegralTerm[i-1] +
00604 RutherfordIntegral(j,fSplineEnergy[i-1],
00605 fSplineEnergy[i] );
00606 }
00607 else
00608 {
00609 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
00610 fEnergyInterval[j+1] );
00611 j++;
00612 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
00613 RutherfordIntegral(j,fEnergyInterval[j],
00614 fSplineEnergy[i] );
00615 }
00616
00617 }
00618 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
00619 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
00620
00621
00622
00623
00624
00625
00626 for(G4int k = 1; k <= fIntervalNumber-1; k++ )
00627 {
00628 for( j = 1; j <= 2; j++ )
00629 {
00630 i = (k-1)*2 + j;
00631 fImPartDielectricConst[i] = fNormalizationCof*
00632 ImPartDielectricConst(k,fSplineEnergy[i]);
00633 fRePartDielectricConst[i] = fNormalizationCof*
00634 RePartDielectricConst(fSplineEnergy[i]);
00635 fIntegralTerm[i] *= fNormalizationCof;
00636
00637 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
00638 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00639 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
00640 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00641 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
00642 }
00643 }
00644
00645 }
00646
00648
00649
00650
00651
00652
00653 void G4PAIxSection::SplainPAI(G4double betaGammaSq)
00654 {
00655 G4int k = 1;
00656 G4int i = 1;
00657
00658 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
00659 {
00660 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
00661 {
00662 k++;
00663 i++;
00664 continue;
00665 }
00666
00667
00668 fSplineNumber++;
00669
00670 for(G4int j = fSplineNumber; j >= i+2; j-- )
00671 {
00672 fSplineEnergy[j] = fSplineEnergy[j-1];
00673 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
00674 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
00675 fIntegralTerm[j] = fIntegralTerm[j-1];
00676
00677 fDifPAIxSection[j] = fDifPAIxSection[j-1];
00678 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
00679 fdNdxMM[j] = fdNdxMM[j-1];
00680 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
00681 fdNdxResonance[j] = fdNdxResonance[j-1];
00682 }
00683 G4double x1 = fSplineEnergy[i];
00684 G4double x2 = fSplineEnergy[i+1];
00685 G4double yy1 = fDifPAIxSection[i];
00686 G4double y2 = fDifPAIxSection[i+1];
00687
00688 G4double en1 = sqrt(x1*x2);
00689 fSplineEnergy[i+1] = en1;
00690
00691
00692
00693
00694 G4double a = log10(y2/yy1)/log10(x2/x1);
00695 G4double b = log10(yy1) - a*log10(x1);
00696 G4double y = a*log10(en1) + b;
00697 y = pow(10.,y);
00698
00699
00700
00701 fImPartDielectricConst[i+1] = fNormalizationCof*
00702 ImPartDielectricConst(k,fSplineEnergy[i+1]);
00703 fRePartDielectricConst[i+1] = fNormalizationCof*
00704 RePartDielectricConst(fSplineEnergy[i+1]);
00705 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
00706 RutherfordIntegral(k,fSplineEnergy[i],
00707 fSplineEnergy[i+1]);
00708
00709 fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
00710 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
00711 fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
00712 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
00713 fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
00714
00715
00716
00717
00718 G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
00719
00720 if( x < 0 )
00721 {
00722 x = -x;
00723 }
00724 if( x > fError && fSplineNumber < fMaxSplineSize-1 )
00725 {
00726 continue;
00727 }
00728 i += 2;
00729
00730 }
00731
00732 }
00733
00734
00736
00737
00738
00739
00740 G4double G4PAIxSection::RutherfordIntegral( G4int k,
00741 G4double x1,
00742 G4double x2 )
00743 {
00744 G4double c1, c2, c3;
00745
00746 c1 = (x2 - x1)/x1/x2;
00747 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
00748 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
00749
00750
00751 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
00752
00753 }
00754
00755
00757
00758
00759
00760
00761 G4double G4PAIxSection::ImPartDielectricConst( G4int k ,
00762 G4double energy1 )
00763 {
00764 G4double energy2,energy3,energy4,result;
00765
00766 energy2 = energy1*energy1;
00767 energy3 = energy2*energy1;
00768 energy4 = energy3*energy1;
00769
00770 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
00771 result *=hbarc/energy1;
00772
00773 return result;
00774
00775 }
00776
00778
00779
00780
00781 G4double G4PAIxSection::GetPhotonRange( G4double energy1 )
00782 {
00783 G4int i;
00784 G4double energy2, energy3, energy4, result, lambda;
00785
00786 energy2 = energy1*energy1;
00787 energy3 = energy2*energy1;
00788 energy4 = energy3*energy1;
00789
00790
00791
00792
00793
00794 for( i = 1; i <= fIntervalNumber; i++ )
00795 {
00796 if( energy1 < fEnergyInterval[i]) break;
00797 }
00798 i--;
00799 if(i == 0) i = 1;
00800
00801 result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
00802
00803 if( result > DBL_MIN ) lambda = 1./result;
00804 else lambda = DBL_MAX;
00805
00806 return lambda;
00807 }
00808
00810
00811
00812
00813
00814 G4double G4PAIxSection::GetElectronRange( G4double energy )
00815 {
00816 G4double range;
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837 G4double cofA = 5.37e-4*g/cm2/keV;
00838 G4double cofB = 0.9815;
00839 G4double cofC = 3.123e-3/keV;
00840
00841
00842 range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
00843
00844
00845 range /= fDensity;
00846
00847 return range;
00848 }
00849
00851
00852
00853
00854
00855
00856 G4double G4PAIxSection::RePartDielectricConst(G4double enb)
00857 {
00858 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
00859 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
00860
00861 x0 = enb;
00862 result = 0;
00863
00864 for(G4int i=1;i<=fIntervalNumber-1;i++)
00865 {
00866 x1 = fEnergyInterval[i];
00867 x2 = fEnergyInterval[i+1];
00868 xx1 = x1 - x0;
00869 xx2 = x2 - x0;
00870 xx12 = xx2/xx1;
00871
00872 if(xx12<0)
00873 {
00874 xx12 = -xx12;
00875 }
00876 xln1 = log(x2/x1);
00877 xln2 = log(xx12);
00878 xln3 = log((x2 + x0)/(x1 + x0));
00879 x02 = x0*x0;
00880 x03 = x02*x0;
00881 x04 = x03*x0;
00882 x05 = x04*x0;
00883 c1 = (x2 - x1)/x1/x2;
00884 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
00885 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
00886
00887 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
00888 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
00889 result -= fA3[i]*c2/2/x02;
00890 result -= fA4[i]*c3/3/x02;
00891
00892 cof1 = fA1[i]/x02 + fA3[i]/x04;
00893 cof2 = fA2[i]/x03 + fA4[i]/x05;
00894
00895 result += 0.5*(cof1 +cof2)*xln2;
00896 result += 0.5*(cof1 - cof2)*xln3;
00897 }
00898 result *= 2*hbarc/pi;
00899
00900 return result;
00901
00902 }
00903
00905
00906
00907
00908
00909
00910 G4double G4PAIxSection::DifPAIxSection( G4int i ,
00911 G4double betaGammaSq )
00912 {
00913 G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
00914
00915 G4double betaBohr = fine_structure_const;
00916
00917
00918
00919 G4double be2 = betaGammaSq/(1 + betaGammaSq);
00920 G4double beta = sqrt(be2);
00921
00922
00923 cof = 1.;
00924 x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
00925
00926 if( betaGammaSq < 0.01 ) x2 = log(be2);
00927 else
00928 {
00929 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
00930 (1/betaGammaSq - fRePartDielectricConst[i]) +
00931 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
00932 }
00933 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
00934 {
00935 x6 = 0.;
00936 }
00937 else
00938 {
00939 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
00940 x5 = -1 - fRePartDielectricConst[i] +
00941 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
00942 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
00943
00944 x7 = atan2(fImPartDielectricConst[i],x3);
00945 x6 = x5 * x7;
00946 }
00947
00948
00949 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
00950
00951
00952
00953 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
00954 fImPartDielectricConst[i]*fImPartDielectricConst[i];
00955
00956 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
00957
00958 if( result < 1.0e-8 ) result = 1.0e-8;
00959
00960 result *= fine_structure_const/be2/pi;
00961
00962
00963
00964 G4double lowCof = 4.0 ;
00965
00966 result *= (1 - exp(-beta/betaBohr/lowCof));
00967
00968
00969
00970
00971
00972
00973
00974 if(fDensity >= 0.1)
00975 {
00976 result /= x8;
00977 }
00978 return result;
00979
00980 }
00981
00983
00984
00985
00986 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i ,
00987 G4double betaGammaSq )
00988 {
00989 G4double logarithm, x3, x5, argument, modul2, dNdxC;
00990 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
00991
00992 cofBetaBohr = 4.0;
00993 betaBohr2 = fine_structure_const*fine_structure_const;
00994 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
00995
00996 be2 = betaGammaSq/(1 + betaGammaSq);
00997 be4 = be2*be2;
00998
00999 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq);
01000 else
01001 {
01002 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
01003 (1/betaGammaSq - fRePartDielectricConst[i]) +
01004 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
01005 logarithm += log(1+1.0/betaGammaSq);
01006 }
01007
01008 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
01009 {
01010 argument = 0.0;
01011 }
01012 else
01013 {
01014 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
01015 x5 = -1.0 - fRePartDielectricConst[i] +
01016 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
01017 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
01018 if( x3 == 0.0 ) argument = 0.5*pi;
01019 else argument = atan2(fImPartDielectricConst[i],x3);
01020 argument *= x5 ;
01021 }
01022 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
01023
01024 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
01025
01026 dNdxC *= fine_structure_const/be2/pi;
01027
01028 dNdxC *= (1-exp(-be4/betaBohr4));
01029
01030 if(fDensity >= 0.1)
01031 {
01032 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
01033 fImPartDielectricConst[i]*fImPartDielectricConst[i];
01034 dNdxC /= modul2;
01035 }
01036 return dNdxC;
01037
01038 }
01039
01041
01042
01043
01044 G4double G4PAIxSection::PAIdNdxMM( G4int i ,
01045 G4double betaGammaSq )
01046 {
01047 G4double logarithm, x3, x5, argument, dNdxC;
01048 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
01049
01050 cofBetaBohr = 4.0;
01051 betaBohr2 = fine_structure_const*fine_structure_const;
01052 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
01053
01054 be2 = betaGammaSq/(1 + betaGammaSq);
01055 be4 = be2*be2;
01056
01057 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq);
01058 else
01059 {
01060 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
01061 (1/betaGammaSq - fRePartDielectricConst[i]) +
01062 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
01063 logarithm += log(1+1.0/betaGammaSq);
01064 }
01065
01066 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
01067 {
01068 argument = 0.0;
01069 }
01070 else
01071 {
01072 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
01073 x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
01074 if( x3 == 0.0 ) argument = 0.5*pi;
01075 else argument = atan2(fImPartDielectricConst[i],x3);
01076 argument *= x5 ;
01077 }
01078 dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
01079
01080 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
01081
01082 dNdxC *= fine_structure_const/be2/pi;
01083
01084 dNdxC *= (1-exp(-be4/betaBohr4));
01085 return dNdxC;
01086
01087 }
01088
01090
01091
01092
01093
01094 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i ,
01095 G4double betaGammaSq )
01096 {
01097 G4double resonance, modul2, dNdxP, cof = 1.;
01098 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
01099
01100
01101 cofBetaBohr = 4.0;
01102 betaBohr2 = fine_structure_const*fine_structure_const;
01103 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
01104
01105 be2 = betaGammaSq/(1 + betaGammaSq);
01106 be4 = be2*be2;
01107
01108 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
01109 resonance *= fImPartDielectricConst[i]/hbarc;
01110
01111
01112 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
01113
01114 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
01115
01116 dNdxP *= fine_structure_const/be2/pi;
01117 dNdxP *= (1-exp(-be4/betaBohr4));
01118
01119 if( fDensity >= 0.1 )
01120 {
01121 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
01122 fImPartDielectricConst[i]*fImPartDielectricConst[i];
01123 dNdxP /= modul2;
01124 }
01125 return dNdxP;
01126
01127 }
01128
01130
01131
01132
01133
01134 G4double G4PAIxSection::PAIdNdxResonance( G4int i ,
01135 G4double betaGammaSq )
01136 {
01137 G4double resonance, modul2, dNdxP;
01138 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
01139
01140 cofBetaBohr = 4.0;
01141 betaBohr2 = fine_structure_const*fine_structure_const;
01142 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
01143
01144 be2 = betaGammaSq/(1 + betaGammaSq);
01145 be4 = be2*be2;
01146
01147 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
01148 resonance *= fImPartDielectricConst[i]/hbarc;
01149
01150
01151 dNdxP = resonance;
01152
01153 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
01154
01155 dNdxP *= fine_structure_const/be2/pi;
01156 dNdxP *= (1-exp(-be4/betaBohr4));
01157
01158 if( fDensity >= 0.1 )
01159 {
01160 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
01161 fImPartDielectricConst[i]*fImPartDielectricConst[i];
01162 dNdxP /= modul2;
01163 }
01164 return dNdxP;
01165
01166 }
01167
01169
01170
01171
01172
01173
01174 void G4PAIxSection::IntegralPAIxSection()
01175 {
01176 fIntegralPAIxSection[fSplineNumber] = 0;
01177 fIntegralPAIdEdx[fSplineNumber] = 0;
01178 fIntegralPAIxSection[0] = 0;
01179 G4int k = fIntervalNumber -1;
01180
01181 for(G4int i = fSplineNumber-1; i >= 1; i--)
01182 {
01183 if(fSplineEnergy[i] >= fEnergyInterval[k])
01184 {
01185 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
01186 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
01187 }
01188 else
01189 {
01190 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
01191 SumOverBorder(i+1,fEnergyInterval[k]);
01192 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
01193 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
01194 k--;
01195 }
01196 }
01197 }
01198
01200
01201
01202
01203
01204
01205 void G4PAIxSection::IntegralCerenkov()
01206 {
01207 G4int i, k;
01208 fIntegralCerenkov[fSplineNumber] = 0;
01209 fIntegralCerenkov[0] = 0;
01210 k = fIntervalNumber -1;
01211
01212 for( i = fSplineNumber-1; i >= 1; i-- )
01213 {
01214 if(fSplineEnergy[i] >= fEnergyInterval[k])
01215 {
01216 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
01217
01218 }
01219 else
01220 {
01221 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
01222 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
01223 k--;
01224
01225 }
01226 }
01227
01228 }
01229
01231
01232
01233
01234
01235
01236 void G4PAIxSection::IntegralMM()
01237 {
01238 G4int i, k;
01239 fIntegralMM[fSplineNumber] = 0;
01240 fIntegralMM[0] = 0;
01241 k = fIntervalNumber -1;
01242
01243 for( i = fSplineNumber-1; i >= 1; i-- )
01244 {
01245 if(fSplineEnergy[i] >= fEnergyInterval[k])
01246 {
01247 fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
01248
01249 }
01250 else
01251 {
01252 fIntegralMM[i] = fIntegralMM[i+1] +
01253 SumOverBordMM(i+1,fEnergyInterval[k]);
01254 k--;
01255
01256 }
01257 }
01258
01259 }
01260
01262
01263
01264
01265
01266
01267 void G4PAIxSection::IntegralPlasmon()
01268 {
01269 fIntegralPlasmon[fSplineNumber] = 0;
01270 fIntegralPlasmon[0] = 0;
01271 G4int k = fIntervalNumber -1;
01272 for(G4int i=fSplineNumber-1;i>=1;i--)
01273 {
01274 if(fSplineEnergy[i] >= fEnergyInterval[k])
01275 {
01276 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
01277 }
01278 else
01279 {
01280 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
01281 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
01282 k--;
01283 }
01284 }
01285
01286 }
01287
01289
01290
01291
01292
01293
01294 void G4PAIxSection::IntegralResonance()
01295 {
01296 fIntegralResonance[fSplineNumber] = 0;
01297 fIntegralResonance[0] = 0;
01298 G4int k = fIntervalNumber -1;
01299 for(G4int i=fSplineNumber-1;i>=1;i--)
01300 {
01301 if(fSplineEnergy[i] >= fEnergyInterval[k])
01302 {
01303 fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
01304 }
01305 else
01306 {
01307 fIntegralResonance[i] = fIntegralResonance[i+1] +
01308 SumOverBordResonance(i+1,fEnergyInterval[k]);
01309 k--;
01310 }
01311 }
01312
01313 }
01314
01316
01317
01318
01319
01320
01321 G4double G4PAIxSection::SumOverInterval( G4int i )
01322 {
01323 G4double x0,x1,y0,yy1,a,b,c,result;
01324
01325 x0 = fSplineEnergy[i];
01326 x1 = fSplineEnergy[i+1];
01327 y0 = fDifPAIxSection[i];
01328 yy1 = fDifPAIxSection[i+1];
01329 c = x1/x0;
01330 a = log10(yy1/y0)/log10(c);
01331
01332 b = y0/pow(x0,a);
01333 a += 1.;
01334 if( std::fabs(a) < 1.e-6 )
01335 {
01336 result = b*log(x1/x0);
01337 }
01338 else
01339 {
01340 result = y0*(x1*pow(c,a-1) - x0)/a;
01341 }
01342 a += 1.;
01343 if( std::fabs(a) < 1.e-6 )
01344 {
01345 fIntegralPAIxSection[0] += b*log(x1/x0);
01346 }
01347 else
01348 {
01349 fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01350 }
01351 return result;
01352
01353 }
01354
01356
01357 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i )
01358 {
01359 G4double x0,x1,y0,yy1,a,b,c,result;
01360
01361 x0 = fSplineEnergy[i];
01362 x1 = fSplineEnergy[i+1];
01363 y0 = fDifPAIxSection[i];
01364 yy1 = fDifPAIxSection[i+1];
01365 c = x1/x0;
01366 a = log10(yy1/y0)/log10(c);
01367
01368 b = y0/pow(x0,a);
01369 a += 2;
01370 if(a == 0)
01371 {
01372 result = b*log(x1/x0);
01373 }
01374 else
01375 {
01376 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01377 }
01378 return result;
01379
01380 }
01381
01383
01384
01385
01386
01387
01388 G4double G4PAIxSection::SumOverInterCerenkov( G4int i )
01389 {
01390 G4double x0,x1,y0,yy1,a,b,c,result;
01391
01392 x0 = fSplineEnergy[i];
01393 x1 = fSplineEnergy[i+1];
01394 y0 = fdNdxCerenkov[i];
01395 yy1 = fdNdxCerenkov[i+1];
01396
01397
01398
01399 c = x1/x0;
01400 a = log10(yy1/y0)/log10(c);
01401 b = y0/pow(x0,a);
01402
01403 a += 1.0;
01404 if(a == 0) result = b*log(c);
01405 else result = y0*(x1*pow(c,a-1) - x0)/a;
01406 a += 1.0;
01407
01408 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
01409 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01410
01411 return result;
01412
01413 }
01414
01416
01417
01418
01419
01420
01421 G4double G4PAIxSection::SumOverInterMM( G4int i )
01422 {
01423 G4double x0,x1,y0,yy1,a,b,c,result;
01424
01425 x0 = fSplineEnergy[i];
01426 x1 = fSplineEnergy[i+1];
01427 y0 = fdNdxMM[i];
01428 yy1 = fdNdxMM[i+1];
01429
01430
01431
01432 c = x1/x0;
01433 a = log10(yy1/y0)/log10(c);
01434 b = y0/pow(x0,a);
01435
01436 a += 1.0;
01437 if(a == 0) result = b*log(c);
01438 else result = y0*(x1*pow(c,a-1) - x0)/a;
01439 a += 1.0;
01440
01441 if( a == 0 ) fIntegralMM[0] += b*log(x1/x0);
01442 else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01443
01444 return result;
01445
01446 }
01447
01449
01450
01451
01452
01453
01454 G4double G4PAIxSection::SumOverInterPlasmon( G4int i )
01455 {
01456 G4double x0,x1,y0,yy1,a,b,c,result;
01457
01458 x0 = fSplineEnergy[i];
01459 x1 = fSplineEnergy[i+1];
01460 y0 = fdNdxPlasmon[i];
01461 yy1 = fdNdxPlasmon[i+1];
01462 c =x1/x0;
01463 a = log10(yy1/y0)/log10(c);
01464
01465 b = y0/pow(x0,a);
01466
01467 a += 1.0;
01468 if(a == 0) result = b*log(x1/x0);
01469 else result = y0*(x1*pow(c,a-1) - x0)/a;
01470 a += 1.0;
01471
01472 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
01473 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01474
01475 return result;
01476
01477 }
01478
01480
01481
01482
01483
01484
01485 G4double G4PAIxSection::SumOverInterResonance( G4int i )
01486 {
01487 G4double x0,x1,y0,yy1,a,b,c,result;
01488
01489 x0 = fSplineEnergy[i];
01490 x1 = fSplineEnergy[i+1];
01491 y0 = fdNdxResonance[i];
01492 yy1 = fdNdxResonance[i+1];
01493 c =x1/x0;
01494 a = log10(yy1/y0)/log10(c);
01495
01496 b = y0/pow(x0,a);
01497
01498 a += 1.0;
01499 if(a == 0) result = b*log(x1/x0);
01500 else result = y0*(x1*pow(c,a-1) - x0)/a;
01501 a += 1.0;
01502
01503 if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
01504 else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
01505
01506 return result;
01507
01508 }
01509
01511
01512
01513
01514
01515 G4double G4PAIxSection::SumOverBorder( G4int i ,
01516 G4double en0 )
01517 {
01518 G4double x0,x1,y0,yy1,a,b,d,e0,result;
01519
01520 e0 = en0;
01521 x0 = fSplineEnergy[i];
01522 x1 = fSplineEnergy[i+1];
01523 y0 = fDifPAIxSection[i];
01524 yy1 = fDifPAIxSection[i+1];
01525
01526
01527 d = e0/x0;
01528 a = log10(yy1/y0)/log10(x1/x0);
01529
01530 b = y0/pow(x0,a);
01531
01532 a += 1.;
01533 if( std::fabs(a) < 1.e-6 )
01534 {
01535 result = b*log(x0/e0);
01536 }
01537 else
01538 {
01539 result = y0*(x0 - e0*pow(d,a-1))/a;
01540 }
01541 a += 1.;
01542 if( std::fabs(a) < 1.e-6 )
01543 {
01544 fIntegralPAIxSection[0] += b*log(x0/e0);
01545 }
01546 else
01547 {
01548 fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01549 }
01550 x0 = fSplineEnergy[i - 1];
01551 x1 = fSplineEnergy[i - 2];
01552 y0 = fDifPAIxSection[i - 1];
01553 yy1 = fDifPAIxSection[i - 2];
01554
01555
01556 d = e0/x0;
01557 a = log10(yy1/y0)/log10(x1/x0);
01558
01559 b = y0/pow(x0,a);
01560 a += 1.;
01561 if( std::fabs(a) < 1.e-6 )
01562 {
01563 result += b*log(e0/x0);
01564 }
01565 else
01566 {
01567 result += y0*(e0*pow(d,a-1) - x0)/a;
01568 }
01569 a += 1.;
01570 if( std::fabs(a) < 1.e-6 )
01571 {
01572 fIntegralPAIxSection[0] += b*log(e0/x0);
01573 }
01574 else
01575 {
01576 fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01577 }
01578 return result;
01579
01580 }
01581
01583
01584 G4double G4PAIxSection::SumOverBorderdEdx( G4int i ,
01585 G4double en0 )
01586 {
01587 G4double x0,x1,y0,yy1,a,b,d,e0,result;
01588
01589 e0 = en0;
01590 x0 = fSplineEnergy[i];
01591 x1 = fSplineEnergy[i+1];
01592 y0 = fDifPAIxSection[i];
01593 yy1 = fDifPAIxSection[i+1];
01594
01595
01596 d = e0/x0;
01597 a = log10(yy1/y0)/log10(x1/x0);
01598
01599 b = y0/pow(x0,a);
01600
01601 a += 2;
01602 if(a == 0)
01603 {
01604 result = b*log(x0/e0);
01605 }
01606 else
01607 {
01608 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01609 }
01610 x0 = fSplineEnergy[i - 1];
01611 x1 = fSplineEnergy[i - 2];
01612 y0 = fDifPAIxSection[i - 1];
01613 yy1 = fDifPAIxSection[i - 2];
01614
01615
01616 d = e0/x0;
01617 a = log10(yy1/y0)/log10(x1/x0);
01618
01619 b = y0/pow(x0,a);
01620 a += 2;
01621 if(a == 0)
01622 {
01623 result += b*log(e0/x0);
01624 }
01625 else
01626 {
01627 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01628 }
01629 return result;
01630
01631 }
01632
01634
01635
01636
01637
01638 G4double G4PAIxSection::SumOverBordCerenkov( G4int i ,
01639 G4double en0 )
01640 {
01641 G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
01642
01643 e0 = en0;
01644 x0 = fSplineEnergy[i];
01645 x1 = fSplineEnergy[i+1];
01646 y0 = fdNdxCerenkov[i];
01647 yy1 = fdNdxCerenkov[i+1];
01648
01649
01650
01651
01652 c = x1/x0;
01653 d = e0/x0;
01654 a = log10(yy1/y0)/log10(c);
01655
01656 b = y0/pow(x0,a);
01657
01658 a += 1.0;
01659 if( a == 0 ) result = b*log(x0/e0);
01660 else result = y0*(x0 - e0*pow(d,a-1))/a;
01661 a += 1.0;
01662
01663 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
01664 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01665
01666
01667
01668 x0 = fSplineEnergy[i - 1];
01669 x1 = fSplineEnergy[i - 2];
01670 y0 = fdNdxCerenkov[i - 1];
01671 yy1 = fdNdxCerenkov[i - 2];
01672
01673
01674
01675
01676 c = x1/x0;
01677 d = e0/x0;
01678 a = log10(yy1/y0)/log10(x1/x0);
01679
01680 b = y0/pow(x0,a);
01681
01682 a += 1.0;
01683 if( a == 0 ) result += b*log(e0/x0);
01684 else result += y0*(e0*pow(d,a-1) - x0 )/a;
01685 a += 1.0;
01686
01687 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
01688 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01689
01690
01691
01692
01693 return result;
01694
01695 }
01696
01698
01699
01700
01701
01702 G4double G4PAIxSection::SumOverBordMM( G4int i ,
01703 G4double en0 )
01704 {
01705 G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
01706
01707 e0 = en0;
01708 x0 = fSplineEnergy[i];
01709 x1 = fSplineEnergy[i+1];
01710 y0 = fdNdxMM[i];
01711 yy1 = fdNdxMM[i+1];
01712
01713
01714
01715
01716 c = x1/x0;
01717 d = e0/x0;
01718 a = log10(yy1/y0)/log10(c);
01719
01720 b = y0/pow(x0,a);
01721
01722 a += 1.0;
01723 if( a == 0 ) result = b*log(x0/e0);
01724 else result = y0*(x0 - e0*pow(d,a-1))/a;
01725 a += 1.0;
01726
01727 if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
01728 else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01729
01730
01731
01732 x0 = fSplineEnergy[i - 1];
01733 x1 = fSplineEnergy[i - 2];
01734 y0 = fdNdxMM[i - 1];
01735 yy1 = fdNdxMM[i - 2];
01736
01737
01738
01739
01740 c = x1/x0;
01741 d = e0/x0;
01742 a = log10(yy1/y0)/log10(x1/x0);
01743
01744 b = y0/pow(x0,a);
01745
01746 a += 1.0;
01747 if( a == 0 ) result += b*log(e0/x0);
01748 else result += y0*(e0*pow(d,a-1) - x0 )/a;
01749 a += 1.0;
01750
01751 if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
01752 else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01753
01754
01755
01756
01757 return result;
01758
01759 }
01760
01762
01763
01764
01765
01766 G4double G4PAIxSection::SumOverBordPlasmon( G4int i ,
01767 G4double en0 )
01768 {
01769 G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
01770
01771 e0 = en0;
01772 x0 = fSplineEnergy[i];
01773 x1 = fSplineEnergy[i+1];
01774 y0 = fdNdxPlasmon[i];
01775 yy1 = fdNdxPlasmon[i+1];
01776
01777 c = x1/x0;
01778 d = e0/x0;
01779 a = log10(yy1/y0)/log10(c);
01780
01781 b = y0/pow(x0,a);
01782
01783 a += 1.0;
01784 if( a == 0 ) result = b*log(x0/e0);
01785 else result = y0*(x0 - e0*pow(d,a-1))/a;
01786 a += 1.0;
01787
01788 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
01789 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01790
01791 x0 = fSplineEnergy[i - 1];
01792 x1 = fSplineEnergy[i - 2];
01793 y0 = fdNdxPlasmon[i - 1];
01794 yy1 = fdNdxPlasmon[i - 2];
01795
01796 c = x1/x0;
01797 d = e0/x0;
01798 a = log10(yy1/y0)/log10(c);
01799
01800 b = y0/pow(x0,a);
01801
01802 a += 1.0;
01803 if( a == 0 ) result += b*log(e0/x0);
01804 else result += y0*(e0*pow(d,a-1) - x0)/a;
01805 a += 1.0;
01806
01807 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
01808 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01809
01810 return result;
01811
01812 }
01813
01815
01816
01817
01818
01819 G4double G4PAIxSection::SumOverBordResonance( G4int i ,
01820 G4double en0 )
01821 {
01822 G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
01823
01824 e0 = en0;
01825 x0 = fSplineEnergy[i];
01826 x1 = fSplineEnergy[i+1];
01827 y0 = fdNdxResonance[i];
01828 yy1 = fdNdxResonance[i+1];
01829
01830 c = x1/x0;
01831 d = e0/x0;
01832 a = log10(yy1/y0)/log10(c);
01833
01834 b = y0/pow(x0,a);
01835
01836 a += 1.0;
01837 if( a == 0 ) result = b*log(x0/e0);
01838 else result = y0*(x0 - e0*pow(d,a-1))/a;
01839 a += 1.0;
01840
01841 if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
01842 else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01843
01844 x0 = fSplineEnergy[i - 1];
01845 x1 = fSplineEnergy[i - 2];
01846 y0 = fdNdxResonance[i - 1];
01847 yy1 = fdNdxResonance[i - 2];
01848
01849 c = x1/x0;
01850 d = e0/x0;
01851 a = log10(yy1/y0)/log10(c);
01852
01853 b = y0/pow(x0,a);
01854
01855 a += 1.0;
01856 if( a == 0 ) result += b*log(e0/x0);
01857 else result += y0*(e0*pow(d,a-1) - x0)/a;
01858 a += 1.0;
01859
01860 if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
01861 else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01862
01863 return result;
01864
01865 }
01866
01868
01869
01870
01871 G4double G4PAIxSection::GetStepEnergyLoss( G4double step )
01872 {
01873 G4long numOfCollisions;
01874 G4double meanNumber, loss = 0.0;
01875
01876
01877
01878 meanNumber = fIntegralPAIxSection[1]*step;
01879 numOfCollisions = G4Poisson(meanNumber);
01880
01881
01882
01883 while(numOfCollisions)
01884 {
01885 loss += GetEnergyTransfer();
01886 numOfCollisions--;
01887 }
01888
01889
01890 return loss;
01891 }
01892
01894
01895
01896
01897 G4double G4PAIxSection::GetEnergyTransfer()
01898 {
01899 G4int iTransfer ;
01900
01901 G4double energyTransfer, position;
01902
01903 position = fIntegralPAIxSection[1]*G4UniformRand();
01904
01905 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
01906 {
01907 if( position >= fIntegralPAIxSection[iTransfer] ) break;
01908 }
01909 if(iTransfer > fSplineNumber) iTransfer--;
01910
01911 energyTransfer = fSplineEnergy[iTransfer];
01912
01913 if(iTransfer > 1)
01914 {
01915 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
01916 }
01917 return energyTransfer;
01918 }
01919
01921
01922
01923
01924 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step )
01925 {
01926 G4long numOfCollisions;
01927 G4double meanNumber, loss = 0.0;
01928
01929
01930
01931 meanNumber = fIntegralCerenkov[1]*step;
01932 numOfCollisions = G4Poisson(meanNumber);
01933
01934
01935
01936 while(numOfCollisions)
01937 {
01938 loss += GetCerenkovEnergyTransfer();
01939 numOfCollisions--;
01940 }
01941
01942
01943 return loss;
01944 }
01945
01947
01948
01949
01950 G4double G4PAIxSection::GetStepMMLoss( G4double step )
01951 {
01952 G4long numOfCollisions;
01953 G4double meanNumber, loss = 0.0;
01954
01955
01956
01957 meanNumber = fIntegralMM[1]*step;
01958 numOfCollisions = G4Poisson(meanNumber);
01959
01960
01961
01962 while(numOfCollisions)
01963 {
01964 loss += GetMMEnergyTransfer();
01965 numOfCollisions--;
01966 }
01967
01968
01969 return loss;
01970 }
01971
01973
01974
01975
01976 G4double G4PAIxSection::GetCerenkovEnergyTransfer()
01977 {
01978 G4int iTransfer ;
01979
01980 G4double energyTransfer, position;
01981
01982 position = fIntegralCerenkov[1]*G4UniformRand();
01983
01984 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
01985 {
01986 if( position >= fIntegralCerenkov[iTransfer] ) break;
01987 }
01988 if(iTransfer > fSplineNumber) iTransfer--;
01989
01990 energyTransfer = fSplineEnergy[iTransfer];
01991
01992 if(iTransfer > 1)
01993 {
01994 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
01995 }
01996 return energyTransfer;
01997 }
01998
02000
02001
02002
02003 G4double G4PAIxSection::GetMMEnergyTransfer()
02004 {
02005 G4int iTransfer ;
02006
02007 G4double energyTransfer, position;
02008
02009 position = fIntegralMM[1]*G4UniformRand();
02010
02011 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
02012 {
02013 if( position >= fIntegralMM[iTransfer] ) break;
02014 }
02015 if(iTransfer > fSplineNumber) iTransfer--;
02016
02017 energyTransfer = fSplineEnergy[iTransfer];
02018
02019 if(iTransfer > 1)
02020 {
02021 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
02022 }
02023 return energyTransfer;
02024 }
02025
02027
02028
02029
02030 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step )
02031 {
02032 G4long numOfCollisions;
02033 G4double meanNumber, loss = 0.0;
02034
02035
02036
02037 meanNumber = fIntegralPlasmon[1]*step;
02038 numOfCollisions = G4Poisson(meanNumber);
02039
02040
02041
02042 while(numOfCollisions)
02043 {
02044 loss += GetPlasmonEnergyTransfer();
02045 numOfCollisions--;
02046 }
02047
02048
02049 return loss;
02050 }
02051
02053
02054
02055
02056 G4double G4PAIxSection::GetPlasmonEnergyTransfer()
02057 {
02058 G4int iTransfer ;
02059
02060 G4double energyTransfer, position;
02061
02062 position = fIntegralPlasmon[1]*G4UniformRand();
02063
02064 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
02065 {
02066 if( position >= fIntegralPlasmon[iTransfer] ) break;
02067 }
02068 if(iTransfer > fSplineNumber) iTransfer--;
02069
02070 energyTransfer = fSplineEnergy[iTransfer];
02071
02072 if(iTransfer > 1)
02073 {
02074 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
02075 }
02076 return energyTransfer;
02077 }
02078
02080
02081
02082
02083 G4double G4PAIxSection::GetStepResonanceLoss( G4double step )
02084 {
02085 G4long numOfCollisions;
02086 G4double meanNumber, loss = 0.0;
02087
02088
02089
02090 meanNumber = fIntegralResonance[1]*step;
02091 numOfCollisions = G4Poisson(meanNumber);
02092
02093
02094
02095 while(numOfCollisions)
02096 {
02097 loss += GetResonanceEnergyTransfer();
02098 numOfCollisions--;
02099 }
02100
02101
02102 return loss;
02103 }
02104
02105
02107
02108
02109
02110 G4double G4PAIxSection::GetResonanceEnergyTransfer()
02111 {
02112 G4int iTransfer ;
02113
02114 G4double energyTransfer, position;
02115
02116 position = fIntegralResonance[1]*G4UniformRand();
02117
02118 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
02119 {
02120 if( position >= fIntegralResonance[iTransfer] ) break;
02121 }
02122 if(iTransfer > fSplineNumber) iTransfer--;
02123
02124 energyTransfer = fSplineEnergy[iTransfer];
02125
02126 if(iTransfer > 1)
02127 {
02128 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
02129 }
02130 return energyTransfer;
02131 }
02132
02133
02135
02136
02137
02138 G4double G4PAIxSection::GetRutherfordEnergyTransfer()
02139 {
02140 G4int iTransfer ;
02141
02142 G4double energyTransfer, position;
02143
02144 position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
02145
02146 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
02147 {
02148 if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
02149 }
02150 if(iTransfer > fSplineNumber) iTransfer--;
02151
02152 energyTransfer = fSplineEnergy[iTransfer];
02153
02154 if(iTransfer > 1)
02155 {
02156 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
02157 }
02158 return energyTransfer;
02159 }
02160
02162
02163
02164 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const
02165 {
02166 G4String head = "G4PAIxSection::" + methodName + "()";
02167 G4ExceptionDescription ed;
02168 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber << G4endl;
02169 G4Exception(head,"pai001",FatalException,ed);
02170 }
02171
02173
02174
02175
02176
02177 G4int G4PAIxSection::fNumberOfGammas = 111;
02178
02179 const G4double G4PAIxSection::fLorentzFactor[112] =
02180 {
02181 0.0,
02182 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
02183 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00,
02184 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
02185 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00,
02186 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
02187 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00,
02188 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
02189 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01,
02190 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
02191 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01,
02192 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
02193 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02,
02194 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
02195 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02,
02196 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
02197 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03,
02198 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
02199 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03,
02200 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
02201 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04,
02202 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
02203 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04,
02204 1.065799e+05
02205 };
02206
02208
02209
02210
02211
02212 const
02213 G4int G4PAIxSection::fRefGammaNumber = 29;
02214
02215
02216
02217
02218
02220