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