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 #include "G4ForwardXrayTR.hh"
00039
00040 #include "globals.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4Poisson.hh"
00044 #include "G4Material.hh"
00045 #include "G4PhysicsTable.hh"
00046 #include "G4PhysicsVector.hh"
00047 #include "G4PhysicsLinearVector.hh"
00048 #include "G4PhysicsLogVector.hh"
00049 #include "G4ProductionCutsTable.hh"
00050 #include "G4GeometryTolerance.hh"
00051
00052
00053
00054 G4int G4ForwardXrayTR::fSympsonNumber = 100;
00055
00056 G4double G4ForwardXrayTR::fTheMinEnergyTR = 1.0*keV;
00057 G4double G4ForwardXrayTR::fTheMaxEnergyTR = 100.0*keV;
00058 G4double G4ForwardXrayTR::fTheMaxAngle = 1.0e-3;
00059 G4double G4ForwardXrayTR::fTheMinAngle = 5.0e-6;
00060 G4int G4ForwardXrayTR::fBinTR = 50;
00061
00062 G4double G4ForwardXrayTR::fMinProtonTkin = 100.0*GeV;
00063 G4double G4ForwardXrayTR::fMaxProtonTkin = 100.0*TeV;
00064 G4int G4ForwardXrayTR::fTotBin = 50;
00065
00066 G4double G4ForwardXrayTR::fPlasmaCof = 4.0*pi*fine_structure_const*
00067 hbarc*hbarc*hbarc/electron_mass_c2;
00068
00069 G4double G4ForwardXrayTR::fCofTR = fine_structure_const/pi;
00070
00071
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 G4ForwardXrayTR::
00083 G4ForwardXrayTR( const G4String& matName1,
00084 const G4String& matName2,
00085 const G4String& processName )
00086 : G4TransitionRadiation(processName)
00087 {
00088 fPtrGamma = 0;
00089 fGammaCutInKineticEnergy = 0;
00090 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma =
00091 fSigma1 = fSigma2 = 0.0;
00092 fAngleDistrTable = 0;
00093 fEnergyDistrTable = 0;
00094 fMatIndex1 = fMatIndex2 = 0;
00095
00096
00097
00098 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
00099 fMaxProtonTkin, fTotBin );
00100 G4int iMat;
00101 const G4ProductionCutsTable* theCoupleTable=
00102 G4ProductionCutsTable::GetProductionCutsTable();
00103 G4int numOfCouples = theCoupleTable->GetTableSize();
00104
00105 G4bool build = true;
00106
00107 for(iMat=0;iMat<numOfCouples;iMat++)
00108 {
00109 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
00110 if( matName1 == couple->GetMaterial()->GetName() )
00111 {
00112 fMatIndex1 = couple->GetIndex();
00113 break;
00114 }
00115 }
00116 if(iMat == numOfCouples)
00117 {
00118 G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR01", JustWarning, "Invalid first material name in G4ForwardXrayTR constructor!");
00119 build = false;
00120 }
00121
00122 if(build) {
00123 for(iMat=0;iMat<numOfCouples;iMat++)
00124 {
00125 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
00126 if( matName2 == couple->GetMaterial()->GetName() )
00127 {
00128 fMatIndex2 = couple->GetIndex();
00129 break;
00130 }
00131 }
00132 if(iMat == numOfCouples)
00133 {
00134 G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR02", JustWarning, "Invalid second material name in G4ForwardXrayTR constructor!");
00135 build = false;
00136 }
00137 }
00138
00139 if(build) { BuildXrayTRtables(); }
00140 }
00141
00143
00144
00145
00146 G4ForwardXrayTR::
00147 G4ForwardXrayTR( const G4String& processName )
00148 : G4TransitionRadiation(processName)
00149 {
00150 fPtrGamma = 0;
00151 fGammaCutInKineticEnergy = 0;
00152 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma =
00153 fSigma1 = fSigma2 = 0.0;
00154 fAngleDistrTable = 0;
00155 fEnergyDistrTable = 0;
00156 fMatIndex1 = fMatIndex2 = 0;
00157
00158
00159
00160 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
00161 fMaxProtonTkin, fTotBin );
00162 }
00163
00164
00166
00167
00168
00169
00170 G4ForwardXrayTR::~G4ForwardXrayTR()
00171 {
00172 delete fAngleDistrTable;
00173 delete fEnergyDistrTable;
00174 delete fProtonEnergyVector;
00175 }
00176
00177 G4double G4ForwardXrayTR::GetMeanFreePath(const G4Track&,
00178 G4double,
00179 G4ForceCondition* condition)
00180 {
00181 *condition = Forced;
00182 return DBL_MAX;
00183 }
00184
00186
00187
00188
00189 void G4ForwardXrayTR::BuildXrayTRtables()
00190 {
00191 G4int iMat, jMat, iTkin, iTR, iPlace;
00192 const G4ProductionCutsTable* theCoupleTable=
00193 G4ProductionCutsTable::GetProductionCutsTable();
00194 G4int numOfCouples = theCoupleTable->GetTableSize();
00195
00196 fGammaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
00197
00198 fAngleDistrTable = new G4PhysicsTable(2*fTotBin);
00199 fEnergyDistrTable = new G4PhysicsTable(2*fTotBin);
00200
00201
00202 for(iMat=0;iMat<numOfCouples;iMat++)
00203 {
00204 if( iMat != fMatIndex1 && iMat != fMatIndex2 ) continue;
00205
00206 for(jMat=0;jMat<numOfCouples;jMat++)
00207 {
00208 if( iMat == jMat || ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
00209 {
00210 continue;
00211 }
00212 else
00213 {
00214 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
00215 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
00216 const G4Material* mat1 = iCouple->GetMaterial();
00217 const G4Material* mat2 = jCouple->GetMaterial();
00218
00219 fSigma1 = fPlasmaCof*(mat1->GetElectronDensity());
00220 fSigma2 = fPlasmaCof*(mat2->GetElectronDensity());
00221
00222
00223
00224 fGammaTkinCut = 0.0;
00225
00226 if(fGammaTkinCut > fTheMinEnergyTR)
00227 {
00228 fMinEnergyTR = fGammaTkinCut;
00229 }
00230 else
00231 {
00232 fMinEnergyTR = fTheMinEnergyTR;
00233 }
00234 if(fGammaTkinCut > fTheMaxEnergyTR)
00235 {
00236 fMaxEnergyTR = 2.0*fGammaTkinCut;
00237 }
00238 else
00239 {
00240 fMaxEnergyTR = fTheMaxEnergyTR;
00241 }
00242 for(iTkin=0;iTkin<fTotBin;iTkin++)
00243 {
00244 G4PhysicsLogVector*
00245 energyVector = new G4PhysicsLogVector( fMinEnergyTR,
00246 fMaxEnergyTR,
00247 fBinTR );
00248
00249 fGamma = 1.0 + (fProtonEnergyVector->
00250 GetLowEdgeEnergy(iTkin)/proton_mass_c2);
00251
00252 fMaxThetaTR = 10000.0/(fGamma*fGamma);
00253
00254 if(fMaxThetaTR > fTheMaxAngle)
00255 {
00256 fMaxThetaTR = fTheMaxAngle;
00257 }
00258 else
00259 {
00260 if(fMaxThetaTR < fTheMinAngle)
00261 {
00262 fMaxThetaTR = fTheMinAngle;
00263 }
00264 }
00265
00266 G4PhysicsLinearVector*
00267 angleVector = new G4PhysicsLinearVector( 0.0,
00268 fMaxThetaTR,
00269 fBinTR );
00270 G4double energySum = 0.0;
00271 G4double angleSum = 0.0;
00272
00273 energyVector->PutValue(fBinTR-1,energySum);
00274 angleVector->PutValue(fBinTR-1,angleSum);
00275
00276 for(iTR=fBinTR-2;iTR>=0;iTR--)
00277 {
00278 energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
00279 energyVector->GetLowEdgeEnergy(iTR+1));
00280
00281 angleSum += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
00282 angleVector->GetLowEdgeEnergy(iTR+1));
00283
00284 energyVector->PutValue(iTR,energySum);
00285 angleVector ->PutValue(iTR,angleSum);
00286 }
00287
00288
00289 if(jMat < iMat)
00290 {
00291 iPlace = fTotBin+iTkin;
00292 }
00293 else
00294 {
00295 iPlace = iTkin;
00296 }
00297 fEnergyDistrTable->insertAt(iPlace,energyVector);
00298 fAngleDistrTable->insertAt(iPlace,angleVector);
00299 }
00300 }
00301 }
00302 }
00303
00304 }
00305
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 G4double
00317 G4ForwardXrayTR::SpectralAngleTRdensity( G4double energy,
00318 G4double varAngle ) const
00319 {
00320 G4double formationLength1, formationLength2;
00321 formationLength1 = 1.0/
00322 (1.0/(fGamma*fGamma)
00323 + fSigma1/(energy*energy)
00324 + varAngle);
00325 formationLength2 = 1.0/
00326 (1.0/(fGamma*fGamma)
00327 + fSigma2/(energy*energy)
00328 + varAngle);
00329 return (varAngle/energy)*(formationLength1 - formationLength2)
00330 *(formationLength1 - formationLength2);
00331
00332 }
00333
00334
00336
00337
00338
00339
00340 G4double G4ForwardXrayTR::AngleDensity( G4double energy,
00341 G4double varAngle ) const
00342 {
00343 G4double x, x2, c, d, f, a2, b2, a4, b4;
00344 G4double cof1, cof2, cof3;
00345 x = 1.0/energy;
00346 x2 = x*x;
00347 c = 1.0/fSigma1;
00348 d = 1.0/fSigma2;
00349 f = (varAngle + 1.0/(fGamma*fGamma));
00350 a2 = c*f;
00351 b2 = d*f;
00352 a4 = a2*a2;
00353 b4 = b2*b2;
00354
00355
00356 cof1 = c*c*(0.5/(a2*(x2 +a2)) +0.5*std::log(x2/(x2 +a2))/a4);
00357 cof3 = d*d*(0.5/(b2*(x2 +b2)) +0.5*std::log(x2/(x2 +b2))/b4);
00358 cof2 = -c*d*(std::log(x2/(x2 +b2))/b2 - std::log(x2/(x2 +a2))/a2)/(a2 - b2) ;
00359 return -varAngle*(cof1 + cof2 + cof3);
00360 }
00361
00363
00364
00365
00366
00367
00368 G4double G4ForwardXrayTR::EnergyInterval( G4double energy1,
00369 G4double energy2,
00370 G4double varAngle ) const
00371 {
00372 return AngleDensity(energy2,varAngle)
00373 - AngleDensity(energy1,varAngle);
00374 }
00375
00377
00378
00379
00380
00381
00382 G4double G4ForwardXrayTR::AngleSum( G4double varAngle1,
00383 G4double varAngle2 ) const
00384 {
00385 G4int i;
00386 G4double h , sumEven = 0.0 , sumOdd = 0.0;
00387 h = 0.5*(varAngle2 - varAngle1)/fSympsonNumber;
00388 for(i=1;i<fSympsonNumber;i++)
00389 {
00390 sumEven += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1 + 2*i*h );
00391 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
00392 varAngle1 + (2*i - 1)*h );
00393 }
00394 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
00395 varAngle1 + (2*fSympsonNumber - 1)*h );
00396
00397 return h*(EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1)
00398 + EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle2)
00399 + 4.0*sumOdd + 2.0*sumEven )/3.0;
00400 }
00401
00403
00404
00405
00406
00407
00408 G4double G4ForwardXrayTR::SpectralDensity( G4double energy,
00409 G4double x ) const
00410 {
00411 G4double a, b;
00412 a = 1.0/(fGamma*fGamma)
00413 + fSigma1/(energy*energy) ;
00414 b = 1.0/(fGamma*fGamma)
00415 + fSigma2/(energy*energy) ;
00416 return ( (a + b)*std::log((x + b)/(x + a))/(a - b)
00417 + a/(x + a) + b/(x + b) )/energy;
00418
00419 }
00420
00422
00423
00424
00425
00426
00427 G4double G4ForwardXrayTR::AngleInterval( G4double energy,
00428 G4double varAngle1,
00429 G4double varAngle2 ) const
00430 {
00431 return SpectralDensity(energy,varAngle2)
00432 - SpectralDensity(energy,varAngle1);
00433 }
00434
00436
00437
00438
00439
00440
00441 G4double G4ForwardXrayTR::EnergySum( G4double energy1,
00442 G4double energy2 ) const
00443 {
00444 G4int i;
00445 G4double h , sumEven = 0.0 , sumOdd = 0.0;
00446 h = 0.5*(energy2 - energy1)/fSympsonNumber;
00447 for(i=1;i<fSympsonNumber;i++)
00448 {
00449 sumEven += AngleInterval(energy1 + 2*i*h,0.0,fMaxThetaTR);
00450 sumOdd += AngleInterval(energy1 + (2*i - 1)*h,0.0,fMaxThetaTR);
00451 }
00452 sumOdd += AngleInterval(energy1 + (2*fSympsonNumber - 1)*h,
00453 0.0,fMaxThetaTR);
00454
00455 return h*( AngleInterval(energy1,0.0,fMaxThetaTR)
00456 + AngleInterval(energy2,0.0,fMaxThetaTR)
00457 + 4.0*sumOdd + 2.0*sumEven )/3.0;
00458 }
00459
00461
00462
00463
00464
00465
00466 G4VParticleChange* G4ForwardXrayTR::PostStepDoIt(const G4Track& aTrack,
00467 const G4Step& aStep)
00468 {
00469 aParticleChange.Initialize(aTrack);
00470
00471 G4int iMat, jMat, iTkin, iPlace, numOfTR, iTR, iTransfer;
00472
00473 G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ;
00474 G4double W, W1, W2, E1, E2;
00475
00476 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
00477 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
00478 G4double tol=0.5*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00479
00480 if (pPostStepPoint->GetStepStatus() != fGeomBoundary)
00481 {
00482 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00483 }
00484 if (aTrack.GetStepLength() <= tol)
00485 {
00486 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00487 }
00488
00489
00490 const G4MaterialCutsCouple* iCouple = pPreStepPoint ->GetPhysicalVolume()->
00491 GetLogicalVolume()->GetMaterialCutsCouple();
00492 const G4MaterialCutsCouple* jCouple = pPostStepPoint ->GetPhysicalVolume()->
00493 GetLogicalVolume()->GetMaterialCutsCouple();
00494 const G4Material* iMaterial = iCouple->GetMaterial();
00495 const G4Material* jMaterial = jCouple->GetMaterial();
00496 iMat = iCouple->GetIndex();
00497 jMat = jCouple->GetIndex();
00498
00499
00500
00501
00502 if ( iMat == jMat
00503 || ( (fMatIndex1 >= 0 && fMatIndex1 >= 0)
00504 && ( iMat != fMatIndex1 && iMat != fMatIndex2 )
00505 && ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
00506
00507 || iMaterial->GetState() == jMaterial->GetState()
00508
00509 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
00510
00511 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
00512 {
00513 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00514 }
00515
00516 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00517 G4double charge = aParticle->GetDefinition()->GetPDGCharge();
00518
00519 if(charge == 0.0)
00520 {
00521 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00522 }
00523
00524
00525 G4double chargeSq = charge*charge;
00526 G4double kinEnergy = aParticle->GetKineticEnergy();
00527 G4double massRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass();
00528 G4double TkinScaled = kinEnergy*massRatio;
00529 for(iTkin=0;iTkin<fTotBin;iTkin++)
00530 {
00531 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
00532 {
00533 break;
00534 }
00535 }
00536 if(jMat < iMat)
00537 {
00538 iPlace = fTotBin + iTkin - 1;
00539 }
00540 else
00541 {
00542 iPlace = iTkin - 1;
00543 }
00544
00545
00546
00547
00548
00549
00550 G4ParticleMomentum particleDir = aParticle->GetMomentumDirection();
00551
00552 if(iTkin == fTotBin)
00553 {
00554
00555
00556
00557
00558 numOfTR = G4Poisson( ( (*(*fEnergyDistrTable)(iPlace))(0) +
00559 (*(*fAngleDistrTable)(iPlace))(0) )
00560 *chargeSq*0.5 );
00561 if(numOfTR == 0)
00562 {
00563 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00564 }
00565 else
00566 {
00567
00568
00569 aParticleChange.SetNumberOfSecondaries(numOfTR);
00570
00571 for(iTR=0;iTR<numOfTR;iTR++)
00572 {
00573 energyPos = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand();
00574 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00575 {
00576 if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break;
00577 }
00578 energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
00579
00580
00581
00582 kinEnergy -= energyTR;
00583 aParticleChange.ProposeEnergy(kinEnergy);
00584
00585 anglePos = (*(*fAngleDistrTable)(iPlace))(0)*G4UniformRand();
00586 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00587 {
00588 if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer)) break;
00589 }
00590 theta = std::sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1));
00591
00592
00593
00594 phi = twopi*G4UniformRand();
00595 dirX = std::sin(theta)*std::cos(phi) ;
00596 dirY = std::sin(theta)*std::sin(phi) ;
00597 dirZ = std::cos(theta) ;
00598 G4ThreeVector directionTR(dirX,dirY,dirZ);
00599 directionTR.rotateUz(particleDir);
00600 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
00601 directionTR,
00602 energyTR );
00603 aParticleChange.AddSecondary(aPhotonTR);
00604 }
00605 }
00606 }
00607 else
00608 {
00609 if(iTkin == 0)
00610 {
00611 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00612 }
00613 else
00614 {
00615 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
00616 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
00617 W = 1.0/(E2 - E1);
00618 W1 = (E2 - TkinScaled)*W;
00619 W2 = (TkinScaled - E1)*W;
00620
00621
00622
00623
00624
00625
00626
00627 numOfTR = G4Poisson((((*(*fEnergyDistrTable)(iPlace))(0)+
00628 (*(*fAngleDistrTable)(iPlace))(0))*W1 +
00629 ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
00630 (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
00631 *chargeSq*0.5 );
00632 if(numOfTR == 0)
00633 {
00634 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00635 }
00636 else
00637 {
00638
00639
00640 aParticleChange.SetNumberOfSecondaries(numOfTR);
00641 for(iTR=0;iTR<numOfTR;iTR++)
00642 {
00643 energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+
00644 (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand();
00645 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00646 {
00647 if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+
00648 (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break;
00649 }
00650 energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+
00651 ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2;
00652
00653
00654
00655 kinEnergy -= energyTR;
00656 aParticleChange.ProposeEnergy(kinEnergy);
00657
00658 anglePos = ((*(*fAngleDistrTable)(iPlace))(0)*W1+
00659 (*(*fAngleDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand();
00660 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00661 {
00662 if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer)*W1+
00663 (*(*fAngleDistrTable)(iPlace + 1))(iTransfer)*W2)) break;
00664 }
00665 theta = std::sqrt(((*fAngleDistrTable)(iPlace)->
00666 GetLowEdgeEnergy(iTransfer-1))*W1+
00667 ((*fAngleDistrTable)(iPlace + 1)->
00668 GetLowEdgeEnergy(iTransfer-1))*W2);
00669
00670
00671
00672 phi = twopi*G4UniformRand();
00673 dirX = std::sin(theta)*std::cos(phi) ;
00674 dirY = std::sin(theta)*std::sin(phi) ;
00675 dirZ = std::cos(theta) ;
00676 G4ThreeVector directionTR(dirX,dirY,dirZ);
00677 directionTR.rotateUz(particleDir);
00678 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
00679 directionTR,
00680 energyTR );
00681 aParticleChange.AddSecondary(aPhotonTR);
00682 }
00683 }
00684 }
00685 }
00686 return &aParticleChange;
00687 }
00688
00690
00691
00692
00693
00694
00695 G4double
00696 G4ForwardXrayTR::GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const
00697 {
00698 G4int iPlace, numOfTR, iTR, iTransfer;
00699 G4double energyTR = 0.0;
00700 G4double energyPos ;
00701 G4double W1, W2;
00702
00703 const G4ProductionCutsTable* theCoupleTable=
00704 G4ProductionCutsTable::GetProductionCutsTable();
00705 G4int numOfCouples = theCoupleTable->GetTableSize();
00706
00707
00708
00709
00710 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
00711 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
00712 const G4Material* iMaterial = iCouple->GetMaterial();
00713 const G4Material* jMaterial = jCouple->GetMaterial();
00714
00715 if ( iMat == jMat
00716
00717 || iMaterial->GetState() == jMaterial->GetState()
00718
00719 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
00720
00721 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
00722
00723 {
00724 return energyTR;
00725 }
00726
00727 if(jMat < iMat)
00728 {
00729 iPlace = (iMat*(numOfCouples - 1) + jMat)*fTotBin + iTkin - 1;
00730 }
00731 else
00732 {
00733 iPlace = (iMat*(numOfCouples - 1) + jMat - 1)*fTotBin + iTkin - 1;
00734 }
00735 G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ;
00736 G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1);
00737
00738 if(iTkin == fTotBin)
00739 {
00740 numOfTR = G4Poisson( (*energyVector1)(0) );
00741 if(numOfTR == 0)
00742 {
00743 return energyTR;
00744 }
00745 else
00746 {
00747 for(iTR=0;iTR<numOfTR;iTR++)
00748 {
00749 energyPos = (*energyVector1)(0)*G4UniformRand();
00750 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00751 {
00752 if(energyPos >= (*energyVector1)(iTransfer)) break;
00753 }
00754 energyTR += energyVector1->GetLowEdgeEnergy(iTransfer);
00755 }
00756 }
00757 }
00758 else
00759 {
00760 if(iTkin == 0)
00761 {
00762 return energyTR;
00763 }
00764 else
00765 {
00766 W1 = 0.5;
00767 W2 = 0.5;
00768 numOfTR = G4Poisson( (*energyVector1)(0)*W1 +
00769 (*energyVector2)(0)*W2 );
00770 if(numOfTR == 0)
00771 {
00772 return energyTR;
00773 }
00774 else
00775 {
00776 G4cout<<"It is still OK in GetEnergyTR(int,int,int)"<<G4endl;
00777 for(iTR=0;iTR<numOfTR;iTR++)
00778 {
00779 energyPos = ((*energyVector1)(0)*W1+
00780 (*energyVector2)(0)*W2)*G4UniformRand();
00781 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
00782 {
00783 if(energyPos >= ((*energyVector1)(iTransfer)*W1+
00784 (*energyVector2)(iTransfer)*W2)) break;
00785 }
00786 energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer))*W1+
00787 (energyVector2->GetLowEdgeEnergy(iTransfer))*W2;
00788
00789 }
00790 }
00791 }
00792 }
00793
00794 return energyTR;
00795 }
00796
00798
00799
00800
00801
00802
00803
00804 G4double
00805 G4ForwardXrayTR::GetThetaTR(G4int, G4int, G4int) const
00806 {
00807 G4double theta = 0.0;
00808
00809 return theta;
00810 }
00811
00812
00813
00814
00815