G4SynchrotronRadiation Class Reference

#include <G4SynchrotronRadiation.hh>

Inheritance diagram for G4SynchrotronRadiation:

G4VDiscreteProcess G4VProcess

Public Member Functions

 G4SynchrotronRadiation (const G4String &pName="SynRad", G4ProcessType type=fElectromagnetic)
virtual ~G4SynchrotronRadiation ()
G4double GetMeanFreePath (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4VParticleChangePostStepDoIt (const G4Track &track, const G4Step &Step)
G4double GetPhotonEnergy (const G4Track &trackData, const G4Step &stepData)
G4double GetRandomEnergySR (G4double, G4double)
G4double InvSynFracInt (G4double x)
G4double Chebyshev (G4double a, G4double b, const G4double c[], G4int n, G4double x)
G4bool IsApplicable (const G4ParticleDefinition &)
void BuildPhysicsTable (const G4ParticleDefinition &)
void PrintInfoDefinition ()

Detailed Description

Definition at line 63 of file G4SynchrotronRadiation.hh.


Constructor & Destructor Documentation

G4SynchrotronRadiation::G4SynchrotronRadiation ( const G4String pName = "SynRad",
G4ProcessType  type = fElectromagnetic 
)

Definition at line 55 of file G4SynchrotronRadiation.cc.

References fSynchrotronRadiation, G4TransportationManager::GetPropagatorInField(), G4TransportationManager::GetTransportationManager(), G4VProcess::SetProcessSubType(), and G4VProcess::verboseLevel.

00056                      :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 }

G4SynchrotronRadiation::~G4SynchrotronRadiation (  )  [virtual]

Definition at line 79 of file G4SynchrotronRadiation.cc.

00080 {}


Member Function Documentation

void G4SynchrotronRadiation::BuildPhysicsTable ( const G4ParticleDefinition  )  [virtual]

Reimplemented from G4VProcess.

Definition at line 421 of file G4SynchrotronRadiation.cc.

References PrintInfoDefinition(), and G4VProcess::verboseLevel.

00422 {
00423   if(0 < verboseLevel && &part==theElectron ) PrintInfoDefinition();
00424 }

G4double G4SynchrotronRadiation::Chebyshev ( G4double  a,
G4double  b,
const G4double  c[],
G4int  n,
G4double  x 
) [inline]

Definition at line 117 of file G4SynchrotronRadiation.hh.

Referenced by InvSynFracInt().

00119 {
00120   G4double y;
00121   G4double y2=2.0*(y=(2.0*x-a-b)/(b-a)); // Change of variable.
00122   G4double d=0,dd=0;
00123   for (G4int j=n-1;j>=1;--j) // Clenshaw's recurrence.
00124   { G4double sv=d;
00125         d=y2*d-dd+c[j];
00126         dd=sv;
00127   }
00128   return y*d-dd+0.5*c[0];
00129 }

G4double G4SynchrotronRadiation::GetMeanFreePath ( const G4Track track,
G4double  previousStepSize,
G4ForceCondition condition 
) [virtual]

Implements G4VDiscreteProcess.

Definition at line 91 of file G4SynchrotronRadiation.cc.

References DBL_MAX, G4PropagatorInField::FindAndSetFieldManager(), G4BestUnit, G4cout, G4endl, G4DynamicParticle::GetDefinition(), G4FieldManager::GetDetectorField(), G4Track::GetDynamicParticle(), G4Field::GetFieldValue(), G4Track::GetGlobalTime(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGCharge(), G4Track::GetPosition(), G4DynamicParticle::GetTotalEnergy(), G4Track::GetVolume(), NotForced, and G4VProcess::verboseLevel.

00094 {
00095   // gives the MeanFreePath in GEANT4 internal units
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         // If the field manager has no field, there is no field !
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 ); // full bending radius
00167           G4double Theta=unitMomentum.theta(FieldValue); // angle between particle and field
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 }

