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
00041
00042 #include "G4SynchrotronRadiationInMat.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4Integrator.hh"
00046 #include "G4EmProcessSubType.hh"
00047
00049
00050
00051
00052
00053 const G4double
00054 G4SynchrotronRadiationInMat::fLambdaConst = std::sqrt(3.0)*electron_mass_c2/
00055 (2.5*fine_structure_const*eplus*c_light) ;
00056
00058
00059
00060
00061
00062 const G4double
00063 G4SynchrotronRadiationInMat::fEnergyConst = 1.5*c_light*c_light*eplus*hbar_Planck/
00064 electron_mass_c2 ;
00065
00067
00068
00069
00070
00071
00072
00073 const G4double
00074 G4SynchrotronRadiationInMat::fIntegralProbabilityOfSR[200] =
00075 {
00076 1.000000e+00, 9.428859e-01, 9.094095e-01, 8.813971e-01, 8.565154e-01,
00077 8.337008e-01, 8.124961e-01, 7.925217e-01, 7.735517e-01, 7.554561e-01,
00078 7.381233e-01, 7.214521e-01, 7.053634e-01, 6.898006e-01, 6.747219e-01,
00079 6.600922e-01, 6.458793e-01, 6.320533e-01, 6.185872e-01, 6.054579e-01,
00080 5.926459e-01, 5.801347e-01, 5.679103e-01, 5.559604e-01, 5.442736e-01,
00081 5.328395e-01, 5.216482e-01, 5.106904e-01, 4.999575e-01, 4.894415e-01,
00082 4.791351e-01, 4.690316e-01, 4.591249e-01, 4.494094e-01, 4.398800e-01,
00083 4.305320e-01, 4.213608e-01, 4.123623e-01, 4.035325e-01, 3.948676e-01,
00084 3.863639e-01, 3.780179e-01, 3.698262e-01, 3.617858e-01, 3.538933e-01,
00085 3.461460e-01, 3.385411e-01, 3.310757e-01, 3.237474e-01, 3.165536e-01,
00086 3.094921e-01, 3.025605e-01, 2.957566e-01, 2.890784e-01, 2.825237e-01,
00087 2.760907e-01, 2.697773e-01, 2.635817e-01, 2.575020e-01, 2.515365e-01,
00088 2.456834e-01, 2.399409e-01, 2.343074e-01, 2.287812e-01, 2.233607e-01,
00089 2.180442e-01, 2.128303e-01, 2.077174e-01, 2.027040e-01, 1.977885e-01,
00090 1.929696e-01, 1.882457e-01, 1.836155e-01, 1.790775e-01, 1.746305e-01,
00091 1.702730e-01, 1.660036e-01, 1.618212e-01, 1.577243e-01, 1.537117e-01,
00092 1.497822e-01, 1.459344e-01, 1.421671e-01, 1.384791e-01, 1.348691e-01,
00093 1.313360e-01, 1.278785e-01, 1.244956e-01, 1.211859e-01, 1.179483e-01,
00094 1.147818e-01, 1.116850e-01, 1.086570e-01, 1.056966e-01, 1.028026e-01,
00095 9.997405e-02, 9.720975e-02, 9.450865e-02, 9.186969e-02, 8.929179e-02,
00096 8.677391e-02, 8.431501e-02, 8.191406e-02, 7.957003e-02, 7.728192e-02,
00097 7.504872e-02, 7.286944e-02, 7.074311e-02, 6.866874e-02, 6.664538e-02,
00098 6.467208e-02, 6.274790e-02, 6.087191e-02, 5.904317e-02, 5.726079e-02,
00099 5.552387e-02, 5.383150e-02, 5.218282e-02, 5.057695e-02, 4.901302e-02,
00100 4.749020e-02, 4.600763e-02, 4.456450e-02, 4.315997e-02, 4.179325e-02,
00101 4.046353e-02, 3.917002e-02, 3.791195e-02, 3.668855e-02, 3.549906e-02,
00102 3.434274e-02, 3.321884e-02, 3.212665e-02, 3.106544e-02, 3.003452e-02,
00103 2.903319e-02, 2.806076e-02, 2.711656e-02, 2.619993e-02, 2.531021e-02,
00104 2.444677e-02, 2.360897e-02, 2.279620e-02, 2.200783e-02, 2.124327e-02,
00105 2.050194e-02, 1.978324e-02, 1.908662e-02, 1.841151e-02, 1.775735e-02,
00106 1.712363e-02, 1.650979e-02, 1.591533e-02, 1.533973e-02, 1.478250e-02,
00107 1.424314e-02, 1.372117e-02, 1.321613e-02, 1.272755e-02, 1.225498e-02,
00108 1.179798e-02, 1.135611e-02, 1.092896e-02, 1.051609e-02, 1.011712e-02,
00109 9.731635e-03, 9.359254e-03, 8.999595e-03, 8.652287e-03, 8.316967e-03,
00110 7.993280e-03, 7.680879e-03, 7.379426e-03, 7.088591e-03, 6.808051e-03,
00111 6.537491e-03, 6.276605e-03, 6.025092e-03, 5.782661e-03, 5.549027e-03,
00112 5.323912e-03, 5.107045e-03, 4.898164e-03, 4.697011e-03, 4.503336e-03,
00113 4.316896e-03, 4.137454e-03, 3.964780e-03, 3.798649e-03, 3.638843e-03,
00114 3.485150e-03, 3.337364e-03, 3.195284e-03, 3.058715e-03, 2.927469e-03,
00115 2.801361e-03, 2.680213e-03, 2.563852e-03, 2.452110e-03, 2.344824e-03
00116 };
00117
00119
00120
00121
00122
00123 G4SynchrotronRadiationInMat::G4SynchrotronRadiationInMat(const G4String& processName,
00124 G4ProcessType type):G4VDiscreteProcess (processName, type),
00125 LowestKineticEnergy (10.*keV),
00126 HighestKineticEnergy (100.*TeV),
00127 TotBin(200),
00128 theGamma (G4Gamma::Gamma() ),
00129 theElectron ( G4Electron::Electron() ),
00130 thePositron ( G4Positron::Positron() ),
00131 GammaCutInKineticEnergy(0),
00132 ElectronCutInKineticEnergy(0),
00133 PositronCutInKineticEnergy(0),
00134 ParticleCutInKineticEnergy(0),
00135 fAlpha(0.0), fRootNumber(80),
00136 fVerboseLevel( verboseLevel )
00137 {
00138 G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
00139
00140 fFieldPropagator = transportMgr->GetPropagatorInField();
00141 SetProcessSubType(fSynchrotronRadiation);
00142 CutInRange = GammaCutInKineticEnergyNow = ElectronCutInKineticEnergyNow =
00143 PositronCutInKineticEnergyNow = ParticleCutInKineticEnergyNow = fKsi =
00144 fPsiGamma = fEta = fOrderAngleK = 0.0;
00145 }
00146
00148
00149
00150
00151
00152 G4SynchrotronRadiationInMat::~G4SynchrotronRadiationInMat()
00153 {}
00154
00155
00156 G4bool
00157 G4SynchrotronRadiationInMat::IsApplicable( const G4ParticleDefinition& particle )
00158 {
00159
00160 return ( ( &particle == (const G4ParticleDefinition *)theElectron ) ||
00161 ( &particle == (const G4ParticleDefinition *)thePositron ) );
00162
00163 }
00164
00166
00167
00168
00169
00170
00171
00172
00173 G4double
00174 G4SynchrotronRadiationInMat::GetMeanFreePath( const G4Track& trackData,
00175 G4double,
00176 G4ForceCondition* condition)
00177 {
00178
00179 G4double MeanFreePath;
00180
00181 const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
00182
00183
00184
00185
00186 *condition = NotForced ;
00187
00188 G4double gamma = aDynamicParticle->GetTotalEnergy()/
00189 aDynamicParticle->GetMass();
00190
00191 G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00192
00193 G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
00194
00195 if ( KineticEnergy < LowestKineticEnergy || gamma < 1.0e3 ) MeanFreePath = DBL_MAX;
00196 else
00197 {
00198
00199 G4ThreeVector FieldValue;
00200 const G4Field* pField = 0;
00201
00202 G4FieldManager* fieldMgr=0;
00203 G4bool fieldExertsForce = false;
00204
00205 if( (particleCharge != 0.0) )
00206 {
00207 fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
00208
00209 if ( fieldMgr != 0 )
00210 {
00211
00212
00213 fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
00214 }
00215 }
00216 if ( fieldExertsForce )
00217 {
00218 pField = fieldMgr->GetDetectorField() ;
00219 G4ThreeVector globPosition = trackData.GetPosition();
00220
00221 G4double globPosVec[4], FieldValueVec[6];
00222
00223 globPosVec[0] = globPosition.x();
00224 globPosVec[1] = globPosition.y();
00225 globPosVec[2] = globPosition.z();
00226 globPosVec[3] = trackData.GetGlobalTime();
00227
00228 pField->GetFieldValue( globPosVec, FieldValueVec );
00229
00230 FieldValue = G4ThreeVector( FieldValueVec[0],
00231 FieldValueVec[1],
00232 FieldValueVec[2] );
00233
00234
00235
00236 G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
00237 G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum) ;
00238 G4double perpB = unitMcrossB.mag() ;
00239 G4double beta = aDynamicParticle->GetTotalMomentum()/
00240 (aDynamicParticle->GetTotalEnergy() );
00241
00242 if( perpB > 0.0 ) MeanFreePath = fLambdaConst*beta/perpB;
00243 else MeanFreePath = DBL_MAX;
00244 }
00245 else MeanFreePath = DBL_MAX;
00246 }
00247 if(fVerboseLevel > 0)
00248 {
00249 G4cout<<"G4SynchrotronRadiationInMat::MeanFreePath = "<<MeanFreePath/m<<" m"<<G4endl;
00250 }
00251 return MeanFreePath;
00252 }
00253
00255
00256
00257
00258 G4VParticleChange*
00259 G4SynchrotronRadiationInMat::PostStepDoIt(const G4Track& trackData,
00260 const G4Step& stepData )
00261
00262 {
00263 aParticleChange.Initialize(trackData);
00264
00265 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
00266
00267 G4double gamma = aDynamicParticle->GetTotalEnergy()/
00268 (aDynamicParticle->GetMass() );
00269
00270 if(gamma <= 1.0e3 )
00271 {
00272 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00273 }
00274 G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00275
00276 G4ThreeVector FieldValue;
00277 const G4Field* pField = 0 ;
00278
00279 G4FieldManager* fieldMgr=0;
00280 G4bool fieldExertsForce = false;
00281
00282 if( (particleCharge != 0.0) )
00283 {
00284 fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
00285 if ( fieldMgr != 0 )
00286 {
00287
00288
00289 fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
00290 }
00291 }
00292 if ( fieldExertsForce )
00293 {
00294 pField = fieldMgr->GetDetectorField() ;
00295 G4ThreeVector globPosition = trackData.GetPosition() ;
00296 G4double globPosVec[4], FieldValueVec[6] ;
00297 globPosVec[0] = globPosition.x() ;
00298 globPosVec[1] = globPosition.y() ;
00299 globPosVec[2] = globPosition.z() ;
00300 globPosVec[3] = trackData.GetGlobalTime();
00301
00302 pField->GetFieldValue( globPosVec, FieldValueVec ) ;
00303 FieldValue = G4ThreeVector( FieldValueVec[0],
00304 FieldValueVec[1],
00305 FieldValueVec[2] );
00306
00307 G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
00308 G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
00309 G4double perpB = unitMcrossB.mag() ;
00310 if(perpB > 0.0)
00311 {
00312
00313
00314 G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
00315
00316 if(fVerboseLevel > 0)
00317 {
00318 G4cout<<"SR photon energy = "<<energyOfSR/keV<<" keV"<<G4endl;
00319 }
00320
00321
00322 if( energyOfSR <= 0.0 )
00323 {
00324 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00325 }
00326 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00327 G4ParticleMomentum
00328 particleDirection = aDynamicParticle->GetMomentumDirection();
00329
00330
00331
00332
00333
00334 G4double cosTheta, sinTheta, fcos, beta;
00335
00336 do
00337 {
00338 cosTheta = 1. - 2.*G4UniformRand();
00339 fcos = (1 + cosTheta*cosTheta)*0.5;
00340 }
00341 while( fcos < G4UniformRand() );
00342
00343 beta = std::sqrt(1. - 1./(gamma*gamma));
00344
00345 cosTheta = (cosTheta + beta)/(1. + beta*cosTheta);
00346
00347 if( cosTheta > 1. ) cosTheta = 1.;
00348 if( cosTheta < -1. ) cosTheta = -1.;
00349
00350 sinTheta = std::sqrt(1. - cosTheta*cosTheta );
00351
00352 G4double Phi = twopi * G4UniformRand() ;
00353
00354 G4double dirx = sinTheta*std::cos(Phi) ,
00355 diry = sinTheta*std::sin(Phi) ,
00356 dirz = cosTheta;
00357
00358 G4ThreeVector gammaDirection ( dirx, diry, dirz);
00359 gammaDirection.rotateUz(particleDirection);
00360
00361
00362
00363
00364
00365
00366
00367 G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
00368 gammaPolarization = gammaPolarization.unit();
00369
00370
00371
00372
00373
00374
00375 G4DynamicParticle* aGamma= new G4DynamicParticle ( G4Gamma::Gamma(),
00376 gammaDirection,
00377 energyOfSR );
00378 aGamma->SetPolarization( gammaPolarization.x(),
00379 gammaPolarization.y(),
00380 gammaPolarization.z() );
00381
00382
00383 aParticleChange.SetNumberOfSecondaries(1);
00384 aParticleChange.AddSecondary(aGamma);
00385
00386
00387
00388 G4double newKinEnergy = kineticEnergy - energyOfSR ;
00389
00390 if (newKinEnergy > 0.)
00391 {
00392 aParticleChange.ProposeMomentumDirection( particleDirection );
00393 aParticleChange.ProposeEnergy( newKinEnergy );
00394 aParticleChange.ProposeLocalEnergyDeposit (0.);
00395 }
00396 else
00397 {
00398 aParticleChange.ProposeEnergy( 0. );
00399 aParticleChange.ProposeLocalEnergyDeposit (0.);
00400 G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00401 if (charge<0.)
00402 {
00403 aParticleChange.ProposeTrackStatus(fStopAndKill) ;
00404 }
00405 else
00406 {
00407 aParticleChange.ProposeTrackStatus(fStopButAlive) ;
00408 }
00409 }
00410 }
00411 else
00412 {
00413 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00414 }
00415 }
00416 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00417 }
00418
00419
00420 G4double
00421 G4SynchrotronRadiationInMat::GetPhotonEnergy( const G4Track& trackData,
00422 const G4Step& )
00423
00424 {
00425 G4int i ;
00426 G4double energyOfSR = -1.0 ;
00427
00428
00429 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
00430
00431 G4double gamma = aDynamicParticle->GetTotalEnergy()/
00432 (aDynamicParticle->GetMass() ) ;
00433
00434 G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00435
00436 G4ThreeVector FieldValue;
00437 const G4Field* pField = 0 ;
00438
00439 G4FieldManager* fieldMgr=0;
00440 G4bool fieldExertsForce = false;
00441
00442 if( (particleCharge != 0.0) )
00443 {
00444 fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
00445 if ( fieldMgr != 0 )
00446 {
00447
00448
00449 fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
00450 }
00451 }
00452 if ( fieldExertsForce )
00453 {
00454 pField = fieldMgr->GetDetectorField();
00455 G4ThreeVector globPosition = trackData.GetPosition();
00456 G4double globPosVec[3], FieldValueVec[3];
00457
00458 globPosVec[0] = globPosition.x();
00459 globPosVec[1] = globPosition.y();
00460 globPosVec[2] = globPosition.z();
00461
00462 pField->GetFieldValue( globPosVec, FieldValueVec );
00463 FieldValue = G4ThreeVector( FieldValueVec[0],
00464 FieldValueVec[1],
00465 FieldValueVec[2] );
00466
00467 G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
00468 G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum) ;
00469 G4double perpB = unitMcrossB.mag();
00470 if( perpB > 0.0 )
00471 {
00472
00473
00474 G4double random = G4UniformRand() ;
00475 for(i=0;i<200;i++)
00476 {
00477 if(random >= fIntegralProbabilityOfSR[i]) break ;
00478 }
00479 energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
00480
00481
00482
00483 if(energyOfSR <= 0.0)
00484 {
00485 return -1.0 ;
00486 }
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501 }
00502 else
00503 {
00504 return -1.0 ;
00505 }
00506 }
00507 return energyOfSR ;
00508 }
00509
00511
00512
00513
00514 G4double G4SynchrotronRadiationInMat::GetRandomEnergySR(G4double gamma, G4double perpB)
00515 {
00516 G4int i, iMax;
00517 G4double energySR, random, position;
00518
00519 iMax = 200;
00520 random = G4UniformRand();
00521
00522 for( i = 0; i < iMax; i++ )
00523 {
00524 if( random >= fIntegralProbabilityOfSR[i] ) break;
00525 }
00526 if(i <= 0 ) position = G4UniformRand();
00527 else if( i>= iMax) position = G4double(iMax);
00528 else position = i + G4UniformRand();
00529
00530
00531
00532
00533 energySR = 0.0001*position*position*fEnergyConst*gamma*gamma*perpB;
00534
00535 if( energySR < 0. ) energySR = 0.;
00536
00537 return energySR;
00538 }
00539
00541
00542
00543
00544 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt( G4double t)
00545 {
00546 G4double result, hypCos2, hypCos=std::cosh(t);
00547
00548 hypCos2 = hypCos*hypCos;
00549 result = std::cosh(5.*t/3.)*std::exp(t-fKsi*hypCos);
00550 result /= hypCos2;
00551 return result;
00552 }
00553
00555
00556
00557
00558
00559
00560 G4double G4SynchrotronRadiationInMat::GetIntProbSR( G4double ksi)
00561 {
00562 if (ksi <= 0.) return 1.0;
00563 fKsi = ksi;
00564 G4int n;
00565 G4double result, a;
00566
00567 a = fAlpha;
00568 n = fRootNumber;
00569
00570 G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
00571
00572 result = integral.Laguerre(this,
00573 &G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt, a, n);
00574
00575 result *= 3./5./pi;
00576
00577 return result;
00578 }
00579
00581
00582
00583
00584 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy( G4double t)
00585 {
00586 G4double result, hypCos=std::cosh(t);
00587
00588 result = std::cosh(5.*t/3.)*std::exp(t - fKsi*hypCos);
00589 result /= hypCos;
00590 return result;
00591 }
00592
00594
00595
00596
00597
00598
00599 G4double G4SynchrotronRadiationInMat::GetEnergyProbSR( G4double ksi)
00600 {
00601 if (ksi <= 0.) return 1.0;
00602 fKsi = ksi;
00603 G4int n;
00604 G4double result, a;
00605
00606 a = fAlpha;
00607 n = fRootNumber;
00608
00609 G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
00610
00611 result = integral.Laguerre(this,
00612 &G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy, a, n);
00613
00614 result *= 9.*std::sqrt(3.)*ksi/8./pi;
00615
00616 return result;
00617 }
00618
00620
00621
00622
00623 G4double G4SynchrotronRadiationInMat::GetIntegrandForAngleK( G4double t)
00624 {
00625 G4double result, hypCos=std::cosh(t);
00626
00627 result = std::cosh(fOrderAngleK*t)*std::exp(t - fEta*hypCos);
00628 result /= hypCos;
00629 return result;
00630 }
00631
00633
00634
00635
00636 G4double G4SynchrotronRadiationInMat::GetAngleK( G4double eta)
00637 {
00638 fEta = eta;
00639 G4int n;
00640 G4double result, a;
00641
00642 a = fAlpha;
00643 n = fRootNumber;
00644
00645 G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
00646
00647 result = integral.Laguerre(this,
00648 &G4SynchrotronRadiationInMat::GetIntegrandForAngleK, a, n);
00649
00650 return result;
00651 }
00652
00654
00655
00656
00657 G4double G4SynchrotronRadiationInMat::GetAngleNumberAtGammaKsi( G4double gpsi)
00658 {
00659 G4double result, funK, funK2, gpsi2 = gpsi*gpsi;
00660
00661 fPsiGamma = gpsi;
00662 fEta = 0.5*fKsi*(1 + gpsi2)*std::sqrt(1 + gpsi2);
00663
00664 fOrderAngleK = 1./3.;
00665 funK = GetAngleK(fEta);
00666 funK2 = funK*funK;
00667
00668 result = gpsi2*funK2/(1 + gpsi2);
00669
00670 fOrderAngleK = 2./3.;
00671 funK = GetAngleK(fEta);
00672 funK2 = funK*funK;
00673
00674 result += funK2;
00675 result *= (1 + gpsi2)*fKsi;
00676
00677 return result;
00678 }
00679
00680
00682