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 ]

Diff markup

Differences between /processes/electromagnetic/xrays/src/G4SynchrotronRadiation.cc (Version 11.3.0) and /processes/electromagnetic/xrays/src/G4SynchrotronRadiation.cc (Version 10.0.p1)


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