Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/xrays/src/G4SynchrotronRadiation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // --------------------------------------------------------------
 27 //      GEANT 4 class implementation file
 28 //
 29 //      History: first implementation,
 30 //      21-5-98 V.Grichine
 31 //      28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation
 32 //      04.03.05, V.Grichine: get local field interface
 33 //      18-05-06 H. Burkhardt: Energy spectrum from function rather than table
 34 //
 35 ///////////////////////////////////////////////////////////////////////////
 36 
 37 #include "G4SynchrotronRadiation.hh"
 38 
 39 #include "G4DipBustGenerator.hh"
 40 #include "G4Electron.hh"
 41 #include "G4EmProcessSubType.hh"
 42 #include "G4Log.hh"
 43 #include "G4LossTableManager.hh"
 44 #include "G4Gamma.hh"
 45 #include "G4PhysicalConstants.hh"
 46 #include "G4PropagatorInField.hh"
 47 #include "G4SystemOfUnits.hh"
 48 #include "G4TransportationManager.hh"
 49 #include "G4UnitsTable.hh"
 50 #include "G4PhysicsModelCatalog.hh"
 51 
 52 ///////////////////////////////////////////////////////////////////////
 53 //  Constructor
 54 G4SynchrotronRadiation::G4SynchrotronRadiation(const G4String& processName,
 55                                                G4ProcessType type)
 56   : G4VDiscreteProcess(processName, type)
 57   , theGamma(G4Gamma::Gamma())
 58 {
 59   G4TransportationManager* transportMgr =
 60     G4TransportationManager::GetTransportationManager();
 61 
 62   fFieldPropagator = transportMgr->GetPropagatorInField();
 63 
 64   secID = G4PhysicsModelCatalog::GetModelID("model_SynRad");
 65   SetProcessSubType(fSynchrotronRadiation);
 66   verboseLevel = 1;
 67   FirstTime    = true;
 68   FirstTime1   = true;
 69   genAngle     = nullptr;
 70   SetAngularGenerator(new G4DipBustGenerator());
 71   theManager = G4LossTableManager::Instance();
 72   theManager->Register(this);
 73 }
 74 
 75 /////////////////////////////////////////////////////////////////////////
 76 // Destructor
 77 G4SynchrotronRadiation::~G4SynchrotronRadiation()
 78 {
 79   delete genAngle;
 80   theManager->DeRegister(this);
 81 }
 82 
 83 /////////////////////////////// METHODS /////////////////////////////////
 84 
 85 void G4SynchrotronRadiation::SetAngularGenerator(G4VEmAngularDistribution* p)
 86 {
 87   if(p != genAngle)
 88   {
 89     delete genAngle;
 90     genAngle = p;
 91   }
 92 }
 93 
 94 G4bool G4SynchrotronRadiation::IsApplicable(
 95   const G4ParticleDefinition& particle)
 96 {
 97   return (particle.GetPDGCharge() != 0.0 && !particle.IsShortLived());
 98 }
 99 
