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