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 9.6)


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