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 #include "G4InitXscPAI.hh"
00044
00045 #include "globals.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4ios.hh"
00049 #include "G4Poisson.hh"
00050 #include "G4Integrator.hh"
00051 #include "G4Material.hh"
00052 #include "G4MaterialCutsCouple.hh"
00053 #include "G4SandiaTable.hh"
00054
00055
00056
00057
00058
00059 const G4double G4InitXscPAI::fDelta = 0.005 ;
00060 const G4int G4InitXscPAI::fPAIbin = 100 ;
00061 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ;
00062
00064
00065
00066
00067
00068 using namespace std;
00069
00070 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC)
00071 : fPAIxscVector(NULL),
00072 fPAIdEdxVector(NULL),
00073 fPAIphotonVector(NULL),
00074 fPAIelectronVector(NULL),
00075 fChCosSqVector(NULL),
00076 fChWidthVector(NULL)
00077 {
00078 G4int i, j, matIndex;
00079
00080 fDensity = matCC->GetMaterial()->GetDensity();
00081 fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
00082 matIndex = matCC->GetMaterial()->GetIndex();
00083
00084 fSandia = new G4SandiaTable(matIndex);
00085 fIntervalNumber = fSandia->GetMaxInterval()-1;
00086
00087 fMatSandiaMatrix = new G4OrderedTable();
00088
00089 for (i = 0; i < fIntervalNumber; i++)
00090 {
00091 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
00092 }
00093 for (i = 0; i < fIntervalNumber; i++)
00094 {
00095 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
00096
00097 for(j = 1; j < 5 ; j++)
00098 {
00099 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
00100 }
00101 }
00102 KillCloseIntervals();
00103 Normalisation();
00104 fBetaGammaSq = fTmax = 0.0;
00105 fIntervalTmax = fCurrentInterval = 0;
00106 }
00107
00108
00109
00110
00112
00113
00114
00115 G4InitXscPAI::~G4InitXscPAI()
00116 {
00117 if(fPAIxscVector) delete fPAIxscVector;
00118 if(fPAIdEdxVector) delete fPAIdEdxVector;
00119 if(fPAIphotonVector) delete fPAIphotonVector;
00120 if(fPAIelectronVector) delete fPAIelectronVector;
00121 if(fChCosSqVector) delete fChCosSqVector;
00122 if(fChWidthVector) delete fChWidthVector;
00123 delete fSandia;
00124 delete fMatSandiaMatrix;
00125 }
00126
00128
00129
00130
00131 void G4InitXscPAI::KillCloseIntervals()
00132 {
00133 G4int i, j, k;
00134 G4double energy1, energy2;
00135
00136 for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
00137 {
00138 energy1 = (*(*fMatSandiaMatrix)[i])[0];
00139 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00140
00141 if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ;
00142 else
00143 {
00144 for(j = i; j < fIntervalNumber-1; j++)
00145 {
00146 for( k = 0; k < 5; k++ )
00147 {
00148 (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
00149 }
00150 }
00151 fIntervalNumber-- ;
00152 i-- ;
00153 }
00154 }
00155
00156 }
00157
00159
00160
00161
00162 void G4InitXscPAI::Normalisation()
00163 {
00164 G4int i, j;
00165 G4double energy1, energy2, cof;
00166
00167 energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
00168 energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
00169
00170
00171 cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
00172
00173 for( i = fIntervalNumber-2; i >= 0; i-- )
00174 {
00175 energy1 = (*(*fMatSandiaMatrix)[i])[0];
00176 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00177
00178 cof += RutherfordIntegral(i,energy1,energy2);
00179
00180 }
00181 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
00182 fNormalizationCof *= fElectronDensity;
00183
00184 fNormalizationCof /= cof;
00185
00186
00187
00188 for (i = 0; i < fIntervalNumber; i++)
00189 {
00190 for(j = 1; j < 5 ; j++)
00191 {
00192 (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
00193 }
00194 }
00195
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
00221
00222
00223
00224
00225 }
00226
00228
00229
00230
00231
00232 G4double G4InitXscPAI::RutherfordIntegral( G4int k,
00233 G4double x1,
00234 G4double x2 )
00235 {
00236 G4double c1, c2, c3, a1, a2, a3, a4 ;
00237
00238 a1 = (*(*fMatSandiaMatrix)[k])[1];
00239 a2 = (*(*fMatSandiaMatrix)[k])[2];
00240 a3 = (*(*fMatSandiaMatrix)[k])[3];
00241 a4 = (*(*fMatSandiaMatrix)[k])[4];
00242
00243 c1 = (x2 - x1)/x1/x2 ;
00244 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
00245 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
00246
00247
00248 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
00249
00250 }
00251
00253
00254
00255
00256 G4double G4InitXscPAI::IntegralTerm(G4double omega)
00257 {
00258 G4int i;
00259 G4double energy1, energy2, result = 0.;
00260
00261 for( i = 0; i <= fIntervalTmax; i++ )
00262 {
00263 if(i == fIntervalTmax)
00264 {
00265 energy1 = (*(*fMatSandiaMatrix)[i])[0];
00266 result += RutherfordIntegral(i,energy1,omega);
00267 }
00268 else
00269 {
00270 if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
00271 {
00272 energy1 = (*(*fMatSandiaMatrix)[i])[0];
00273 result += RutherfordIntegral(i,energy1,omega);
00274 break;
00275 }
00276 else
00277 {
00278 energy1 = (*(*fMatSandiaMatrix)[i])[0];
00279 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00280 result += RutherfordIntegral(i,energy1,energy2);
00281 }
00282 }
00283
00284 }
00285 return result;
00286 }
00287
00288
00290
00291
00292
00293
00294 G4double G4InitXscPAI::ImPartDielectricConst( G4int k ,
00295 G4double energy1 )
00296 {
00297 G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
00298
00299 a1 = (*(*fMatSandiaMatrix)[k])[1];
00300 a2 = (*(*fMatSandiaMatrix)[k])[2];
00301 a3 = (*(*fMatSandiaMatrix)[k])[3];
00302 a4 = (*(*fMatSandiaMatrix)[k])[4];
00303
00304 energy2 = energy1*energy1;
00305 energy3 = energy2*energy1;
00306 energy4 = energy3*energy1;
00307
00308 result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
00309 result *= hbarc/energy1 ;
00310
00311 return result ;
00312
00313 }
00314
00316
00317
00318
00319
00320 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int k ,
00321 G4double omega )
00322 {
00323 G4double eIm2, eRe2, result;
00324
00325 result = ImPartDielectricConst(k,omega);
00326 eIm2 = result*result;
00327
00328 result = RePartDielectricConst(omega);
00329 eRe2 = result*result;
00330
00331 result = eIm2 + eRe2;
00332
00333 return result ;
00334 }
00335
00336
00338
00339
00340
00341
00342
00343 G4double G4InitXscPAI::RePartDielectricConst(G4double enb)
00344 {
00345 G4int i;
00346 G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
00347 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
00348
00349 x0 = enb ;
00350 result = 0 ;
00351
00352 for( i = 0; i < fIntervalNumber-1; i++)
00353 {
00354 x1 = (*(*fMatSandiaMatrix)[i])[0];
00355 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
00356
00357 a1 = (*(*fMatSandiaMatrix)[i])[1];
00358 a2 = (*(*fMatSandiaMatrix)[i])[2];
00359 a3 = (*(*fMatSandiaMatrix)[i])[3];
00360 a4 = (*(*fMatSandiaMatrix)[i])[4];
00361
00362 if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta )
00363 {
00364 if(x0 >= x1) x0 = x1*(1+fDelta);
00365 else x0 = x1*(1-fDelta);
00366 }
00367 if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta )
00368 {
00369 if(x0 >= x2) x0 = x2*(1+fDelta);
00370 else x0 = x2*(1-fDelta);
00371 }
00372 xx1 = x1 - x0 ;
00373 xx2 = x2 - x0 ;
00374 xx12 = xx2/xx1 ;
00375
00376 if( xx12 < 0 ) xx12 = -xx12;
00377
00378 xln1 = log(x2/x1) ;
00379 xln2 = log(xx12) ;
00380 xln3 = log((x2 + x0)/(x1 + x0)) ;
00381
00382 x02 = x0*x0 ;
00383 x03 = x02*x0 ;
00384 x04 = x03*x0 ;
00385 x05 = x04*x0;
00386
00387 c1 = (x2 - x1)/x1/x2 ;
00388 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
00389 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
00390
00391 result -= (a1/x02 + a3/x04)*xln1 ;
00392 result -= (a2/x02 + a4/x04)*c1 ;
00393 result -= a3*c2/2/x02 ;
00394 result -= a4*c3/3/x02 ;
00395
00396 cof1 = a1/x02 + a3/x04 ;
00397 cof2 = a2/x03 + a4/x05 ;
00398
00399 result += 0.5*(cof1 +cof2)*xln2 ;
00400 result += 0.5*(cof1 - cof2)*xln3 ;
00401 }
00402 result *= 2*hbarc/pi ;
00403
00404 return result ;
00405
00406 }
00407
00409
00410
00411
00412
00413
00414 G4double G4InitXscPAI::DifPAIxSection( G4double omega )
00415 {
00416 G4int i = fCurrentInterval;
00417 G4double betaGammaSq = fBetaGammaSq;
00418 G4double integralTerm = IntegralTerm(omega);
00419 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
00420 G4double epsilonRe = RePartDielectricConst(omega);
00421 G4double epsilonIm = ImPartDielectricConst(i,omega);
00422 G4double be4 ;
00423 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
00424 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
00425 be2 = betaGammaSq/(1 + betaGammaSq) ;
00426 be4 = be2*be2 ;
00427
00428 cof = 1 ;
00429 x1 = log(2*electron_mass_c2/omega) ;
00430
00431 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
00432 else
00433 {
00434 x2 = -log( (1/betaGammaSq - epsilonRe)*
00435 (1/betaGammaSq - epsilonRe) +
00436 epsilonIm*epsilonIm )/2 ;
00437 }
00438 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
00439 {
00440 x6=0 ;
00441 }
00442 else
00443 {
00444 x3 = -epsilonRe + 1/betaGammaSq ;
00445 x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
00446 epsilonIm*epsilonIm) ;
00447
00448 x7 = atan2(epsilonIm,x3) ;
00449 x6 = x5 * x7 ;
00450 }
00451
00452
00453 x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
00454
00455 x8 = (1 + epsilonRe)*(1 + epsilonRe) +
00456 epsilonIm*epsilonIm;
00457
00458 result = (x4 + cof*integralTerm/omega/omega) ;
00459 if(result < 1.0e-8) result = 1.0e-8 ;
00460 result *= fine_structure_const/be2/pi ;
00461
00462
00463 result *= (1-exp(-be4/betaBohr4)) ;
00464 if(fDensity >= fSolidDensity)
00465 {
00466 result /= x8 ;
00467 }
00468 return result ;
00469
00470 }
00471
00473
00474
00475
00476
00477 G4double G4InitXscPAI::DifPAIdEdx( G4double omega )
00478 {
00479 G4double dEdx = omega*DifPAIxSection(omega);
00480 return dEdx;
00481 }
00482
00484
00485
00486
00487 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega )
00488 {
00489 G4int i = fCurrentInterval;
00490 G4double betaGammaSq = fBetaGammaSq;
00491 G4double epsilonRe = RePartDielectricConst(omega);
00492 G4double epsilonIm = ImPartDielectricConst(i,omega);
00493
00494 G4double logarithm, x3, x5, argument, modul2, dNdxC ;
00495 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
00496
00497
00498 cofBetaBohr = 4.0 ;
00499 betaBohr2 = fine_structure_const*fine_structure_const ;
00500 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
00501
00502 be2 = betaGammaSq/(1 + betaGammaSq) ;
00503 be4 = be2*be2 ;
00504
00505 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ;
00506 else
00507 {
00508 logarithm = -log( (1/betaGammaSq - epsilonRe)*
00509 (1/betaGammaSq - epsilonRe) +
00510 epsilonIm*epsilonIm )*0.5 ;
00511 logarithm += log(1+1.0/betaGammaSq) ;
00512 }
00513
00514 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
00515 {
00516 argument = 0.0 ;
00517 }
00518 else
00519 {
00520 x3 = -epsilonRe + 1.0/betaGammaSq ;
00521 x5 = -1.0 - epsilonRe +
00522 be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
00523 epsilonIm*epsilonIm) ;
00524 if( x3 == 0.0 ) argument = 0.5*pi;
00525 else argument = atan2(epsilonIm,x3) ;
00526 argument *= x5 ;
00527 }
00528 dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
00529
00530 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
00531
00532 dNdxC *= fine_structure_const/be2/pi ;
00533
00534 dNdxC *= (1-exp(-be4/betaBohr4)) ;
00535
00536 if(fDensity >= fSolidDensity)
00537 {
00538 modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
00539 epsilonIm*epsilonIm;
00540 dNdxC /= modul2 ;
00541 }
00542 return dNdxC ;
00543
00544 }
00545
00547
00548
00549
00550
00551 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega )
00552 {
00553 G4int i = fCurrentInterval;
00554 G4double betaGammaSq = fBetaGammaSq;
00555 G4double integralTerm = IntegralTerm(omega);
00556 G4double epsilonRe = RePartDielectricConst(omega);
00557 G4double epsilonIm = ImPartDielectricConst(i,omega);
00558
00559 G4double cof, resonance, modul2, dNdxP ;
00560 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
00561
00562 cof = 1 ;
00563 cofBetaBohr = 4.0 ;
00564 betaBohr2 = fine_structure_const*fine_structure_const ;
00565 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
00566
00567 be2 = betaGammaSq/(1 + betaGammaSq) ;
00568 be4 = be2*be2 ;
00569
00570 resonance = log(2*electron_mass_c2*be2/omega) ;
00571 resonance *= epsilonIm/hbarc ;
00572
00573
00574 dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
00575
00576 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
00577
00578 dNdxP *= fine_structure_const/be2/pi ;
00579 dNdxP *= (1-exp(-be4/betaBohr4)) ;
00580
00581 if( fDensity >= fSolidDensity )
00582 {
00583 modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
00584 epsilonIm*epsilonIm;
00585 dNdxP /= modul2 ;
00586 }
00587 return dNdxP ;
00588
00589 }
00590
00592
00593
00594
00595
00596
00597 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax)
00598 {
00599 G4int i,k,i1,i2;
00600 G4double energy1, energy2, result = 0.;
00601
00602 fBetaGammaSq = bg2;
00603 fTmax = Tmax;
00604
00605 if(fPAIxscVector) delete fPAIxscVector;
00606
00607 fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00608 fPAIxscVector->PutValue(fPAIbin-1,result);
00609
00610 for( i = fIntervalNumber - 1; i >= 0; i-- )
00611 {
00612 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00613 }
00614 if (i < 0) i = 0;
00615
00616 fIntervalTmax = i;
00617
00618 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00619
00620 for( k = fPAIbin - 2; k >= 0; k-- )
00621 {
00622 energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
00623 energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
00624
00625 for( i = fIntervalTmax; i >= 0; i-- )
00626 {
00627 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00628 }
00629 if(i < 0) i = 0;
00630 i2 = i;
00631
00632 for( i = fIntervalTmax; i >= 0; i-- )
00633 {
00634 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00635 }
00636 if(i < 0) i = 0;
00637 i1 = i;
00638
00639 if( i1 == i2 )
00640 {
00641 fCurrentInterval = i1;
00642 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
00643 energy1,energy2);
00644 fPAIxscVector->PutValue(k,result);
00645 }
00646 else
00647 {
00648 for( i = i2; i >= i1; i-- )
00649 {
00650 fCurrentInterval = i;
00651
00652 if( i==i2 ) result += integral.Legendre10(this,
00653 &G4InitXscPAI::DifPAIxSection,
00654 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00655
00656 else if( i == i1 ) result += integral.Legendre10(this,
00657 &G4InitXscPAI::DifPAIxSection,energy1,
00658 (*(*fMatSandiaMatrix)[i+1])[0]);
00659
00660 else result += integral.Legendre10(this,
00661 &G4InitXscPAI::DifPAIxSection,
00662 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00663 }
00664 fPAIxscVector->PutValue(k,result);
00665 }
00666
00667 }
00668 return ;
00669 }
00670
00671
00673
00674
00675
00676
00677
00678 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax)
00679 {
00680 G4int i,k,i1,i2;
00681 G4double energy1, energy2, result = 0.;
00682
00683 fBetaGammaSq = bg2;
00684 fTmax = Tmax;
00685
00686 if(fPAIdEdxVector) delete fPAIdEdxVector;
00687
00688 fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00689 fPAIdEdxVector->PutValue(fPAIbin-1,result);
00690
00691 for( i = fIntervalNumber - 1; i >= 0; i-- )
00692 {
00693 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00694 }
00695 if (i < 0) i = 0;
00696
00697 fIntervalTmax = i;
00698
00699 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00700
00701 for( k = fPAIbin - 2; k >= 0; k-- )
00702 {
00703 energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
00704 energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
00705
00706 for( i = fIntervalTmax; i >= 0; i-- )
00707 {
00708 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00709 }
00710 if(i < 0) i = 0;
00711 i2 = i;
00712
00713 for( i = fIntervalTmax; i >= 0; i-- )
00714 {
00715 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00716 }
00717 if(i < 0) i = 0;
00718 i1 = i;
00719
00720 if( i1 == i2 )
00721 {
00722 fCurrentInterval = i1;
00723 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
00724 energy1,energy2);
00725 fPAIdEdxVector->PutValue(k,result);
00726 }
00727 else
00728 {
00729 for( i = i2; i >= i1; i-- )
00730 {
00731 fCurrentInterval = i;
00732
00733 if( i==i2 ) result += integral.Legendre10(this,
00734 &G4InitXscPAI::DifPAIdEdx,
00735 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00736
00737 else if( i == i1 ) result += integral.Legendre10(this,
00738 &G4InitXscPAI::DifPAIdEdx,energy1,
00739 (*(*fMatSandiaMatrix)[i+1])[0]);
00740
00741 else result += integral.Legendre10(this,
00742 &G4InitXscPAI::DifPAIdEdx,
00743 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00744 }
00745 fPAIdEdxVector->PutValue(k,result);
00746 }
00747
00748 }
00749 return ;
00750 }
00751
00753
00754
00755
00756
00757
00758 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax)
00759 {
00760 G4int i,k,i1,i2;
00761 G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
00762
00763 fBetaGammaSq = bg2;
00764 fTmax = Tmax;
00765 beta2 = bg2/(1+bg2);
00766
00767 if(fPAIphotonVector) delete fPAIphotonVector;
00768 if(fChCosSqVector) delete fChCosSqVector;
00769 if(fChWidthVector) delete fChWidthVector;
00770
00771 fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00772 fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00773 fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00774
00775 fPAIphotonVector->PutValue(fPAIbin-1,result);
00776 fChCosSqVector->PutValue(fPAIbin-1,1.);
00777 fChWidthVector->PutValue(fPAIbin-1,1e-7);
00778
00779 for( i = fIntervalNumber - 1; i >= 0; i-- )
00780 {
00781 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00782 }
00783 if (i < 0) i = 0;
00784
00785 fIntervalTmax = i;
00786
00787 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00788
00789 for( k = fPAIbin - 2; k >= 0; k-- )
00790 {
00791 energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
00792 energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
00793
00794 for( i = fIntervalTmax; i >= 0; i-- )
00795 {
00796 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00797 }
00798 if(i < 0) i = 0;
00799 i2 = i;
00800
00801 for( i = fIntervalTmax; i >= 0; i-- )
00802 {
00803 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00804 }
00805 if(i < 0) i = 0;
00806 i1 = i;
00807
00808 module2 = ModuleSqDielectricConst(i1,energy1);
00809 cos2 = RePartDielectricConst(energy1)/module2/beta2;
00810 width = ImPartDielectricConst(i1,energy1)/module2/beta2;
00811
00812 fChCosSqVector->PutValue(k,cos2);
00813 fChWidthVector->PutValue(k,width);
00814
00815 if( i1 == i2 )
00816 {
00817 fCurrentInterval = i1;
00818 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
00819 energy1,energy2);
00820 fPAIphotonVector->PutValue(k,result);
00821
00822 }
00823 else
00824 {
00825 for( i = i2; i >= i1; i-- )
00826 {
00827 fCurrentInterval = i;
00828
00829 if( i==i2 ) result += integral.Legendre10(this,
00830 &G4InitXscPAI::PAIdNdxCherenkov,
00831 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00832
00833 else if( i == i1 ) result += integral.Legendre10(this,
00834 &G4InitXscPAI::PAIdNdxCherenkov,energy1,
00835 (*(*fMatSandiaMatrix)[i+1])[0]);
00836
00837 else result += integral.Legendre10(this,
00838 &G4InitXscPAI::PAIdNdxCherenkov,
00839 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00840 }
00841 fPAIphotonVector->PutValue(k,result);
00842 }
00843
00844 }
00845 return;
00846 }
00847
00849
00850
00851
00852
00853
00854 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax)
00855 {
00856 G4int i,k,i1,i2;
00857 G4double energy1, energy2, result = 0.;
00858
00859 fBetaGammaSq = bg2;
00860 fTmax = Tmax;
00861
00862 if(fPAIelectronVector) delete fPAIelectronVector;
00863
00864 fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00865 fPAIelectronVector->PutValue(fPAIbin-1,result);
00866
00867 for( i = fIntervalNumber - 1; i >= 0; i-- )
00868 {
00869 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00870 }
00871 if (i < 0) i = 0;
00872
00873 fIntervalTmax = i;
00874
00875 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00876
00877 for( k = fPAIbin - 2; k >= 0; k-- )
00878 {
00879 energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
00880 energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
00881
00882 for( i = fIntervalTmax; i >= 0; i-- )
00883 {
00884 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00885 }
00886 if(i < 0) i = 0;
00887 i2 = i;
00888
00889 for( i = fIntervalTmax; i >= 0; i-- )
00890 {
00891 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00892 }
00893 if(i < 0) i = 0;
00894 i1 = i;
00895
00896 if( i1 == i2 )
00897 {
00898 fCurrentInterval = i1;
00899 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
00900 energy1,energy2);
00901 fPAIelectronVector->PutValue(k,result);
00902 }
00903 else
00904 {
00905 for( i = i2; i >= i1; i-- )
00906 {
00907 fCurrentInterval = i;
00908
00909 if( i==i2 ) result += integral.Legendre10(this,
00910 &G4InitXscPAI::PAIdNdxPlasmon,
00911 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00912
00913 else if( i == i1 ) result += integral.Legendre10(this,
00914 &G4InitXscPAI::PAIdNdxPlasmon,energy1,
00915 (*(*fMatSandiaMatrix)[i+1])[0]);
00916
00917 else result += integral.Legendre10(this,
00918 &G4InitXscPAI::PAIdNdxPlasmon,
00919 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00920 }
00921 fPAIelectronVector->PutValue(k,result);
00922 }
00923
00924 }
00925 return;
00926 }
00927
00928
00930
00931
00932
00933 G4double G4InitXscPAI::GetPhotonLambda( G4double omega )
00934 {
00935 G4int i ;
00936 G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
00937
00938 omega2 = omega*omega ;
00939 omega3 = omega2*omega ;
00940 omega4 = omega2*omega2 ;
00941
00942 for(i = 0; i < fIntervalNumber;i++)
00943 {
00944 if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
00945 }
00946 if( i == 0 )
00947 {
00948 G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
00949 }
00950 else i-- ;
00951
00952 a1 = (*(*fMatSandiaMatrix)[i])[1];
00953 a2 = (*(*fMatSandiaMatrix)[i])[2];
00954 a3 = (*(*fMatSandiaMatrix)[i])[3];
00955 a4 = (*(*fMatSandiaMatrix)[i])[4];
00956
00957 lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
00958
00959 return lambda ;
00960 }
00961
00963
00964
00965
00967
00968
00969
00970 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step )
00971 {
00972 G4double loss = 0.0 ;
00973 loss *= step;
00974
00975 return loss ;
00976 }
00977
00979
00980
00981
00982 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step )
00983 {
00984 G4double loss = 0.0 ;
00985 loss *= step;
00986
00987 return loss ;
00988 }
00989
00991
00992
00993
00994 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step )
00995 {
00996
00997
00998 G4double loss = 0.0 ;
00999 loss *= step;
01000 return loss ;
01001 }
01002
01003
01004
01005
01006
01008