ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SynchrotronRadiationInMat.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SynchrotronRadiationInMat.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 // 19-05-06, V.Ivanchenko rename from G4SynchrotronRadiation
37 //
38 //
40 
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Integrator.hh"
45 #include "G4EmProcessSubType.hh"
46 
48 //
49 // Constant for calculation of mean free path
50 //
51 
52 const G4double
55 
57 //
58 // Constant for calculation of characterictic energy
59 //
60 
61 const G4double
64 
66 //
67 // Array of integral probability of synchrotron photons:
68 //
69 // the corresponding energy = 0.0001*i*i*(characteristic energy)
70 //
71 
72 const G4double
74 {
75  1.000000e+00, 9.428859e-01, 9.094095e-01, 8.813971e-01, 8.565154e-01,
76  8.337008e-01, 8.124961e-01, 7.925217e-01, 7.735517e-01, 7.554561e-01,
77  7.381233e-01, 7.214521e-01, 7.053634e-01, 6.898006e-01, 6.747219e-01,
78  6.600922e-01, 6.458793e-01, 6.320533e-01, 6.185872e-01, 6.054579e-01,
79  5.926459e-01, 5.801347e-01, 5.679103e-01, 5.559604e-01, 5.442736e-01,
80  5.328395e-01, 5.216482e-01, 5.106904e-01, 4.999575e-01, 4.894415e-01,
81  4.791351e-01, 4.690316e-01, 4.591249e-01, 4.494094e-01, 4.398800e-01,
82  4.305320e-01, 4.213608e-01, 4.123623e-01, 4.035325e-01, 3.948676e-01,
83  3.863639e-01, 3.780179e-01, 3.698262e-01, 3.617858e-01, 3.538933e-01,
84  3.461460e-01, 3.385411e-01, 3.310757e-01, 3.237474e-01, 3.165536e-01,
85  3.094921e-01, 3.025605e-01, 2.957566e-01, 2.890784e-01, 2.825237e-01,
86  2.760907e-01, 2.697773e-01, 2.635817e-01, 2.575020e-01, 2.515365e-01,
87  2.456834e-01, 2.399409e-01, 2.343074e-01, 2.287812e-01, 2.233607e-01,
88  2.180442e-01, 2.128303e-01, 2.077174e-01, 2.027040e-01, 1.977885e-01,
89  1.929696e-01, 1.882457e-01, 1.836155e-01, 1.790775e-01, 1.746305e-01,
90  1.702730e-01, 1.660036e-01, 1.618212e-01, 1.577243e-01, 1.537117e-01,
91  1.497822e-01, 1.459344e-01, 1.421671e-01, 1.384791e-01, 1.348691e-01,
92  1.313360e-01, 1.278785e-01, 1.244956e-01, 1.211859e-01, 1.179483e-01,
93  1.147818e-01, 1.116850e-01, 1.086570e-01, 1.056966e-01, 1.028026e-01,
94  9.997405e-02, 9.720975e-02, 9.450865e-02, 9.186969e-02, 8.929179e-02,
95  8.677391e-02, 8.431501e-02, 8.191406e-02, 7.957003e-02, 7.728192e-02,
96  7.504872e-02, 7.286944e-02, 7.074311e-02, 6.866874e-02, 6.664538e-02,
97  6.467208e-02, 6.274790e-02, 6.087191e-02, 5.904317e-02, 5.726079e-02,
98  5.552387e-02, 5.383150e-02, 5.218282e-02, 5.057695e-02, 4.901302e-02,
99  4.749020e-02, 4.600763e-02, 4.456450e-02, 4.315997e-02, 4.179325e-02,
100  4.046353e-02, 3.917002e-02, 3.791195e-02, 3.668855e-02, 3.549906e-02,
101  3.434274e-02, 3.321884e-02, 3.212665e-02, 3.106544e-02, 3.003452e-02,
102  2.903319e-02, 2.806076e-02, 2.711656e-02, 2.619993e-02, 2.531021e-02,
103  2.444677e-02, 2.360897e-02, 2.279620e-02, 2.200783e-02, 2.124327e-02,
104  2.050194e-02, 1.978324e-02, 1.908662e-02, 1.841151e-02, 1.775735e-02,
105  1.712363e-02, 1.650979e-02, 1.591533e-02, 1.533973e-02, 1.478250e-02,
106  1.424314e-02, 1.372117e-02, 1.321613e-02, 1.272755e-02, 1.225498e-02,
107  1.179798e-02, 1.135611e-02, 1.092896e-02, 1.051609e-02, 1.011712e-02,
108  9.731635e-03, 9.359254e-03, 8.999595e-03, 8.652287e-03, 8.316967e-03,
109  7.993280e-03, 7.680879e-03, 7.379426e-03, 7.088591e-03, 6.808051e-03,
110  6.537491e-03, 6.276605e-03, 6.025092e-03, 5.782661e-03, 5.549027e-03,
111  5.323912e-03, 5.107045e-03, 4.898164e-03, 4.697011e-03, 4.503336e-03,
112  4.316896e-03, 4.137454e-03, 3.964780e-03, 3.798649e-03, 3.638843e-03,
113  3.485150e-03, 3.337364e-03, 3.195284e-03, 3.058715e-03, 2.927469e-03,
114  2.801361e-03, 2.680213e-03, 2.563852e-03, 2.452110e-03, 2.344824e-03
115 };
116 
118 //
119 // Constructor
120 //
121 
123  G4ProcessType type):G4VDiscreteProcess (processName, type),
124  LowestKineticEnergy (10.*keV),
125  theGamma (G4Gamma::Gamma() ),
126  theElectron ( G4Electron::Electron() ),
127  thePositron ( G4Positron::Positron() ),
128  fAlpha(0.0), fRootNumber(80),
129  fVerboseLevel( verboseLevel )
130 {
132 
133  fFieldPropagator = transportMgr->GetPropagatorInField();
137  fPsiGamma = fEta = fOrderAngleK = 0.0;
138 }
139 
141 //
142 // Destructor
143 //
144 
146 {}
147 
148 
149 G4bool
151 {
152  return ( ( &particle == (const G4ParticleDefinition *)theElectron ) ||
153  ( &particle == (const G4ParticleDefinition *)thePositron ));
154 }
155 
157 {
158  return fLambdaConst;
159 }
160 
162 {
163  return fEnergyConst;
164 }
165 
167 //
168 //
169 // Production of synchrotron X-ray photon
170 // GEANT4 internal units.
171 //
172 
173 
174 G4double
176  G4double,
178 {
179  // gives the MeanFreePath in GEANT4 internal units
180  G4double MeanFreePath;
181 
182  const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
183  // G4Material* aMaterial = trackData.GetMaterial();
184 
185  //G4bool isOutRange ;
186 
187  *condition = NotForced ;
188 
189  G4double gamma = aDynamicParticle->GetTotalEnergy()/
190  aDynamicParticle->GetMass();
191 
192  G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
193 
194  G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
195 
196  if ( KineticEnergy < LowestKineticEnergy || gamma < 1.0e3 ) MeanFreePath = DBL_MAX;
197  else
198  {
199 
200  G4ThreeVector FieldValue;
201  const G4Field* pField = nullptr;
202 
203  G4FieldManager* fieldMgr=nullptr;
204  G4bool fieldExertsForce = false;
205 
206  if( (particleCharge != 0.0) )
207  {
208  fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
209 
210  if ( fieldMgr != nullptr )
211  {
212  // If the field manager has no field, there is no field !
213 
214  fieldExertsForce = ( fieldMgr->GetDetectorField() != nullptr );
215  }
216  }
217  if ( fieldExertsForce )
218  {
219  pField = fieldMgr->GetDetectorField() ;
220  G4ThreeVector globPosition = trackData.GetPosition();
221 
222  G4double globPosVec[4], FieldValueVec[6];
223 
224  globPosVec[0] = globPosition.x();
225  globPosVec[1] = globPosition.y();
226  globPosVec[2] = globPosition.z();
227  globPosVec[3] = trackData.GetGlobalTime();
228 
229  pField->GetFieldValue( globPosVec, FieldValueVec );
230 
231  FieldValue = G4ThreeVector( FieldValueVec[0],
232  FieldValueVec[1],
233  FieldValueVec[2] );
234 
235 
236 
237  G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
238  G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum) ;
239  G4double perpB = unitMcrossB.mag() ;
240  G4double beta = aDynamicParticle->GetTotalMomentum()/
241  (aDynamicParticle->GetTotalEnergy() );
242 
243  if( perpB > 0.0 ) MeanFreePath = fLambdaConst*beta/perpB;
244  else MeanFreePath = DBL_MAX;
245  }
246  else MeanFreePath = DBL_MAX;
247  }
248  if(fVerboseLevel > 0)
249  {
250  G4cout<<"G4SynchrotronRadiationInMat::MeanFreePath = "<<MeanFreePath/m<<" m"<<G4endl;
251  }
252  return MeanFreePath;
253 }
254 
256 //
257 //
258 
261  const G4Step& stepData )
262 
263 {
264  aParticleChange.Initialize(trackData);
265 
266  const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
267 
268  G4double gamma = aDynamicParticle->GetTotalEnergy()/
269  (aDynamicParticle->GetMass() );
270 
271  if(gamma <= 1.0e3 )
272  {
273  return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
274  }
275  G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
276 
277  G4ThreeVector FieldValue;
278  const G4Field* pField = nullptr ;
279 
280  G4FieldManager* fieldMgr=nullptr;
281  G4bool fieldExertsForce = false;
282 
283  if( (particleCharge != 0.0) )
284  {
285  fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
286  if ( fieldMgr != nullptr )
287  {
288  // If the field manager has no field, there is no field !
289 
290  fieldExertsForce = ( fieldMgr->GetDetectorField() != nullptr );
291  }
292  }
293  if ( fieldExertsForce )
294  {
295  pField = fieldMgr->GetDetectorField() ;
296  G4ThreeVector globPosition = trackData.GetPosition() ;
297  G4double globPosVec[4], FieldValueVec[6] ;
298  globPosVec[0] = globPosition.x() ;
299  globPosVec[1] = globPosition.y() ;
300  globPosVec[2] = globPosition.z() ;
301  globPosVec[3] = trackData.GetGlobalTime();
302 
303  pField->GetFieldValue( globPosVec, FieldValueVec ) ;
304  FieldValue = G4ThreeVector( FieldValueVec[0],
305  FieldValueVec[1],
306  FieldValueVec[2] );
307 
308  G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
309  G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
310  G4double perpB = unitMcrossB.mag() ;
311  if(perpB > 0.0)
312  {
313  // M-C of synchrotron photon energy
314 
315  G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
316 
317  if(fVerboseLevel > 0)
318  {
319  G4cout<<"SR photon energy = "<<energyOfSR/keV<<" keV"<<G4endl;
320  }
321  // check against insufficient energy
322 
323  if( energyOfSR <= 0.0 )
324  {
325  return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
326  }
327  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
329  particleDirection = aDynamicParticle->GetMomentumDirection();
330 
331  // M-C of its direction, simplified dipole busted approach
332 
333  // G4double Teta = G4UniformRand()/gamma ; // Very roughly
334 
335  G4double cosTheta, sinTheta, fcos, beta;
336 
337  do
338  {
339  cosTheta = 1. - 2.*G4UniformRand();
340  fcos = (1 + cosTheta*cosTheta)*0.5;
341  }
342  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
343  while( fcos < G4UniformRand() );
344 
345  beta = std::sqrt(1. - 1./(gamma*gamma));
346 
347  cosTheta = (cosTheta + beta)/(1. + beta*cosTheta);
348 
349  if( cosTheta > 1. ) cosTheta = 1.;
350  if( cosTheta < -1. ) cosTheta = -1.;
351 
352  sinTheta = std::sqrt(1. - cosTheta*cosTheta );
353 
354  G4double Phi = twopi * G4UniformRand() ;
355 
356  G4double dirx = sinTheta*std::cos(Phi) ,
357  diry = sinTheta*std::sin(Phi) ,
358  dirz = cosTheta;
359 
360  G4ThreeVector gammaDirection ( dirx, diry, dirz);
361  gammaDirection.rotateUz(particleDirection);
362 
363  // polarization of new gamma
364 
365  // G4double sx = std::cos(Teta)*std::cos(Phi);
366  // G4double sy = std::cos(Teta)*std::sin(Phi);
367  // G4double sz = -std::sin(Teta);
368 
369  G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
370  gammaPolarization = gammaPolarization.unit();
371 
372  // (sx, sy, sz);
373  // gammaPolarization.rotateUz(particleDirection);
374 
375  // create G4DynamicParticle object for the SR photon
376 
378  gammaDirection,
379  energyOfSR );
380  aGamma->SetPolarization( gammaPolarization.x(),
381  gammaPolarization.y(),
382  gammaPolarization.z() );
383 
384 
386  aParticleChange.AddSecondary(aGamma);
387 
388  // Update the incident particle
389 
390  G4double newKinEnergy = kineticEnergy - energyOfSR ;
391 
392  if (newKinEnergy > 0.)
393  {
394  aParticleChange.ProposeMomentumDirection( particleDirection );
395  aParticleChange.ProposeEnergy( newKinEnergy );
397  }
398  else
399  {
402  G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();
403  if (charge<0.)
404  {
406  }
407  else
408  {
410  }
411  }
412  }
413  else
414  {
415  return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
416  }
417  }
418  return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
419 }
420 
421 
422 G4double
424  const G4Step& )
425 
426 {
427  G4int i ;
428  G4double energyOfSR = -1.0 ;
429  //G4Material* aMaterial=trackData.GetMaterial() ;
430 
431  const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
432 
433  G4double gamma = aDynamicParticle->GetTotalEnergy()/
434  (aDynamicParticle->GetMass() ) ;
435 
436  G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
437 
438  G4ThreeVector FieldValue;
439  const G4Field* pField = nullptr ;
440 
441  G4FieldManager* fieldMgr=nullptr;
442  G4bool fieldExertsForce = false;
443 
444  if( (particleCharge != 0.0) )
445  {
446  fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
447  if ( fieldMgr != nullptr )
448  {
449  // If the field manager has no field, there is no field !
450 
451  fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
452  }
453  }
454  if ( fieldExertsForce )
455  {
456  pField = fieldMgr->GetDetectorField();
457  G4ThreeVector globPosition = trackData.GetPosition();
458  G4double globPosVec[3], FieldValueVec[3];
459 
460  globPosVec[0] = globPosition.x();
461  globPosVec[1] = globPosition.y();
462  globPosVec[2] = globPosition.z();
463 
464  pField->GetFieldValue( globPosVec, FieldValueVec );
465  FieldValue = G4ThreeVector( FieldValueVec[0],
466  FieldValueVec[1],
467  FieldValueVec[2] );
468 
469  G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection();
470  G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum) ;
471  G4double perpB = unitMcrossB.mag();
472  if( perpB > 0.0 )
473  {
474  // M-C of synchrotron photon energy
475 
476  G4double random = G4UniformRand() ;
477  for(i=0;i<200;i++)
478  {
479  if(random >= fIntegralProbabilityOfSR[i]) break ;
480  }
481  energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
482 
483  // check against insufficient energy
484 
485  if(energyOfSR <= 0.0)
486  {
487  return -1.0 ;
488  }
489  //G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
490  //G4ParticleMomentum
491  //particleDirection = aDynamicParticle->GetMomentumDirection();
492 
493  // Gamma production cut in this material
494  //G4double
495  //gammaEnergyCut = (G4Gamma::GetCutsInEnergy())[aMaterial->GetIndex()];
496 
497  // SR photon has energy more than the current material cut
498  // M-C of its direction
499 
500  //G4double Teta = G4UniformRand()/gamma ; // Very roughly
501 
502  //G4double Phi = twopi * G4UniformRand() ;
503  }
504  else
505  {
506  return -1.0 ;
507  }
508  }
509  return energyOfSR ;
510 }
511 
513 //
514 //
515 
517 {
518  G4int i, iMax;
519  G4double energySR, random, position;
520 
521  iMax = 200;
522  random = G4UniformRand();
523 
524  for( i = 0; i < iMax; i++ )
525  {
526  if( random >= fIntegralProbabilityOfSR[i] ) break;
527  }
528  if(i <= 0 ) position = G4UniformRand(); // 0.
529  else if( i>= iMax) position = G4double(iMax);
530  else position = i + G4UniformRand(); // -1
531  //
532  // it was in initial implementation:
533  // energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
534 
535  energySR = 0.0001*position*position*fEnergyConst*gamma*gamma*perpB;
536 
537  if( energySR < 0. ) energySR = 0.;
538 
539  return energySR;
540 }
541 
543 //
544 // return
545 
547 {
548  G4double result, hypCos2, hypCos=std::cosh(t);
549 
550  hypCos2 = hypCos*hypCos;
551  result = std::cosh(5.*t/3.)*std::exp(t-fKsi*hypCos); // fKsi > 0. !
552  result /= hypCos2;
553  return result;
554 }
555 
557 //
558 // return the probability to emit SR photon with relative energy
559 // energy/energy_c >= ksi
560 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
561 
563 {
564  if (ksi <= 0.) return 1.0;
565  fKsi = ksi; // should be > 0. !
566  G4int n;
567  G4double result, a;
568 
569  a = fAlpha; // always = 0.
570  n = fRootNumber; // around default = 80
571 
573 
574  result = integral.Laguerre(this,
576 
577  result *= 3./5./pi;
578 
579  return result;
580 }
581 
583 //
584 // return an auxiliary function for K_5/3 integral representation
585 
587 {
588  G4double result, hypCos=std::cosh(t);
589 
590  result = std::cosh(5.*t/3.)*std::exp(t - fKsi*hypCos); // fKsi > 0. !
591  result /= hypCos;
592  return result;
593 }
594 
596 //
597 // return the probability to emit SR photon energy with relative energy
598 // energy/energy_c >= ksi
599 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
600 
602 {
603  if (ksi <= 0.) return 1.0;
604  fKsi = ksi; // should be > 0. !
605  G4int n;
606  G4double result, a;
607 
608  a = fAlpha; // always = 0.
609  n = fRootNumber; // around default = 80
610 
612 
613  result = integral.Laguerre(this,
615 
616  result *= 9.*std::sqrt(3.)*ksi/8./pi;
617 
618  return result;
619 }
620 
622 //
623 //
624 
626 {
627  G4double result, hypCos=std::cosh(t);
628 
629  result = std::cosh(fOrderAngleK*t)*std::exp(t - fEta*hypCos); // fEta > 0. !
630  result /= hypCos;
631  return result;
632 }
633 
635 //
636 // Return K 1/3 or 2/3 for angular distribution
637 
639 {
640  fEta = eta; // should be > 0. !
641  G4int n;
642  G4double result, a;
643 
644  a = fAlpha; // always = 0.
645  n = fRootNumber; // around default = 80
646 
648 
649  result = integral.Laguerre(this,
651 
652  return result;
653 }
654 
656 //
657 // Relative angle diff distribution for given fKsi, which is set externally
658 
660 {
661  G4double result, funK, funK2, gpsi2 = gpsi*gpsi;
662 
663  fPsiGamma = gpsi;
664  fEta = 0.5*fKsi*(1. + gpsi2)*std::sqrt(1. + gpsi2);
665 
666  fOrderAngleK = 1./3.;
667  funK = GetAngleK(fEta);
668  funK2 = funK*funK;
669 
670  result = gpsi2*funK2/(1. + gpsi2);
671 
672  fOrderAngleK = 2./3.;
673  funK = GetAngleK(fEta);
674  funK2 = funK*funK;
675 
676  result += funK2;
677  result *= (1. + gpsi2)*fKsi;
678 
679  return result;
680 }
681 
682 
684