G4double G4SynchrotronRadiation::GetPhotonEnergy ( const G4Track trackData,
const G4Step stepData 
)

G4double G4SynchrotronRadiation::GetRandomEnergySR ( G4double  ,
G4double   
)

Definition at line 398 of file G4SynchrotronRadiation.cc.

References G4BestUnit, G4cout, G4endl, G4UniformRand, InvSynFracInt(), and G4VProcess::verboseLevel.

Referenced by PostStepDoIt().

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; // mean photon energy
00406     G4double E_rms=std::sqrt(211./675.)*Ecr; // rms of photon energy distribution
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 }

G4double G4SynchrotronRadiation::InvSynFracInt ( G4double  x  ) 

Definition at line 339 of file G4SynchrotronRadiation.cc.

References Chebyshev().

Referenced by GetRandomEnergySR().

00341 {
00342   // from 0 to 0.7
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   //   from 0.7 to 0.9132260271183847
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    // Chebyshev with exp/log  scale
00364    // a = -Log[1 - SynFracInt[1]]; b = -Log[1 - SynFracInt[7]];
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 }

G4bool G4SynchrotronRadiation::IsApplicable ( const G4ParticleDefinition  )  [inline, virtual]

Reimplemented from G4VProcess.

Definition at line 109 of file G4SynchrotronRadiation.hh.

00110 {
00111 
00112   return ( ( &particle == (const G4ParticleDefinition *)theElectron ) ||
00113            ( &particle == (const G4ParticleDefinition *)thePositron )    );
00114 }

G4VParticleChange * G4SynchrotronRadiation::PostStepDoIt ( const G4Track track,
const G4Step Step 
) [virtual]

Reimplemented from G4VDiscreteProcess.

Definition at line 192 of file G4SynchrotronRadiation.cc.

References G4ParticleChange::AddSecondary(), G4VProcess::aParticleChange, fcos(), G4PropagatorInField::FindAndSetFieldManager(), G4UniformRand, G4DynamicParticle::GetDefinition(), G4FieldManager::GetDetectorField(), G4Track::GetDynamicParticle(), G4Field::GetFieldValue(), G4Track::GetGlobalTime(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGCharge(), G4Track::GetPosition(), GetRandomEnergySR(), G4DynamicParticle::GetTotalEnergy(), G4Track::GetVolume(), G4ParticleChange::Initialize(), G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::SetNumberOfSecondaries(), and G4DynamicParticle::SetPolarization().

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       // If the field manager has no field, there is no field !
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       // M-C of synchrotron photon energy
00246 
00247       G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
00248 
00249       // check against insufficient energy
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       // M-C of its direction, simplified dipole boosted approach
00260 
00261       // G4double Teta, fteta;  // = G4UniformRand()/gamma;    // Very roughly
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       // polarization of new gamma
00291 
00292       // G4double sx =  std::cos(Teta)*std::cos(Phi);
00293       // G4double sy =  std::cos(Teta)*std::sin(Phi);
00294       // G4double sz = -std::sin(Teta);
00295 
00296       G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
00297       gammaPolarization = gammaPolarization.unit();
00298 
00299       // (sx, sy, sz);
00300       // gammaPolarization.rotateUz(particleDirection);
00301 
00302       // create G4DynamicParticle object for the SR photon
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       // Update the incident particle
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 }

void G4SynchrotronRadiation::PrintInfoDefinition (  ) 

Definition at line 426 of file G4SynchrotronRadiation.cc.

References G4cout, G4endl, and G4VProcess::GetProcessName().

Referenced by BuildPhysicsTable().

00427 {
00428   G4String comments ="Incoherent Synchrotron Radiation\n";
00429   G4cout << G4endl << GetProcessName() << ":  " << comments
00430          << "      good description for long magnets at all energies" << G4endl;
00431 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:28 2013 for Geant4 by  doxygen 1.4.7