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 #include "G4PAIySection.hh"
00048
00049 #include "globals.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4ios.hh"
00053 #include "G4Poisson.hh"
00054 #include "G4Material.hh"
00055 #include "G4MaterialCutsCouple.hh"
00056 #include "G4SandiaTable.hh"
00057
00058 using namespace std;
00059
00060
00061
00062 const G4double G4PAIySection::fDelta = 0.005;
00063 const G4double G4PAIySection::fError = 0.005;
00064
00065 const G4int G4PAIySection::fMaxSplineSize = 500;
00066
00067
00069
00070
00071
00072
00073 G4PAIySection::G4PAIySection()
00074 {
00075 fSandia = 0;
00076 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
00077 fIntervalNumber = fSplineNumber = 0;
00078 fVerbose = 0;
00079 for(G4int i=0; i<500; ++i) {
00080 fSplineEnergy[i] = 0.0;
00081 fRePartDielectricConst[i] = 0.0;
00082 fImPartDielectricConst[i] = 0.0;
00083 fIntegralTerm[i] = 0.0;
00084 fDifPAIySection[i] = 0.0;
00085 fdNdxCerenkov[i] = 0.0;
00086 fdNdxPlasmon[i] = 0.0;
00087 fIntegralPAIySection[i] = 0.0;
00088 fIntegralPAIdEdx[i] = 0.0;
00089 fIntegralCerenkov[i] = 0.0;
00090 fIntegralPlasmon[i] = 0.0;
00091 for(G4int j=0; j<112; ++j) { fPAItable[i][j] = 0.0; }
00092 }
00093 }
00094
00096
00097
00098
00099 G4PAIySection::~G4PAIySection()
00100 {}
00101
00103
00104
00105
00106 void G4PAIySection::Initialize( const G4Material* material,
00107 G4double maxEnergyTransfer,
00108 G4double betaGammaSq)
00109 {
00110 G4int i, j;
00111 G4double energy;
00112
00113 fDensity = material->GetDensity();
00114 fElectronDensity = material->GetElectronDensity();
00115
00116
00117 fSandia = material->GetSandiaTable();
00118
00119 fIntervalNumber = fSandia->GetMaxInterval();
00120
00121 fIntervalNumber--;
00122
00123 for( i = 1; i <= fIntervalNumber; i++ )
00124 {
00125 energy = fSandia->GetSandiaMatTablePAI(i-1,0);
00126
00127 if( energy >= maxEnergyTransfer || i > fIntervalNumber )
00128 {
00129 fEnergyInterval[i] = maxEnergyTransfer;
00130 fIntervalNumber = i;
00131 break;
00132 }
00133 fEnergyInterval[i] = energy;
00134 fA1[i] = fSandia->GetSandiaMatTablePAI(i-1,1);
00135 fA2[i] = fSandia->GetSandiaMatTablePAI(i-1,2);
00136 fA3[i] = fSandia->GetSandiaMatTablePAI(i-1,3);
00137 fA4[i] = fSandia->GetSandiaMatTablePAI(i-1,4);
00138
00139 if( fVerbose > 0 && std::fabs( betaGammaSq - 8. ) < 0.4 )
00140 {
00141 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<" keV \t"<<fA1[i]<<"\t"<<fA2[i] <<"\t"<<fA3[i] <<"\t"<<fA4[i]<<G4endl;
00142 }
00143 }
00144
00145
00146 if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
00147 {
00148 fIntervalNumber++;
00149 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
00150 fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
00151 fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
00152 fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
00153 fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
00154 }
00155
00156
00157 for( i = 1; i < fIntervalNumber; i++ )
00158 {
00159
00160
00161 if(fEnergyInterval[i+1]-fEnergyInterval[i] <
00162 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
00163 {
00164 for( j = i; j < fIntervalNumber; j++ )
00165 {
00166 fEnergyInterval[j] = fEnergyInterval[j+1];
00167 fA1[j] = fA1[j+1];
00168 fA2[j] = fA2[j+1];
00169 fA3[j] = fA3[j+1];
00170 fA4[j] = fA4[j+1];
00171 }
00172 fIntervalNumber--;
00173 i--;
00174 }
00175 }
00176 if( fVerbose > 0 && std::fabs( betaGammaSq - 8. ) < 0.4 )
00177 {
00178 G4cout<<"Sandia cofs in G4PAIySection::Initialize(), mat = "<<material->GetName()<<G4endl;
00179 G4cout<<"for bg2 = "<<betaGammaSq<<", Tmax = "<< maxEnergyTransfer/keV<<" keV"<<G4endl;
00180 G4cout<<"energy \t"<<"a1 \t"<<"a2 \t"<<"a3 \t"<<"a4 \t"<<G4endl;
00181
00182 for( j = 1; j < fIntervalNumber; j++ )
00183 {
00184 G4cout<<j<<"\t"<<fEnergyInterval[j]/keV<<" keV \t"<<fA1[j]<<"\t"<<fA2[j] <<"\t"<<fA3[j] <<"\t"<<fA4[j]<<G4endl;
00185 }
00186 }
00187
00188
00189
00190 G4double betaGammaSqRef =
00191 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
00192
00193 ComputeLowEnergyCof(material);
00194 NormShift(betaGammaSqRef);
00195 SplainPAI(betaGammaSqRef);
00196
00197
00198
00199 for( i = 1; i <= fSplineNumber; i++ )
00200 {
00201 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
00202 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00203 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00204 }
00205 IntegralPAIySection();
00206 IntegralCerenkov();
00207 IntegralPlasmon();
00208 }
00209
00211
00212
00213
00214
00215 void G4PAIySection::ComputeLowEnergyCof(const G4Material* material)
00216 {
00217 G4int i, numberOfElements = material->GetNumberOfElements();
00218 G4double sumZ = 0., sumCof = 0.;
00219
00220 const G4double p0 = 1.20923e+00;
00221 const G4double p1 = 3.53256e-01;
00222 const G4double p2 = -1.45052e-03;
00223
00224 G4double* thisMaterialZ = new G4double[numberOfElements];
00225 G4double* thisMaterialCof = new G4double[numberOfElements];
00226
00227 for( i = 0; i < numberOfElements; i++ )
00228 {
00229 thisMaterialZ[i] = material->GetElement(i)->GetZ();
00230 sumZ += thisMaterialZ[i];
00231 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
00232 }
00233 for( i = 0; i < numberOfElements; i++ )
00234 {
00235 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
00236 }
00237 fLowEnergyCof = sumCof;
00238 delete [] thisMaterialZ;
00239 delete [] thisMaterialCof;
00240
00241 }
00242
00243
00244
00245
00247
00248
00249
00250
00251 void G4PAIySection::InitPAI()
00252 {
00253 G4int i;
00254 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
00255 fLorentzFactor[fRefGammaNumber] - 1;
00256
00257
00258
00259 NormShift(betaGammaSq);
00260 SplainPAI(betaGammaSq);
00261
00262 IntegralPAIySection();
00263 IntegralCerenkov();
00264 IntegralPlasmon();
00265
00266 for(i = 0; i<=fSplineNumber; i++)
00267 {
00268 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i];
00269 if(i != 0)
00270 {
00271 fPAItable[i][0] = fSplineEnergy[i];
00272 }
00273 }
00274 fPAItable[0][0] = fSplineNumber;
00275
00276 for(G4int j = 1; j < 112; j++)
00277 {
00278 if( j == fRefGammaNumber ) continue;
00279
00280 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
00281
00282 for(i = 1; i <= fSplineNumber; i++)
00283 {
00284 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
00285 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00286 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00287 }
00288 IntegralPAIySection();
00289 IntegralCerenkov();
00290 IntegralPlasmon();
00291
00292 for(i = 0; i <= fSplineNumber; i++)
00293 {
00294 fPAItable[i][j] = fIntegralPAIySection[i];
00295 }
00296 }
00297
00298 }
00299
00301
00302
00303
00304
00305 void G4PAIySection::NormShift(G4double betaGammaSq)
00306 {
00307 G4int i, j;
00308
00309 for( i = 1; i <= fIntervalNumber-1; i++ )
00310 {
00311 for( j = 1; j <= 2; j++ )
00312 {
00313 fSplineNumber = (i-1)*2 + j;
00314
00315 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
00316 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
00317
00318
00319 }
00320 }
00321 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
00322
00323 j = 1;
00324
00325 for(i=2;i<=fSplineNumber;i++)
00326 {
00327 if(fSplineEnergy[i]<fEnergyInterval[j+1])
00328 {
00329 fIntegralTerm[i] = fIntegralTerm[i-1] +
00330 RutherfordIntegral(j,fSplineEnergy[i-1],
00331 fSplineEnergy[i] );
00332 }
00333 else
00334 {
00335 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
00336 fEnergyInterval[j+1] );
00337 j++;
00338 fIntegralTerm[i] = fIntegralTerm[i-1] + x +
00339 RutherfordIntegral(j,fEnergyInterval[j],
00340 fSplineEnergy[i] );
00341 }
00342
00343 }
00344 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
00345 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
00346
00347
00348
00349
00350
00351
00352 for(G4int k=1;k<=fIntervalNumber-1;k++)
00353 {
00354 for(j=1;j<=2;j++)
00355 {
00356 i = (k-1)*2 + j;
00357 fImPartDielectricConst[i] = fNormalizationCof*
00358 ImPartDielectricConst(k,fSplineEnergy[i]);
00359 fRePartDielectricConst[i] = fNormalizationCof*
00360 RePartDielectricConst(fSplineEnergy[i]);
00361 fIntegralTerm[i] *= fNormalizationCof;
00362
00363 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq);
00364 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
00365 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
00366 }
00367 }
00368
00369 }
00370
00372
00373
00374
00375
00376
00377 void G4PAIySection::SplainPAI(G4double betaGammaSq)
00378 {
00379 G4int k = 1;
00380 G4int i = 1;
00381
00382 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
00383 {
00384 if(fSplineEnergy[i+1] > fEnergyInterval[k+1])
00385 {
00386 k++;
00387 i++;
00388 continue;
00389 }
00390
00391
00392 fSplineNumber++;
00393
00394 for(G4int j = fSplineNumber; j >= i+2; j-- )
00395 {
00396 fSplineEnergy[j] = fSplineEnergy[j-1];
00397 fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
00398 fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
00399 fIntegralTerm[j] = fIntegralTerm[j-1];
00400
00401 fDifPAIySection[j] = fDifPAIySection[j-1];
00402 fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
00403 fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
00404 }
00405 G4double x1 = fSplineEnergy[i];
00406 G4double x2 = fSplineEnergy[i+1];
00407 G4double yy1 = fDifPAIySection[i];
00408 G4double y2 = fDifPAIySection[i+1];
00409
00410 G4double en1 = sqrt(x1*x2);
00411 fSplineEnergy[i+1] = en1;
00412
00413
00414
00415
00416 G4double a = log10(y2/yy1)/log10(x2/x1);
00417 G4double b = log10(yy1) - a*log10(x1);
00418 G4double y = a*log10(en1) + b;
00419 y = pow(10.,y);
00420
00421
00422
00423 fImPartDielectricConst[i+1] = fNormalizationCof*
00424 ImPartDielectricConst(k,fSplineEnergy[i+1]);
00425 fRePartDielectricConst[i+1] = fNormalizationCof*
00426 RePartDielectricConst(fSplineEnergy[i+1]);
00427 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
00428 RutherfordIntegral(k,fSplineEnergy[i],
00429 fSplineEnergy[i+1]);
00430
00431 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq);
00432 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
00433 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
00434
00435
00436
00437
00438 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y);
00439
00440 if( x < 0 )
00441 {
00442 x = -x;
00443 }
00444 if( x > fError && fSplineNumber < fMaxSplineSize-1 )
00445 {
00446 continue;
00447 }
00448 i += 2;
00449
00450 }
00451
00452 }
00453
00454
00456
00457
00458
00459
00460 G4double G4PAIySection::RutherfordIntegral( G4int k,
00461 G4double x1,
00462 G4double x2 )
00463 {
00464 G4double c1, c2, c3;
00465
00466 c1 = (x2 - x1)/x1/x2;
00467 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
00468 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
00469
00470
00471 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
00472
00473 }
00474
00475
00477
00478
00479
00480
00481 G4double G4PAIySection::ImPartDielectricConst( G4int k ,
00482 G4double energy1 )
00483 {
00484 G4double energy2,energy3,energy4,result;
00485
00486 energy2 = energy1*energy1;
00487 energy3 = energy2*energy1;
00488 energy4 = energy3*energy1;
00489
00490 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
00491 result *=hbarc/energy1;
00492
00493 return result;
00494
00495 }
00496
00497
00499
00500
00501
00502
00503
00504 G4double G4PAIySection::RePartDielectricConst(G4double enb)
00505 {
00506 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
00507 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
00508
00509 x0 = enb;
00510 result = 0;
00511
00512 for(G4int i=1;i<=fIntervalNumber-1;i++)
00513 {
00514 x1 = fEnergyInterval[i];
00515 x2 = fEnergyInterval[i+1];
00516 xx1 = x1 - x0;
00517 xx2 = x2 - x0;
00518 xx12 = xx2/xx1;
00519
00520 if(xx12<0)
00521 {
00522 xx12 = -xx12;
00523 }
00524 xln1 = log(x2/x1);
00525 xln2 = log(xx12);
00526 xln3 = log((x2 + x0)/(x1 + x0));
00527 x02 = x0*x0;
00528 x03 = x02*x0;
00529 x04 = x03*x0;
00530 x05 = x04*x0;
00531 c1 = (x2 - x1)/x1/x2;
00532 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
00533 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
00534
00535 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
00536 result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
00537 result -= fA3[i]*c2/2/x02;
00538 result -= fA4[i]*c3/3/x02;
00539
00540 cof1 = fA1[i]/x02 + fA3[i]/x04;
00541 cof2 = fA2[i]/x03 + fA4[i]/x05;
00542
00543 result += 0.5*(cof1 +cof2)*xln2;
00544 result += 0.5*(cof1 - cof2)*xln3;
00545 }
00546 result *= 2*hbarc/pi;
00547
00548 return result;
00549
00550 }
00551
00553
00554
00555
00556
00557
00558 G4double G4PAIySection::DifPAIySection( G4int i ,
00559 G4double betaGammaSq )
00560 {
00561 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
00562
00563
00564 G4double betaBohr = fine_structure_const;
00565
00566
00567 be2 = betaGammaSq/(1 + betaGammaSq);
00568
00569 beta = sqrt(be2);
00570 cof = 1;
00571 x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
00572
00573 if( betaGammaSq < 0.01 ) x2 = log(be2);
00574 else
00575 {
00576 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
00577 (1/betaGammaSq - fRePartDielectricConst[i]) +
00578 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
00579 }
00580 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
00581 {
00582 x6=0;
00583 }
00584 else
00585 {
00586 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
00587 x5 = -1 - fRePartDielectricConst[i] +
00588 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
00589 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
00590
00591 x7 = atan2(fImPartDielectricConst[i],x3);
00592 x6 = x5 * x7;
00593 }
00594
00595
00596 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
00597
00598 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
00599 fImPartDielectricConst[i]*fImPartDielectricConst[i];
00600
00601 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
00602 if(result < 1.0e-8) result = 1.0e-8;
00603 result *= fine_structure_const/be2/pi;
00604
00605
00606 G4double lowCof = fLowEnergyCof;
00607
00608 result *= (1 - exp(-beta/betaBohr/lowCof));
00609
00610
00611
00612
00613
00614 if(x8 > 0.)
00615 {
00616 result /= x8;
00617 }
00618 return result;
00619
00620 }
00621
00623
00624
00625
00626 G4double G4PAIySection::PAIdNdxCerenkov( G4int i ,
00627 G4double betaGammaSq )
00628 {
00629 G4double logarithm, x3, x5, argument, modul2, dNdxC;
00630 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
00631
00632
00633 cofBetaBohr = 4.0;
00634 betaBohr2 = fine_structure_const*fine_structure_const;
00635 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
00636
00637 be2 = betaGammaSq/(1 + betaGammaSq);
00638 be4 = be2*be2;
00639
00640 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq);
00641 else
00642 {
00643 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
00644 (1/betaGammaSq - fRePartDielectricConst[i]) +
00645 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
00646 logarithm += log(1+1.0/betaGammaSq);
00647 }
00648
00649 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
00650 {
00651 argument = 0.0;
00652 }
00653 else
00654 {
00655 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
00656 x5 = -1.0 - fRePartDielectricConst[i] +
00657 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
00658 fImPartDielectricConst[i]*fImPartDielectricConst[i]);
00659 if( x3 == 0.0 ) argument = 0.5*pi;
00660 else argument = atan2(fImPartDielectricConst[i],x3);
00661 argument *= x5 ;
00662 }
00663 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
00664
00665 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
00666
00667 dNdxC *= fine_structure_const/be2/pi;
00668
00669 dNdxC *= (1-exp(-be4/betaBohr4));
00670
00671
00672
00673 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
00674 fImPartDielectricConst[i]*fImPartDielectricConst[i];
00675 if(modul2 > 0.)
00676 {
00677 dNdxC /= modul2;
00678 }
00679 return dNdxC;
00680
00681 }
00682
00684
00685
00686
00687
00688 G4double G4PAIySection::PAIdNdxPlasmon( G4int i ,
00689 G4double betaGammaSq )
00690 {
00691 G4double cof, resonance, modul2, dNdxP;
00692 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
00693
00694 cof = 1;
00695 cofBetaBohr = 4.0;
00696 betaBohr2 = fine_structure_const*fine_structure_const;
00697 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
00698
00699 be2 = betaGammaSq/(1 + betaGammaSq);
00700 be4 = be2*be2;
00701
00702 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
00703 resonance *= fImPartDielectricConst[i]/hbarc;
00704
00705
00706 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
00707
00708 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
00709
00710 dNdxP *= fine_structure_const/be2/pi;
00711 dNdxP *= (1-exp(-be4/betaBohr4));
00712
00713
00714
00715 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
00716 fImPartDielectricConst[i]*fImPartDielectricConst[i];
00717 if(modul2 > 0.)
00718 {
00719 dNdxP /= modul2;
00720 }
00721 return dNdxP;
00722
00723 }
00724
00726
00727
00728
00729
00730
00731 void G4PAIySection::IntegralPAIySection()
00732 {
00733 fIntegralPAIySection[fSplineNumber] = 0;
00734 fIntegralPAIdEdx[fSplineNumber] = 0;
00735 fIntegralPAIySection[0] = 0;
00736 G4int k = fIntervalNumber -1;
00737
00738 for(G4int i = fSplineNumber-1; i >= 1; i--)
00739 {
00740 if(fSplineEnergy[i] >= fEnergyInterval[k])
00741 {
00742 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i);
00743 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
00744 }
00745 else
00746 {
00747 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] +
00748 SumOverBorder(i+1,fEnergyInterval[k]);
00749 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
00750 SumOverBorderdEdx(i+1,fEnergyInterval[k]);
00751 k--;
00752 }
00753 }
00754 }
00755
00757
00758
00759
00760
00761
00762 void G4PAIySection::IntegralCerenkov()
00763 {
00764 G4int i, k;
00765 fIntegralCerenkov[fSplineNumber] = 0;
00766 fIntegralCerenkov[0] = 0;
00767 k = fIntervalNumber -1;
00768
00769 for( i = fSplineNumber-1; i >= 1; i-- )
00770 {
00771 if(fSplineEnergy[i] >= fEnergyInterval[k])
00772 {
00773 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
00774
00775 }
00776 else
00777 {
00778 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
00779 SumOverBordCerenkov(i+1,fEnergyInterval[k]);
00780 k--;
00781
00782 }
00783 }
00784
00785 }
00786
00788
00789
00790
00791
00792
00793 void G4PAIySection::IntegralPlasmon()
00794 {
00795 fIntegralPlasmon[fSplineNumber] = 0;
00796 fIntegralPlasmon[0] = 0;
00797 G4int k = fIntervalNumber -1;
00798 for(G4int i=fSplineNumber-1;i>=1;i--)
00799 {
00800 if(fSplineEnergy[i] >= fEnergyInterval[k])
00801 {
00802 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
00803 }
00804 else
00805 {
00806 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
00807 SumOverBordPlasmon(i+1,fEnergyInterval[k]);
00808 k--;
00809 }
00810 }
00811
00812 }
00813
00815
00816
00817
00818
00819
00820 G4double G4PAIySection::SumOverInterval( G4int i )
00821 {
00822 G4double x0,x1,y0,yy1,a,b,c,result;
00823
00824 x0 = fSplineEnergy[i];
00825 x1 = fSplineEnergy[i+1];
00826 y0 = fDifPAIySection[i];
00827 yy1 = fDifPAIySection[i+1];
00828 c = x1/x0;
00829 a = log10(yy1/y0)/log10(c);
00830
00831 b = y0/pow(x0,a);
00832 a += 1;
00833 if(a == 0)
00834 {
00835 result = b*log(x1/x0);
00836 }
00837 else
00838 {
00839 result = y0*(x1*pow(c,a-1) - x0)/a;
00840 }
00841 a++;
00842 if(a == 0)
00843 {
00844 fIntegralPAIySection[0] += b*log(x1/x0);
00845 }
00846 else
00847 {
00848 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
00849 }
00850 return result;
00851
00852 }
00853
00855
00856 G4double G4PAIySection::SumOverIntervaldEdx( G4int i )
00857 {
00858 G4double x0,x1,y0,yy1,a,b,c,result;
00859
00860 x0 = fSplineEnergy[i];
00861 x1 = fSplineEnergy[i+1];
00862 y0 = fDifPAIySection[i];
00863 yy1 = fDifPAIySection[i+1];
00864 c = x1/x0;
00865 a = log10(yy1/y0)/log10(c);
00866
00867 b = y0/pow(x0,a);
00868 a += 2;
00869 if(a == 0)
00870 {
00871 result = b*log(x1/x0);
00872 }
00873 else
00874 {
00875 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
00876 }
00877 return result;
00878
00879 }
00880
00882
00883
00884
00885
00886
00887 G4double G4PAIySection::SumOverInterCerenkov( G4int i )
00888 {
00889 G4double x0,x1,y0,yy1,a,c,result;
00890
00891 x0 = fSplineEnergy[i];
00892 x1 = fSplineEnergy[i+1];
00893 y0 = fdNdxCerenkov[i];
00894 yy1 = fdNdxCerenkov[i+1];
00895
00896
00897
00898 c = x1/x0;
00899 a = log10(yy1/y0)/log10(c);
00900 G4double b = 0.0;
00901 if(a < 20.) b = y0/pow(x0,a);
00902
00903 a += 1.0;
00904 if(a == 0) result = b*log(c);
00905 else result = y0*(x1*pow(c,a-1) - x0)/a;
00906 a += 1.0;
00907
00908 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
00909 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
00910
00911 return result;
00912
00913 }
00914
00916
00917
00918
00919
00920
00921 G4double G4PAIySection::SumOverInterPlasmon( G4int i )
00922 {
00923 G4double x0,x1,y0,yy1,a,c,result;
00924
00925 x0 = fSplineEnergy[i];
00926 x1 = fSplineEnergy[i+1];
00927 y0 = fdNdxPlasmon[i];
00928 yy1 = fdNdxPlasmon[i+1];
00929 c = x1/x0;
00930 a = log10(yy1/y0)/log10(c);
00931
00932 G4double b = 0.0;
00933 if(a < 20.) b = y0/pow(x0,a);
00934
00935 a += 1.0;
00936 if(a == 0) result = b*log(x1/x0);
00937 else result = y0*(x1*pow(c,a-1) - x0)/a;
00938 a += 1.0;
00939
00940 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
00941 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
00942
00943 return result;
00944
00945 }
00946
00948
00949
00950
00951
00952 G4double G4PAIySection::SumOverBorder( G4int i ,
00953 G4double en0 )
00954 {
00955 G4double x0,x1,y0,yy1,a,d,e0,result;
00956
00957 e0 = en0;
00958 x0 = fSplineEnergy[i];
00959 x1 = fSplineEnergy[i+1];
00960 y0 = fDifPAIySection[i];
00961 yy1 = fDifPAIySection[i+1];
00962
00963
00964 d = e0/x0;
00965 a = log10(yy1/y0)/log10(x1/x0);
00966
00967 G4double b = 0.0;
00968 if(a < 20.) b = y0/pow(x0,a);
00969
00970 a += 1;
00971 if(a == 0)
00972 {
00973 result = b*log(x0/e0);
00974 }
00975 else
00976 {
00977 result = y0*(x0 - e0*pow(d,a-1))/a;
00978 }
00979 a++;
00980 if(a == 0)
00981 {
00982 fIntegralPAIySection[0] += b*log(x0/e0);
00983 }
00984 else
00985 {
00986 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
00987 }
00988 x0 = fSplineEnergy[i - 1];
00989 x1 = fSplineEnergy[i - 2];
00990 y0 = fDifPAIySection[i - 1];
00991 yy1 = fDifPAIySection[i - 2];
00992
00993
00994 d = e0/x0;
00995 a = log10(yy1/y0)/log10(x1/x0);
00996
00997 b = y0/pow(x0,a);
00998 a += 1;
00999 if(a == 0)
01000 {
01001 result += b*log(e0/x0);
01002 }
01003 else
01004 {
01005 result += y0*(e0*pow(d,a-1) - x0)/a;
01006 }
01007 a++;
01008 if(a == 0)
01009 {
01010 fIntegralPAIySection[0] += b*log(e0/x0);
01011 }
01012 else
01013 {
01014 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01015 }
01016 return result;
01017
01018 }
01019
01021
01022 G4double G4PAIySection::SumOverBorderdEdx( G4int i ,
01023 G4double en0 )
01024 {
01025 G4double x0,x1,y0,yy1,a,d,e0,result;
01026
01027 e0 = en0;
01028 x0 = fSplineEnergy[i];
01029 x1 = fSplineEnergy[i+1];
01030 y0 = fDifPAIySection[i];
01031 yy1 = fDifPAIySection[i+1];
01032
01033
01034 d = e0/x0;
01035 a = log10(yy1/y0)/log10(x1/x0);
01036
01037 G4double b = 0.0;
01038 if(a < 20.) b = y0/pow(x0,a);
01039
01040 a += 2;
01041 if(a == 0)
01042 {
01043 result = b*log(x0/e0);
01044 }
01045 else
01046 {
01047 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01048 }
01049 x0 = fSplineEnergy[i - 1];
01050 x1 = fSplineEnergy[i - 2];
01051 y0 = fDifPAIySection[i - 1];
01052 yy1 = fDifPAIySection[i - 2];
01053
01054
01055 d = e0/x0;
01056 a = log10(yy1/y0)/log10(x1/x0);
01057
01058 if(a < 20.) b = y0/pow(x0,a);
01059
01060 a += 2;
01061 if(a == 0)
01062 {
01063 result += b*log(e0/x0);
01064 }
01065 else
01066 {
01067 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01068 }
01069 return result;
01070
01071 }
01072
01074
01075
01076
01077
01078 G4double G4PAIySection::SumOverBordCerenkov( G4int i ,
01079 G4double en0 )
01080 {
01081 G4double x0,x1,y0,yy1,a,e0,c,d,result;
01082
01083 e0 = en0;
01084 x0 = fSplineEnergy[i];
01085 x1 = fSplineEnergy[i+1];
01086 y0 = fdNdxCerenkov[i];
01087 yy1 = fdNdxCerenkov[i+1];
01088
01089
01090
01091
01092 c = x1/x0;
01093 d = e0/x0;
01094 a = log10(yy1/y0)/log10(c);
01095
01096 G4double b = 0.0;
01097 if(a < 20.) b = y0/pow(x0,a);
01098
01099 a += 1.0;
01100 if( a == 0 ) result = b*log(x0/e0);
01101 else result = y0*(x0 - e0*pow(d,a-1))/a;
01102 a += 1.0;
01103
01104 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
01105 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01106
01107
01108
01109 x0 = fSplineEnergy[i - 1];
01110 x1 = fSplineEnergy[i - 2];
01111 y0 = fdNdxCerenkov[i - 1];
01112 yy1 = fdNdxCerenkov[i - 2];
01113
01114
01115
01116
01117 c = x1/x0;
01118 d = e0/x0;
01119 a = log10(yy1/y0)/log10(x1/x0);
01120
01121
01122 if(a < 20.) b = y0/pow(x0,a);
01123
01124 if(a > 20.0) b = 0.0;
01125 else b = y0/pow(x0,a);
01126
01127
01128
01129 a += 1.0;
01130 if( a == 0 ) result += b*log(e0/x0);
01131 else result += y0*(e0*pow(d,a-1) - x0 )/a;
01132 a += 1.0;
01133
01134
01135 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
01136 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01137
01138
01139
01140 return result;
01141
01142 }
01143
01145
01146
01147
01148
01149 G4double G4PAIySection::SumOverBordPlasmon( G4int i ,
01150 G4double en0 )
01151 {
01152 G4double x0,x1,y0,yy1,a,c,d,e0,result;
01153
01154 e0 = en0;
01155 x0 = fSplineEnergy[i];
01156 x1 = fSplineEnergy[i+1];
01157 y0 = fdNdxPlasmon[i];
01158 yy1 = fdNdxPlasmon[i+1];
01159
01160 c = x1/x0;
01161 d = e0/x0;
01162 a = log10(yy1/y0)/log10(c);
01163
01164 G4double b = 0.0;
01165 if(a < 20.) b = y0/pow(x0,a);
01166
01167 a += 1.0;
01168 if( a == 0 ) result = b*log(x0/e0);
01169 else result = y0*(x0 - e0*pow(d,a-1))/a;
01170 a += 1.0;
01171
01172 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
01173 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
01174
01175 x0 = fSplineEnergy[i - 1];
01176 x1 = fSplineEnergy[i - 2];
01177 y0 = fdNdxPlasmon[i - 1];
01178 yy1 = fdNdxPlasmon[i - 2];
01179
01180 c = x1/x0;
01181 d = e0/x0;
01182 a = log10(yy1/y0)/log10(c);
01183
01184 if(a < 20.) b = y0/pow(x0,a);
01185
01186 a += 1.0;
01187 if( a == 0 ) result += b*log(e0/x0);
01188 else result += y0*(e0*pow(d,a-1) - x0)/a;
01189 a += 1.0;
01190
01191 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
01192 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
01193
01194 return result;
01195
01196 }
01197
01199
01200
01201
01202 G4double G4PAIySection::GetStepEnergyLoss( G4double step )
01203 {
01204 G4int iTransfer ;
01205 G4long numOfCollisions;
01206 G4double loss = 0.0;
01207 G4double meanNumber, position;
01208
01209
01210
01211
01212
01213 meanNumber = fIntegralPAIySection[1]*step;
01214 numOfCollisions = G4Poisson(meanNumber);
01215
01216
01217
01218 while(numOfCollisions)
01219 {
01220 position = fIntegralPAIySection[1]*G4UniformRand();
01221
01222 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
01223 {
01224 if( position >= fIntegralPAIySection[iTransfer] ) break;
01225 }
01226 loss += fSplineEnergy[iTransfer] ;
01227 numOfCollisions--;
01228 }
01229
01230
01231 return loss;
01232 }
01233
01235
01236
01237
01238 G4double G4PAIySection::GetStepCerenkovLoss( G4double step )
01239 {
01240 G4int iTransfer ;
01241 G4long numOfCollisions;
01242 G4double loss = 0.0;
01243 G4double meanNumber, position;
01244
01245
01246
01247
01248
01249 meanNumber = fIntegralCerenkov[1]*step;
01250 numOfCollisions = G4Poisson(meanNumber);
01251
01252
01253
01254 while(numOfCollisions)
01255 {
01256 position = fIntegralCerenkov[1]*G4UniformRand();
01257
01258 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
01259 {
01260 if( position >= fIntegralCerenkov[iTransfer] ) break;
01261 }
01262 loss += fSplineEnergy[iTransfer] ;
01263 numOfCollisions--;
01264 }
01265
01266
01267 return loss;
01268 }
01269
01271
01272
01273
01274 G4double G4PAIySection::GetStepPlasmonLoss( G4double step )
01275 {
01276 G4int iTransfer ;
01277 G4long numOfCollisions;
01278 G4double loss = 0.0;
01279 G4double meanNumber, position;
01280
01281
01282
01283
01284
01285 meanNumber = fIntegralPlasmon[1]*step;
01286 numOfCollisions = G4Poisson(meanNumber);
01287
01288
01289
01290 while(numOfCollisions)
01291 {
01292 position = fIntegralPlasmon[1]*G4UniformRand();
01293
01294 for( iTransfer=1; iTransfer<=fSplineNumber; iTransfer++ )
01295 {
01296 if( position >= fIntegralPlasmon[iTransfer] ) break;
01297 }
01298 loss += fSplineEnergy[iTransfer] ;
01299 numOfCollisions--;
01300 }
01301
01302
01303 return loss;
01304 }
01305
01307
01308
01309 void G4PAIySection::CallError(G4int i, const G4String& methodName) const
01310 {
01311 G4String head = "G4PAIySection::" + methodName + "()";
01312 G4ExceptionDescription ed;
01313 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber << G4endl;
01314 G4Exception(head,"pai001",FatalException,ed);
01315 }
01316
01318
01319
01320
01321
01322 G4int G4PAIySection::fNumberOfGammas = 111;
01323
01324 const G4double G4PAIySection::fLorentzFactor[112] =
01325 {
01326 0.0,
01327 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00,
01328 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00,
01329 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00,
01330 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00,
01331 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00,
01332 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00,
01333 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00,
01334 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01,
01335 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01,
01336 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01,
01337 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01,
01338 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02,
01339 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02,
01340 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02,
01341 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03,
01342 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03,
01343 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03,
01344 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03,
01345 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04,
01346 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04,
01347 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04,
01348 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04,
01349 1.065799e+05
01350 };
01351
01353
01354
01355
01356
01357 const
01358 G4int G4PAIySection::fRefGammaNumber = 29;
01359
01360
01361
01362
01363
01365