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
00043
00044 #include "G4SynchrotronRadiation.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4UnitsTable.hh"
00048 #include "G4EmProcessSubType.hh"
00049
00051
00052
00053
00054
00055 G4SynchrotronRadiation::G4SynchrotronRadiation(const G4String& processName,
00056 G4ProcessType type):G4VDiscreteProcess (processName, type),
00057 theGamma (G4Gamma::Gamma() ),
00058 theElectron ( G4Electron::Electron() ),
00059 thePositron ( G4Positron::Positron() )
00060 {
00061 G4TransportationManager* transportMgr =
00062 G4TransportationManager::GetTransportationManager();
00063
00064 fFieldPropagator = transportMgr->GetPropagatorInField();
00065
00066 fLambdaConst = std::sqrt(3.0)*electron_mass_c2/
00067 (2.5*fine_structure_const*eplus*c_light);
00068 fEnergyConst = 1.5*c_light*c_light*eplus*hbar_Planck/electron_mass_c2 ;
00069
00070 SetProcessSubType(fSynchrotronRadiation);
00071 verboseLevel=1;
00072 }
00073
00075
00076
00077
00078
00079 G4SynchrotronRadiation::~G4SynchrotronRadiation()
00080 {}
00081
00083
00084
00085
00086
00087
00088
00089
00090 G4double
00091 G4SynchrotronRadiation::GetMeanFreePath( const G4Track& trackData,
00092 G4double,
00093 G4ForceCondition* condition)
00094 {
00095
00096 G4double MeanFreePath;
00097
00098 const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
00099
00100 *condition = NotForced;
00101
00102 G4double gamma = aDynamicParticle->GetTotalEnergy()/
00103 aDynamicParticle->GetMass();
00104
00105 G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00106
00107 if ( gamma < 1.0e3 ) MeanFreePath = DBL_MAX;
00108 else
00109 {
00110
00111 G4ThreeVector FieldValue;
00112 const G4Field* pField = 0;
00113
00114 G4FieldManager* fieldMgr=0;
00115 G4bool fieldExertsForce = false;
00116
00117 if( (particleCharge != 0.0) )
00118 {
00119 fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
00120
00121 if ( fieldMgr != 0 )
00122 {
00123
00124
00125 fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
00126 }
00127 }
00128 if ( fieldExertsForce )
00129 {
00130 pField = fieldMgr->GetDetectorField();
00131 G4ThreeVector globPosition = trackData.GetPosition();
00132
00133 G4double globPosVec[4], FieldValueVec[6];
00134
00135 globPosVec[0] = globPosition.x();
00136 globPosVec[1] = globPosition.y();
00137 globPosVec[2] = globPosition.z();
00138 globPosVec[3] = trackData.GetGlobalTime();
00139
00140 pField->GetFieldValue( globPosVec, FieldValueVec );
00141
00142 FieldValue = G4ThreeVector( FieldValueVec[0],
00143 FieldValueVec[1],
00144 FieldValueVec[2] );
00145
00146
00147
00148 G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
00149 G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
00150 G4double perpB = unitMcrossB.mag();
00151
00152 if( perpB > 0.0 ) MeanFreePath = fLambdaConst/perpB;
00153 else MeanFreePath = DBL_MAX;
00154
00155 static G4bool FirstTime=true;
00156 if(verboseLevel > 0 && FirstTime)
00157 {
00158 G4cout << "G4SynchrotronRadiation::GetMeanFreePath :" << '\n'
00159 << " MeanFreePath = " << G4BestUnit(MeanFreePath, "Length")
00160 << G4endl;
00161 if(verboseLevel > 1)
00162 {
00163 G4ThreeVector pvec=aDynamicParticle->GetMomentum();
00164 G4double Btot=FieldValue.getR();
00165 G4double ptot=pvec.getR();
00166 G4double rho= ptot / (MeV * c_light * Btot );
00167 G4double Theta=unitMomentum.theta(FieldValue);
00168 G4cout
00169 << " B = " << Btot/tesla << " Tesla"
00170 << " perpB = " << perpB/tesla << " Tesla"
00171 << " Theta = " << Theta << " std::sin(Theta)=" << std::sin(Theta) << '\n'
00172 << " ptot = " << G4BestUnit(ptot,"Energy")
00173 << " rho = " << G4BestUnit(rho,"Length")
00174 << G4endl;
00175 }
00176 FirstTime=false;
00177 }
00178 }
00179 else MeanFreePath = DBL_MAX;
00180
00181
00182 }
00183
00184 return MeanFreePath;
00185 }
00186
00188
00189
00190
00191 G4VParticleChange*
00192 G4SynchrotronRadiation::PostStepDoIt(const G4Track& trackData,
00193 const G4Step& stepData )
00194
00195 {
00196 aParticleChange.Initialize(trackData);
00197
00198 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
00199
00200 G4double gamma = aDynamicParticle->GetTotalEnergy()/
00201 (aDynamicParticle->GetMass() );
00202
00203 if(gamma <= 1.0e3 )
00204 {
00205 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00206 }
00207 G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
00208
00209 G4ThreeVector FieldValue;
00210 const G4Field* pField = 0;
00211
00212 G4FieldManager* fieldMgr=0;
00213 G4bool fieldExertsForce = false;
00214
00215 if( (particleCharge != 0.0) )
00216 {
00217 fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
00218 if ( fieldMgr != 0 )
00219 {
00220
00221
00222 fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
00223 }
00224 }
00225 if ( fieldExertsForce )
00226 {
00227 pField = fieldMgr->GetDetectorField();
00228 G4ThreeVector globPosition = trackData.GetPosition();
00229 G4double globPosVec[4], FieldValueVec[6];
00230 globPosVec[0] = globPosition.x();
00231 globPosVec[1] = globPosition.y();
00232 globPosVec[2] = globPosition.z();
00233 globPosVec[3] = trackData.GetGlobalTime();
00234
00235 pField->GetFieldValue( globPosVec, FieldValueVec );
00236 FieldValue = G4ThreeVector( FieldValueVec[0],
00237 FieldValueVec[1],
00238 FieldValueVec[2] );
00239
00240 G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
00241 G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
00242 G4double perpB = unitMcrossB.mag();
00243 if(perpB > 0.0)
00244 {
00245
00246
00247 G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
00248
00249
00250
00251 if( energyOfSR <= 0.0 )
00252 {
00253 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00254 }
00255 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00256 G4ParticleMomentum
00257 particleDirection = aDynamicParticle->GetMomentumDirection();
00258
00259
00260
00261
00262
00263 G4double cosTheta, sinTheta, fcos, beta;
00264
00265 do
00266 {
00267 cosTheta = 1. - 2.*G4UniformRand();
00268 fcos = (1 + cosTheta*cosTheta)*0.5;
00269 }
00270 while( fcos < G4UniformRand() );
00271
00272 beta = std::sqrt(1. - 1./(gamma*gamma));
00273
00274 cosTheta = (cosTheta + beta)/(1. + beta*cosTheta);
00275
00276 if( cosTheta > 1. ) cosTheta = 1.;
00277 if( cosTheta < -1. ) cosTheta = -1.;
00278
00279 sinTheta = std::sqrt(1. - cosTheta*cosTheta );
00280
00281 G4double Phi = twopi * G4UniformRand();
00282
00283 G4double dirx = sinTheta*std::cos(Phi) ,
00284 diry = sinTheta*std::sin(Phi) ,
00285 dirz = cosTheta;
00286
00287 G4ThreeVector gammaDirection ( dirx, diry, dirz);
00288 gammaDirection.rotateUz(particleDirection);
00289
00290
00291
00292
00293
00294
00295
00296 G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
00297 gammaPolarization = gammaPolarization.unit();
00298
00299
00300
00301
00302
00303
00304 G4DynamicParticle* aGamma= new G4DynamicParticle ( theGamma,
00305 gammaDirection,
00306 energyOfSR );
00307 aGamma->SetPolarization( gammaPolarization.x(),
00308 gammaPolarization.y(),
00309 gammaPolarization.z() );
00310
00311
00312 aParticleChange.SetNumberOfSecondaries(1);
00313 aParticleChange.AddSecondary(aGamma);
00314
00315
00316
00317 G4double newKinEnergy = kineticEnergy - energyOfSR;
00318 aParticleChange.ProposeLocalEnergyDeposit (0.);
00319
00320 if (newKinEnergy > 0.)
00321 {
00322 aParticleChange.ProposeMomentumDirection( particleDirection );
00323 aParticleChange.ProposeEnergy( newKinEnergy );
00324 }
00325 else
00326 {
00327 aParticleChange.ProposeEnergy( 0. );
00328 }
00329 }
00330 }
00331 return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
00332 }
00333
00334
00336
00337
00338
00339 G4double G4SynchrotronRadiation::InvSynFracInt(G4double x)
00340
00341 {
00342
00343 const G4double aa1=0 ,aa2=0.7;
00344 const G4int ncheb1=27;
00345 static const G4double cheb1[] =
00346 { 1.22371665676046468821,0.108956475422163837267,0.0383328524358594396134,0.00759138369340257753721,
00347 0.00205712048644963340914,0.000497810783280019308661,0.000130743691810302187818,0.0000338168760220395409734,
00348 8.97049680900520817728e-6,2.38685472794452241466e-6,6.41923109149104165049e-7,1.73549898982749277843e-7,
00349 4.72145949240790029153e-8,1.29039866111999149636e-8,3.5422080787089834182e-9,9.7594757336403784905e-10,
00350 2.6979510184976065731e-10,7.480422622550977077e-11,2.079598176402699913e-11,5.79533622220841193e-12,
00351 1.61856011449276096e-12,4.529450993473807e-13,1.2698603951096606e-13,3.566117394511206e-14,1.00301587494091e-14,
00352 2.82515346447219e-15,7.9680747949792e-16};
00353
00354 const G4double aa3=0.9132260271183847;
00355 const G4int ncheb2=27;
00356 static const G4double cheb2[] =
00357 { 1.1139496701107756,0.3523967429328067,0.0713849171926623,0.01475818043595387,0.003381255637322462,
00358 0.0008228057599452224,0.00020785506681254216,0.00005390169253706556,0.000014250571923902464,3.823880733161044e-6,
00359 1.0381966089136036e-6,2.8457557457837253e-7,7.86223332179956e-8,2.1866609342508474e-8,6.116186259857143e-9,
00360 1.7191233618437565e-9,4.852755117740807e-10,1.3749966961763457e-10,3.908961987062447e-11,1.1146253766895824e-11,
00361 3.1868887323415814e-12,9.134319791300977e-13,2.6211077371181566e-13,7.588643377757906e-14,2.1528376972619e-14,
00362 6.030906040404772e-15,1.9549163926819867e-15};
00363
00364
00365 const G4double aa4=2.4444485538746025480,aa5=9.3830728608909477079;
00366 const G4int ncheb3=28;
00367 static const G4double cheb3[] =
00368 { 1.2292683840435586977,0.160353449247864455879,-0.0353559911947559448721,0.00776901561223573936985,
00369 -0.00165886451971685133259,0.000335719118906954279467,-0.0000617184951079161143187,9.23534039743246708256e-6,
00370 -6.06747198795168022842e-7,-3.07934045961999778094e-7,1.98818772614682367781e-7,-8.13909971567720135413e-8,
00371 2.84298174969641838618e-8,-9.12829766621316063548e-9,2.77713868004820551077e-9,-8.13032767247834023165e-10,
00372 2.31128525568385247392e-10,-6.41796873254200220876e-11,1.74815310473323361543e-11,-4.68653536933392363045e-12,
00373 1.24016595805520752748e-12,-3.24839432979935522159e-13,8.44601465226513952994e-14,-2.18647276044246803998e-14,
00374 5.65407548745690689978e-15,-1.46553625917463067508e-15,3.82059606377570462276e-16,-1.00457896653436912508e-16};
00375 const G4double aa6=33.122936966163038145;
00376 const G4int ncheb4=27;
00377 static const G4double cheb4[] =
00378 {1.69342658227676741765,0.0742766400841232319225,-0.019337880608635717358,0.00516065527473364110491,
00379 -0.00139342012990307729473,0.000378549864052022522193,-0.000103167085583785340215,0.0000281543441271412178337,
00380 -7.68409742018258198651e-6,2.09543221890204537392e-6,-5.70493140367526282946e-7,1.54961164548564906446e-7,
00381 -4.19665599629607704794e-8,1.13239680054166507038e-8,-3.04223563379021441863e-9,8.13073745977562957997e-10,
00382 -2.15969415476814981374e-10,5.69472105972525594811e-11,-1.48844799572430829499e-11,3.84901514438304484973e-12,
00383 -9.82222575944247161834e-13,2.46468329208292208183e-13,-6.04953826265982691612e-14,1.44055805710671611984e-14,
00384 -3.28200813577388740722e-15,6.96566359173765367675e-16,-1.294122794852896275e-16};
00385
00386 if(x<aa2) return x*x*x*Chebyshev(aa1,aa2,cheb1,ncheb1,x);
00387 else if(x<aa3) return Chebyshev(aa2,aa3,cheb2,ncheb2,x);
00388 else if(x<1-0.0000841363)
00389 { G4double y=-std::log(1-x);
00390 return y*Chebyshev(aa4,aa5,cheb3,ncheb3,y);
00391 }
00392 else
00393 { G4double y=-std::log(1-x);
00394 return y*Chebyshev(aa5,aa6,cheb4,ncheb4,y);
00395 }
00396 }
00397
00398 G4double G4SynchrotronRadiation::GetRandomEnergySR(G4double gamma, G4double perpB)
00399 {
00400
00401 G4double Ecr=fEnergyConst*gamma*gamma*perpB;
00402
00403 static G4bool FirstTime=true;
00404 if(verboseLevel > 0 && FirstTime)
00405 { G4double Emean=8./(15.*std::sqrt(3.))*Ecr;
00406 G4double E_rms=std::sqrt(211./675.)*Ecr;
00407 G4int prec = G4cout.precision();
00408 G4cout << "G4SynchrotronRadiation::GetRandomEnergySR :" << '\n' << std::setprecision(4)
00409 << " Ecr = " << G4BestUnit(Ecr,"Energy") << '\n'
00410 << " Emean = " << G4BestUnit(Emean,"Energy") << '\n'
00411 << " E_rms = " << G4BestUnit(E_rms,"Energy") << G4endl;
00412 FirstTime=false;
00413 G4cout.precision(prec);
00414 }
00415
00416 G4double energySR=Ecr*InvSynFracInt(G4UniformRand());
00417 return energySR;
00418 }
00419
00420
00421 void G4SynchrotronRadiation::BuildPhysicsTable(const G4ParticleDefinition& part)
00422 {
00423 if(0 < verboseLevel && &part==theElectron ) PrintInfoDefinition();
00424 }
00425
00426 void G4SynchrotronRadiation::PrintInfoDefinition()
00427 {
00428 G4String comments ="Incoherent Synchrotron Radiation\n";
00429 G4cout << G4endl << GetProcessName() << ": " << comments
00430 << " good description for long magnets at all energies" << G4endl;
00431 }
00432