100 /////////////////////////////////////////////////////////////////////////
101 // Production of synchrotron X-ray photon
102 // Geant4 internal units.
103 G4double G4SynchrotronRadiation::GetMeanFreePath(const G4Track& trackData,
104                                                  G4double,
105                                                  G4ForceCondition* condition)
106 {
107   // gives the MeanFreePath in Geant4 internal units
108   G4double MeanFreePath = DBL_MAX;
109 
110   const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
111 
112   *condition = NotForced;
113 
114   G4double gamma =
115     aDynamicParticle->GetTotalEnergy() / aDynamicParticle->GetMass();
116 
117   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
118 
119   if(gamma < 1.0e3 || 0.0 == particleCharge)
120   {
121     MeanFreePath = DBL_MAX;
122   }
123   else
124   {
125     G4ThreeVector FieldValue;
126     const G4Field* pField   = nullptr;
127     G4bool fieldExertsForce = false;
128 
129     G4FieldManager* fieldMgr =
130       fFieldPropagator->FindAndSetFieldManager(trackData.GetVolume());
131 
132     if(fieldMgr != nullptr)
133     {
134       // If the field manager has no field, there is no field !
135       fieldExertsForce = (fieldMgr->GetDetectorField() != nullptr);
136     }
137 
138     if(fieldExertsForce)
139     {
140       pField                     = fieldMgr->GetDetectorField();
141       G4ThreeVector globPosition = trackData.GetPosition();
142 
143       G4double globPosVec[4], FieldValueVec[6];
144 
145       globPosVec[0] = globPosition.x();
146       globPosVec[1] = globPosition.y();
147       globPosVec[2] = globPosition.z();
148       globPosVec[3] = trackData.GetGlobalTime();
149 
150       pField->GetFieldValue(globPosVec, FieldValueVec);
151 
152       FieldValue =
153         G4ThreeVector(FieldValueVec[0], FieldValueVec[1], FieldValueVec[2]);
154 
155       G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
156       G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum);
157       G4double perpB             = unitMcrossB.mag();
158 
159       static const G4double fLambdaConst =
160         std::sqrt(3.0) * eplus / (2.5 * fine_structure_const * c_light);
161 
162       if(perpB > 0.0)
163       {
164         MeanFreePath = fLambdaConst *
165                        aDynamicParticle->GetDefinition()->GetPDGMass() /
166                        (perpB * particleCharge * particleCharge);
167       }
168       if(verboseLevel > 0 && FirstTime)
169       {
170         G4cout << "G4SynchrotronRadiation::GetMeanFreePath "
171                << " for particle "
172                << aDynamicParticle->GetDefinition()->GetParticleName() << ":"
173                << '\n'
174                << "  MeanFreePath = " << G4BestUnit(MeanFreePath, "Length")
175                << G4endl;
176         if(verboseLevel > 1)
177         {
178           G4ThreeVector pvec = aDynamicParticle->GetMomentum();
179           G4double Btot      = FieldValue.getR();
180           G4double ptot      = pvec.getR();
181           G4double rho       = ptot / (MeV * c_light * Btot);
182           // full bending radius
183           G4double Theta = unitMomentum.theta(FieldValue);
184           // angle between particle and field
185           G4cout << "  B = " << Btot / tesla << " Tesla"
186                  << "  perpB = " << perpB / tesla << " Tesla"
187                  << "  Theta = " << Theta
188                  << " std::sin(Theta)=" << std::sin(Theta) << '\n'
189                  << "  ptot  = " << G4BestUnit(ptot, "Energy")
190                  << "  rho   = " << G4BestUnit(rho, "Length") << G4endl;
191         }
192         FirstTime = false;
193       }
194     }
195   }
196   return MeanFreePath;
197 }
198 
199 ///////////////////////////////////////////////////////////////////////////////
200 G4VParticleChange* G4SynchrotronRadiation::PostStepDoIt(
201   const G4Track& trackData, const G4Step& stepData)
202 
203 {
204   aParticleChange.Initialize(trackData);
205 
206   const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
207 
208   G4double gamma = aDynamicParticle->GetTotalEnergy() /
209                    (aDynamicParticle->GetDefinition()->GetPDGMass());
210 
211   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
212   if(gamma <= 1.0e3 || 0.0 == particleCharge)
213   {
214     return G4VDiscreteProcess::PostStepDoIt(trackData, stepData);
215   }
216 
217   G4ThreeVector FieldValue;
218   const G4Field* pField = nullptr;
219 
220   G4bool fieldExertsForce = false;
221   G4FieldManager* fieldMgr =
222     fFieldPropagator->FindAndSetFieldManager(trackData.GetVolume());
223 
224   if(fieldMgr != nullptr)
225   {
226     // If the field manager has no field, there is no field !
227     fieldExertsForce = (fieldMgr->GetDetectorField() != nullptr);
228   }
229 
230   if(fieldExertsForce)
231   {
232     pField                     = fieldMgr->GetDetectorField();
233     G4ThreeVector globPosition = trackData.GetPosition();
234     G4double globPosVec[4], FieldValueVec[6];
235     globPosVec[0] = globPosition.x();
236     globPosVec[1] = globPosition.y();
237     globPosVec[2] = globPosition.z();
238     globPosVec[3] = trackData.GetGlobalTime();
239 
240     pField->GetFieldValue(globPosVec, FieldValueVec);
241     FieldValue =
242       G4ThreeVector(FieldValueVec[0], FieldValueVec[1], FieldValueVec[2]);
243 
244     G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
245     G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum);
246     G4double perpB             = unitMcrossB.mag();
247     if(perpB > 0.0)
248     {
249       // M-C of synchrotron photon energy
250       G4double energyOfSR = GetRandomEnergySR(
251         gamma, perpB, aDynamicParticle->GetDefinition()->GetPDGMass());
252 
253       // check against insufficient energy
254       if(energyOfSR <= 0.0)
255       {
256         return G4VDiscreteProcess::PostStepDoIt(trackData, stepData);
257       }
258       G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
259       G4ThreeVector gammaDirection =
260         genAngle->SampleDirection(aDynamicParticle, energyOfSR, 1, nullptr);
261 
262       G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
263       gammaPolarization               = gammaPolarization.unit();
264 
265       // create G4DynamicParticle object for the SR photon
266       auto aGamma =
267         new G4DynamicParticle(theGamma, gammaDirection, energyOfSR);
268       aGamma->SetPolarization(gammaPolarization.x(), gammaPolarization.y(),
269                               gammaPolarization.z());
270 
271       aParticleChange.SetNumberOfSecondaries(1);
272 
273       // Update the incident particle
274       G4double newKinEnergy = kineticEnergy - energyOfSR;
275 
276       if(newKinEnergy > 0.)
277       {
278         aParticleChange.ProposeEnergy(newKinEnergy);
279       }
280       else
281       {
282         aParticleChange.ProposeEnergy(0.);
283       }
284 
285       // Create the G4Track
286       G4Track* aSecondaryTrack = new G4Track(aGamma, trackData.GetGlobalTime(), trackData.GetPosition());
287       aSecondaryTrack->SetTouchableHandle(stepData.GetPostStepPoint()->GetTouchableHandle());
288       aSecondaryTrack->SetParentID(trackData.GetTrackID());
289       aSecondaryTrack->SetCreatorModelID(secID);
290       aParticleChange.AddSecondary(aSecondaryTrack);
291 
292     }
293   }
294   return G4VDiscreteProcess::PostStepDoIt(trackData, stepData);
295 }
296 
297 ///////////////////////////////////////////////////////////////////////////////
298 G4double G4SynchrotronRadiation::InvSynFracInt(G4double x)
299 // direct generation
300 {
301   // from 0 to 0.7
302   static constexpr G4double aa1           = 0;
303   static constexpr G4double aa2           = 0.7;
304   static constexpr G4int ncheb1           = 27;
305   static constexpr G4double cheb1[ncheb1] = {
306     1.22371665676046468821,     0.108956475422163837267,
307     0.0383328524358594396134,   0.00759138369340257753721,
308     0.00205712048644963340914,  0.000497810783280019308661,
309     0.000130743691810302187818, 0.0000338168760220395409734,
310     8.97049680900520817728e-6,  2.38685472794452241466e-6,
311     6.41923109149104165049e-7,  1.73549898982749277843e-7,
312     4.72145949240790029153e-8,  1.29039866111999149636e-8,
313     3.5422080787089834182e-9,   9.7594757336403784905e-10,
314     2.6979510184976065731e-10,  7.480422622550977077e-11,
315     2.079598176402699913e-11,   5.79533622220841193e-12,
316     1.61856011449276096e-12,    4.529450993473807e-13,
317     1.2698603951096606e-13,     3.566117394511206e-14,
318     1.00301587494091e-14,       2.82515346447219e-15,
319     7.9680747949792e-16
320   };
321   //   from 0.7 to 0.9132260271183847
322   static constexpr G4double aa3           = 0.9132260271183847;
323   static constexpr G4int ncheb2           = 27;
324   static constexpr G4double cheb2[ncheb2] = {
325     1.1139496701107756,     0.3523967429328067,     0.0713849171926623,
326     0.01475818043595387,    0.003381255637322462,   0.0008228057599452224,
327     0.00020785506681254216, 0.00005390169253706556, 0.000014250571923902464,
328     3.823880733161044e-6,   1.0381966089136036e-6,  2.8457557457837253e-7,
329     7.86223332179956e-8,    2.1866609342508474e-8,  6.116186259857143e-9,
330     1.7191233618437565e-9,  4.852755117740807e-10,  1.3749966961763457e-10,
331     3.908961987062447e-11,  1.1146253766895824e-11, 3.1868887323415814e-12,
332     9.134319791300977e-13,  2.6211077371181566e-13, 7.588643377757906e-14,
333     2.1528376972619e-14,    6.030906040404772e-15,  1.9549163926819867e-15
334   };
335   // Chebyshev with exp/log  scale
336   // a = -Log[1 - SynFracInt[1]]; b = -Log[1 - SynFracInt[7]];
337   static constexpr G4double aa4           = 2.4444485538746025480;
338   static constexpr G4double aa5           = 9.3830728608909477079;
339   static constexpr G4int ncheb3           = 28;
340   static constexpr G4double cheb3[ncheb3] = {
341     1.2292683840435586977,        0.160353449247864455879,
342     -0.0353559911947559448721,    0.00776901561223573936985,
343     -0.00165886451971685133259,   0.000335719118906954279467,
344     -0.0000617184951079161143187, 9.23534039743246708256e-6,
345     -6.06747198795168022842e-7,   -3.07934045961999778094e-7,
346     1.98818772614682367781e-7,    -8.13909971567720135413e-8,
347     2.84298174969641838618e-8,    -9.12829766621316063548e-9,
348     2.77713868004820551077e-9,    -8.13032767247834023165e-10,
349     2.31128525568385247392e-10,   -6.41796873254200220876e-11,
350     1.74815310473323361543e-11,   -4.68653536933392363045e-12,
351     1.24016595805520752748e-12,   -3.24839432979935522159e-13,
352     8.44601465226513952994e-14,   -2.18647276044246803998e-14,
353     5.65407548745690689978e-15,   -1.46553625917463067508e-15,
354     3.82059606377570462276e-16,   -1.00457896653436912508e-16
355   };
356   static constexpr G4double aa6           = 33.122936966163038145;
357   static constexpr G4int ncheb4           = 27;
358   static constexpr G4double cheb4[ncheb4] = {
359     1.69342658227676741765,      0.0742766400841232319225,
360     -0.019337880608635717358,    0.00516065527473364110491,
361     -0.00139342012990307729473,  0.000378549864052022522193,
362     -0.000103167085583785340215, 0.0000281543441271412178337,
363     -7.68409742018258198651e-6,  2.09543221890204537392e-6,
364     -5.70493140367526282946e-7,  1.54961164548564906446e-7,
365     -4.19665599629607704794e-8,  1.13239680054166507038e-8,
366     -3.04223563379021441863e-9,  8.13073745977562957997e-10,
367     -2.15969415476814981374e-10, 5.69472105972525594811e-11,
368     -1.48844799572430829499e-11, 3.84901514438304484973e-12,
369     -9.82222575944247161834e-13, 2.46468329208292208183e-13,
370     -6.04953826265982691612e-14, 1.44055805710671611984e-14,
371     -3.28200813577388740722e-15, 6.96566359173765367675e-16,
372     -1.294122794852896275e-16
373   };
374 
375   if(x < aa2)
376     return x * x * x * Chebyshev(aa1, aa2, cheb1, ncheb1, x);
377   else if(x < aa3)
378     return Chebyshev(aa2, aa3, cheb2, ncheb2, x);
379   else if(x < 1 - 0.0000841363)
380   {
381     G4double y = -G4Log(1 - x);
382     return y * Chebyshev(aa4, aa5, cheb3, ncheb3, y);
383   }
384   else
385   {
386     G4double y = -G4Log(1 - x);
387     return y * Chebyshev(aa5, aa6, cheb4, ncheb4, y);
388   }
389 }
390 
391 G4double G4SynchrotronRadiation::GetRandomEnergySR(G4double gamma,
392                                                    G4double perpB,
393                                                    G4double mass_c2)
394 {
395   static const G4double fEnergyConst =
396     1.5 * c_light * c_light * eplus * hbar_Planck;
397   G4double Ecr = fEnergyConst * gamma * gamma * perpB / mass_c2;
398 
399   if(verboseLevel > 0 && FirstTime1)
400   {
401     // mean and rms of photon energy
402     G4double Emean = 8. / (15. * std::sqrt(3.)) * Ecr;
403     G4double E_rms = std::sqrt(211. / 675.) * Ecr;
404     G4long prec     = G4cout.precision();
405     G4cout << "G4SynchrotronRadiation::GetRandomEnergySR :" << '\n'
406            << std::setprecision(4) << "  Ecr   = " << G4BestUnit(Ecr, "Energy")
407            << '\n'
408            << "  Emean = " << G4BestUnit(Emean, "Energy") << '\n'
409            << "  E_rms = " << G4BestUnit(E_rms, "Energy") << G4endl;
410     FirstTime1 = false;
411     G4cout.precision(prec);
412   }
413 
414   G4double energySR = Ecr * InvSynFracInt(G4UniformRand());
415   return energySR;
416 }
417 
418 ///////////////////////////////////////////////////////////////////////////////
419 void G4SynchrotronRadiation::BuildPhysicsTable(const G4ParticleDefinition& part)
420 {
421   if(0 < verboseLevel && &part == G4Electron::Electron())
422     ProcessDescription(G4cout);
423   // same for all particles, print only for one (electron)
424 }
425 
426 ///////////////////////////////////////////////////////////////////////////////
427 void G4SynchrotronRadiation::ProcessDescription(std::ostream& out) const
428 {
429   out << GetProcessName()
430       << ":  Incoherent Synchrotron Radiation\n"
431          "Good description for long magnets at all energies.\n";
432 }
433