Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/xrays/src/G4SynchrotronRadiationInMat.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/G4SynchrotronRadiationInMat.cc (Version 11.3.0) and /processes/electromagnetic/xrays/src/G4SynchrotronRadiationInMat.cc (Version 9.0.p2)


  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: G4SynchrotronRadiationInMat.cc,v 1.2 2006/06/29 19:56:17 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-08-03-patch-01 $
                                                   >>  29 //
 26 // -------------------------------------------     30 // --------------------------------------------------------------
 27 //      GEANT 4 class implementation file          31 //      GEANT 4 class implementation file
                                                   >>  32 //      CERN Geneva Switzerland
 28 //                                                 33 //
 29 //      History: first implementation,             34 //      History: first implementation,
 30 //      21-5-98 V.Grichine                         35 //      21-5-98 V.Grichine
 31 //      28-05-01, V.Ivanchenko minor changes t     36 //      28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation
 32 //      04.03.05, V.Grichine: get local field      37 //      04.03.05, V.Grichine: get local field interface
 33 //      19-05-06, V.Ivanchenko rename from G4S     38 //      19-05-06, V.Ivanchenko rename from G4SynchrotronRadiation
 34 //                                                 39 //
                                                   >>  40 //
 35 //////////////////////////////////////////////     41 ///////////////////////////////////////////////////////////////////////////
 36                                                    42 
 37 #include "G4SynchrotronRadiationInMat.hh"          43 #include "G4SynchrotronRadiationInMat.hh"
 38                                                << 
 39 #include "G4EmProcessSubType.hh"               << 
 40 #include "G4Field.hh"                          << 
 41 #include "G4FieldManager.hh"                   << 
 42 #include "G4Integrator.hh"                         44 #include "G4Integrator.hh"
 43 #include "G4PhysicalConstants.hh"              << 
 44 #include "G4PropagatorInField.hh"              << 
 45 #include "G4SystemOfUnits.hh"                  << 
 46 #include "G4PhysicsModelCatalog.hh"            << 
 47                                                << 
 48 const G4double G4SynchrotronRadiationInMat::fI << 
 49   1.000000e+00, 9.428859e-01, 9.094095e-01, 8. << 
 50   8.337008e-01, 8.124961e-01, 7.925217e-01, 7. << 
 51   7.381233e-01, 7.214521e-01, 7.053634e-01, 6. << 
 52   6.600922e-01, 6.458793e-01, 6.320533e-01, 6. << 
 53   5.926459e-01, 5.801347e-01, 5.679103e-01, 5. << 
 54   5.328395e-01, 5.216482e-01, 5.106904e-01, 4. << 
 55   4.791351e-01, 4.690316e-01, 4.591249e-01, 4. << 
 56   4.305320e-01, 4.213608e-01, 4.123623e-01, 4. << 
 57   3.863639e-01, 3.780179e-01, 3.698262e-01, 3. << 
 58   3.461460e-01, 3.385411e-01, 3.310757e-01, 3. << 
 59   3.094921e-01, 3.025605e-01, 2.957566e-01, 2. << 
 60   2.760907e-01, 2.697773e-01, 2.635817e-01, 2. << 
 61   2.456834e-01, 2.399409e-01, 2.343074e-01, 2. << 
 62   2.180442e-01, 2.128303e-01, 2.077174e-01, 2. << 
 63   1.929696e-01, 1.882457e-01, 1.836155e-01, 1. << 
 64   1.702730e-01, 1.660036e-01, 1.618212e-01, 1. << 
 65   1.497822e-01, 1.459344e-01, 1.421671e-01, 1. << 
 66   1.313360e-01, 1.278785e-01, 1.244956e-01, 1. << 
 67   1.147818e-01, 1.116850e-01, 1.086570e-01, 1. << 
 68   9.997405e-02, 9.720975e-02, 9.450865e-02, 9. << 
 69   8.677391e-02, 8.431501e-02, 8.191406e-02, 7. << 
 70   7.504872e-02, 7.286944e-02, 7.074311e-02, 6. << 
 71   6.467208e-02, 6.274790e-02, 6.087191e-02, 5. << 
 72   5.552387e-02, 5.383150e-02, 5.218282e-02, 5. << 
 73   4.749020e-02, 4.600763e-02, 4.456450e-02, 4. << 
 74   4.046353e-02, 3.917002e-02, 3.791195e-02, 3. << 
 75   3.434274e-02, 3.321884e-02, 3.212665e-02, 3. << 
 76   2.903319e-02, 2.806076e-02, 2.711656e-02, 2. << 
 77   2.444677e-02, 2.360897e-02, 2.279620e-02, 2. << 
 78   2.050194e-02, 1.978324e-02, 1.908662e-02, 1. << 
 79   1.712363e-02, 1.650979e-02, 1.591533e-02, 1. << 
 80   1.424314e-02, 1.372117e-02, 1.321613e-02, 1. << 
 81   1.179798e-02, 1.135611e-02, 1.092896e-02, 1. << 
 82   9.731635e-03, 9.359254e-03, 8.999595e-03, 8. << 
 83   7.993280e-03, 7.680879e-03, 7.379426e-03, 7. << 
 84   6.537491e-03, 6.276605e-03, 6.025092e-03, 5. << 
 85   5.323912e-03, 5.107045e-03, 4.898164e-03, 4. << 
 86   4.316896e-03, 4.137454e-03, 3.964780e-03, 3. << 
 87   3.485150e-03, 3.337364e-03, 3.195284e-03, 3. << 
 88   2.801361e-03, 2.680213e-03, 2.563852e-03, 2. << 
 89 };                                             << 
 90                                                    45 
 91 ////////////////////////////////////////////// <<  46 using namespace std;
                                                   >>  47 
                                                   >>  48 ////////////////////////////////////////////////////////////////////
                                                   >>  49 //
                                                   >>  50 // Constant for calculation of mean free path
                                                   >>  51 //
                                                   >>  52 
                                                   >>  53 const G4double
                                                   >>  54 G4SynchrotronRadiationInMat::fLambdaConst = sqrt(3.0)*electron_mass_c2/
                                                   >>  55                                        (2.5*fine_structure_const*eplus*c_light) ;
                                                   >>  56 
                                                   >>  57 /////////////////////////////////////////////////////////////////////
                                                   >>  58 //
                                                   >>  59 // Constant for calculation of characterictic energy
                                                   >>  60 //
                                                   >>  61 
                                                   >>  62 const G4double
                                                   >>  63 G4SynchrotronRadiationInMat::fEnergyConst = 1.5*c_light*c_light*eplus*hbar_Planck/
                                                   >>  64                                        electron_mass_c2  ;
                                                   >>  65 
                                                   >>  66 ////////////////////////////////////////////////////////////////////
                                                   >>  67 //
                                                   >>  68 // Array of integral probability of synchrotron photons:
                                                   >>  69 //
                                                   >>  70 // the corresponding energy = 0.0001*i*i*(characteristic energy)
                                                   >>  71 //
                                                   >>  72 
                                                   >>  73 const G4double
                                                   >>  74 G4SynchrotronRadiationInMat::fIntegralProbabilityOfSR[200] =
                                                   >>  75 {
                                                   >>  76   1.000000e+00, 9.428859e-01, 9.094095e-01, 8.813971e-01, 8.565154e-01,
                                                   >>  77   8.337008e-01, 8.124961e-01, 7.925217e-01, 7.735517e-01, 7.554561e-01,
                                                   >>  78   7.381233e-01, 7.214521e-01, 7.053634e-01, 6.898006e-01, 6.747219e-01,
                                                   >>  79   6.600922e-01, 6.458793e-01, 6.320533e-01, 6.185872e-01, 6.054579e-01,
                                                   >>  80   5.926459e-01, 5.801347e-01, 5.679103e-01, 5.559604e-01, 5.442736e-01,
                                                   >>  81   5.328395e-01, 5.216482e-01, 5.106904e-01, 4.999575e-01, 4.894415e-01,
                                                   >>  82   4.791351e-01, 4.690316e-01, 4.591249e-01, 4.494094e-01, 4.398800e-01,
                                                   >>  83   4.305320e-01, 4.213608e-01, 4.123623e-01, 4.035325e-01, 3.948676e-01,
                                                   >>  84   3.863639e-01, 3.780179e-01, 3.698262e-01, 3.617858e-01, 3.538933e-01,
                                                   >>  85   3.461460e-01, 3.385411e-01, 3.310757e-01, 3.237474e-01, 3.165536e-01,
                                                   >>  86   3.094921e-01, 3.025605e-01, 2.957566e-01, 2.890784e-01, 2.825237e-01,
                                                   >>  87   2.760907e-01, 2.697773e-01, 2.635817e-01, 2.575020e-01, 2.515365e-01,
                                                   >>  88   2.456834e-01, 2.399409e-01, 2.343074e-01, 2.287812e-01, 2.233607e-01,
                                                   >>  89   2.180442e-01, 2.128303e-01, 2.077174e-01, 2.027040e-01, 1.977885e-01,
                                                   >>  90   1.929696e-01, 1.882457e-01, 1.836155e-01, 1.790775e-01, 1.746305e-01,
                                                   >>  91   1.702730e-01, 1.660036e-01, 1.618212e-01, 1.577243e-01, 1.537117e-01,
                                                   >>  92   1.497822e-01, 1.459344e-01, 1.421671e-01, 1.384791e-01, 1.348691e-01,
                                                   >>  93   1.313360e-01, 1.278785e-01, 1.244956e-01, 1.211859e-01, 1.179483e-01,
                                                   >>  94   1.147818e-01, 1.116850e-01, 1.086570e-01, 1.056966e-01, 1.028026e-01,
                                                   >>  95   9.997405e-02, 9.720975e-02, 9.450865e-02, 9.186969e-02, 8.929179e-02,
                                                   >>  96   8.677391e-02, 8.431501e-02, 8.191406e-02, 7.957003e-02, 7.728192e-02,
                                                   >>  97   7.504872e-02, 7.286944e-02, 7.074311e-02, 6.866874e-02, 6.664538e-02,
                                                   >>  98   6.467208e-02, 6.274790e-02, 6.087191e-02, 5.904317e-02, 5.726079e-02,
                                                   >>  99   5.552387e-02, 5.383150e-02, 5.218282e-02, 5.057695e-02, 4.901302e-02,
                                                   >> 100   4.749020e-02, 4.600763e-02, 4.456450e-02, 4.315997e-02, 4.179325e-02,
                                                   >> 101   4.046353e-02, 3.917002e-02, 3.791195e-02, 3.668855e-02, 3.549906e-02,
                                                   >> 102   3.434274e-02, 3.321884e-02, 3.212665e-02, 3.106544e-02, 3.003452e-02,
                                                   >> 103   2.903319e-02, 2.806076e-02, 2.711656e-02, 2.619993e-02, 2.531021e-02,
                                                   >> 104   2.444677e-02, 2.360897e-02, 2.279620e-02, 2.200783e-02, 2.124327e-02,
                                                   >> 105   2.050194e-02, 1.978324e-02, 1.908662e-02, 1.841151e-02, 1.775735e-02,
                                                   >> 106   1.712363e-02, 1.650979e-02, 1.591533e-02, 1.533973e-02, 1.478250e-02,
                                                   >> 107   1.424314e-02, 1.372117e-02, 1.321613e-02, 1.272755e-02, 1.225498e-02,
                                                   >> 108   1.179798e-02, 1.135611e-02, 1.092896e-02, 1.051609e-02, 1.011712e-02,
                                                   >> 109   9.731635e-03, 9.359254e-03, 8.999595e-03, 8.652287e-03, 8.316967e-03,
                                                   >> 110   7.993280e-03, 7.680879e-03, 7.379426e-03, 7.088591e-03, 6.808051e-03,
                                                   >> 111   6.537491e-03, 6.276605e-03, 6.025092e-03, 5.782661e-03, 5.549027e-03,
                                                   >> 112   5.323912e-03, 5.107045e-03, 4.898164e-03, 4.697011e-03, 4.503336e-03,
                                                   >> 113   4.316896e-03, 4.137454e-03, 3.964780e-03, 3.798649e-03, 3.638843e-03,
                                                   >> 114   3.485150e-03, 3.337364e-03, 3.195284e-03, 3.058715e-03, 2.927469e-03,
                                                   >> 115   2.801361e-03, 2.680213e-03, 2.563852e-03, 2.452110e-03, 2.344824e-03
                                                   >> 116 };
                                                   >> 117  
                                                   >> 118 ///////////////////////////////////////////////////////////////////////    
                                                   >> 119 //
 92 //  Constructor                                   120 //  Constructor
 93 G4SynchrotronRadiationInMat::G4SynchrotronRadi << 121 //
 94   const G4String& processName, G4ProcessType t << 122 
 95   : G4VDiscreteProcess(processName, type)      << 123 G4SynchrotronRadiationInMat::G4SynchrotronRadiationInMat(const G4String& processName,
 96   , theGamma(G4Gamma::Gamma())                 << 124   G4ProcessType type):G4VDiscreteProcess (processName, type),
 97   , theElectron(G4Electron::Electron())        << 125   LowestKineticEnergy (10.*keV),
 98   , thePositron(G4Positron::Positron())        << 126   HighestKineticEnergy (100.*TeV),
 99   , LowestKineticEnergy(10. * keV)             << 127   TotBin(200),
100   , fAlpha(0.0)                                << 128   theGamma (G4Gamma::Gamma() ),
101   , fRootNumber(80)                            << 129   theElectron ( G4Electron::Electron() ),
102   , fVerboseLevel(verboseLevel)                << 130   thePositron ( G4Positron::Positron() ), fAlpha(0.0), fRootNumber(80),
                                                   >> 131   fVerboseLevel( verboseLevel )
103 {                                                 132 {
104   G4TransportationManager* transportMgr =      << 133   G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
105     G4TransportationManager::GetTransportation << 
106                                                   134 
107   fFieldPropagator = transportMgr->GetPropagat    135   fFieldPropagator = transportMgr->GetPropagatorInField();
108   secID = G4PhysicsModelCatalog::GetModelID("m << 
109   SetProcessSubType(fSynchrotronRadiation);    << 
110   CutInRange = GammaCutInKineticEnergyNow = El << 
111     PositronCutInKineticEnergyNow = ParticleCu << 
112       fPsiGamma = fEta = fOrderAngleK = 0.0;   << 
113 }                                              << 
114                                                   136 
                                                   >> 137 }
                                                   >> 138  
115 //////////////////////////////////////////////    139 /////////////////////////////////////////////////////////////////////////
                                                   >> 140 //
116 // Destructor                                     141 // Destructor
117 G4SynchrotronRadiationInMat::~G4SynchrotronRad << 142 //
118                                                << 143  
119 G4bool G4SynchrotronRadiationInMat::IsApplicab << 144 G4SynchrotronRadiationInMat::~G4SynchrotronRadiationInMat()
120   const G4ParticleDefinition& particle)        << 
121 {                                                 145 {
122   return ((&particle == (const G4ParticleDefin << 146      ;
123           (&particle == (const G4ParticleDefin << 
124 }                                                 147 }
                                                   >> 148  
                                                   >> 149  
                                                   >> 150 /////////////////////////////// METHODS /////////////////////////////////
                                                   >> 151 //
                                                   >> 152 //
                                                   >> 153 // Production of synchrotron X-ray photon
                                                   >> 154 // GEANT4 internal units.
                                                   >> 155 //
125                                                   156 
126 G4double G4SynchrotronRadiationInMat::GetLambd << 
127                                                << 
128 G4double G4SynchrotronRadiationInMat::GetEnerg << 
129                                                   157 
130 // Production of synchrotron X-ray photon      << 158 G4double 
131 // Geant4 internal units.                      << 159 G4SynchrotronRadiationInMat::GetMeanFreePath( const G4Track& trackData,
132 G4double G4SynchrotronRadiationInMat::GetMeanF << 160                                                G4double,
133   const G4Track& trackData, G4double, G4ForceC << 161                                                G4ForceCondition* condition)
134 {                                                 162 {
135   // gives the MeanFreePath in GEANT4 internal    163   // gives the MeanFreePath in GEANT4 internal units
136   G4double MeanFreePath;                          164   G4double MeanFreePath;
137                                                   165 
138   const G4DynamicParticle* aDynamicParticle =     166   const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
                                                   >> 167   // G4Material* aMaterial = trackData.GetMaterial();
139                                                   168 
140   *condition = NotForced;                      << 169   //G4bool isOutRange ;
                                                   >> 170  
                                                   >> 171   *condition = NotForced ;
141                                                   172 
142   G4double gamma =                             << 173   G4double gamma = aDynamicParticle->GetTotalEnergy()/
143     aDynamicParticle->GetTotalEnergy() / aDyna << 174                    aDynamicParticle->GetMass();
144                                                   175 
145   G4double particleCharge = aDynamicParticle->    176   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
146                                                   177 
147   G4double KineticEnergy = aDynamicParticle->G    178   G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
148                                                   179 
149   if(KineticEnergy < LowestKineticEnergy || ga << 180   if ( KineticEnergy <  LowestKineticEnergy || gamma < 1.0e3 )  MeanFreePath = DBL_MAX; 
150     MeanFreePath = DBL_MAX;                    << 
151   else                                            181   else
152   {                                               182   {
153     G4ThreeVector FieldValue;                  << 
154     const G4Field* pField = nullptr;           << 
155                                                   183 
156     G4FieldManager* fieldMgr = nullptr;        << 184     G4ThreeVector  FieldValue;
157     G4bool fieldExertsForce  = false;          << 185     const G4Field*   pField = 0;
                                                   >> 186 
                                                   >> 187     G4FieldManager* fieldMgr=0;
                                                   >> 188     G4bool          fieldExertsForce = false;
158                                                   189 
159     if((particleCharge != 0.0))                << 190     if( (particleCharge != 0.0) )
160     {                                             191     {
161       fieldMgr =                               << 192       fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
162         fFieldPropagator->FindAndSetFieldManag << 
163                                                   193 
164       if(fieldMgr != nullptr)                  << 194       if ( fieldMgr != 0 ) 
165       {                                           195       {
166         // If the field manager has no field,     196         // If the field manager has no field, there is no field !
167         fieldExertsForce = (fieldMgr->GetDetec << 197 
                                                   >> 198         fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
168       }                                           199       }
169     }                                             200     }
170     if(fieldExertsForce)                       << 201     if ( fieldExertsForce )
171     {                                          << 202     {  
172       pField                     = fieldMgr->G << 203       pField = fieldMgr->GetDetectorField() ;
173       G4ThreeVector globPosition = trackData.G << 204       G4ThreeVector  globPosition = trackData.GetPosition();
174                                                   205 
175       G4double globPosVec[4], FieldValueVec[6] << 206       G4double  globPosVec[3], FieldValueVec[3];
176                                                   207 
177       globPosVec[0] = globPosition.x();           208       globPosVec[0] = globPosition.x();
178       globPosVec[1] = globPosition.y();           209       globPosVec[1] = globPosition.y();
179       globPosVec[2] = globPosition.z();           210       globPosVec[2] = globPosition.z();
180       globPosVec[3] = trackData.GetGlobalTime( << 
181                                                   211 
182       pField->GetFieldValue(globPosVec, FieldV << 212       pField->GetFieldValue( globPosVec, FieldValueVec );
183                                                   213 
184       FieldValue =                             << 214       FieldValue = G4ThreeVector( FieldValueVec[0], 
185         G4ThreeVector(FieldValueVec[0], FieldV << 215                                    FieldValueVec[1], 
                                                   >> 216                                    FieldValueVec[2]   );
                                                   >> 217 
                                                   >> 218        
                                                   >> 219 
                                                   >> 220       G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
                                                   >> 221       G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum) ;
                                                   >> 222       G4double perpB             = unitMcrossB.mag() ;
                                                   >> 223       G4double beta              = aDynamicParticle->GetTotalMomentum()/
                                                   >> 224                                    (aDynamicParticle->GetTotalEnergy()    );
186                                                   225 
187       G4ThreeVector unitMomentum = aDynamicPar << 226       if( perpB > 0.0 ) MeanFreePath = fLambdaConst*beta/perpB;
188       G4ThreeVector unitMcrossB  = FieldValue. << 227       else              MeanFreePath = DBL_MAX;       
189       G4double perpB             = unitMcrossB << 
190       G4double beta              = aDynamicPar << 
191                       (aDynamicParticle->GetTo << 
192                                                << 
193       if(perpB > 0.0)                          << 
194         MeanFreePath = fLambdaConst * beta / p << 
195       else                                     << 
196         MeanFreePath = DBL_MAX;                << 
197     }                                             228     }
198     else                                       << 229     else  MeanFreePath = DBL_MAX;    
199       MeanFreePath = DBL_MAX;                  << 
200   }                                               230   }
201   if(fVerboseLevel > 0)                           231   if(fVerboseLevel > 0)
202   {                                               232   {
203     G4cout << "G4SynchrotronRadiationInMat::Me << 233     G4cout<<"G4SynchrotronRadiationInMat::MeanFreePath = "<<MeanFreePath/m<<" m"<<G4endl;
204            << " m" << G4endl;                  << 
205   }                                               234   }
206   return MeanFreePath;                         << 235   return MeanFreePath; 
207 }                                              << 236 } 
208                                                   237 
209 //////////////////////////////////////////////    238 ////////////////////////////////////////////////////////////////////////////////
210 G4VParticleChange* G4SynchrotronRadiationInMat << 239 //
211   const G4Track& trackData, const G4Step& step << 240 //
                                                   >> 241 
                                                   >> 242 G4VParticleChange*
                                                   >> 243 G4SynchrotronRadiationInMat::PostStepDoIt(const G4Track& trackData,
                                                   >> 244                                      const G4Step&  stepData      )
212                                                   245 
213 {                                                 246 {
214   aParticleChange.Initialize(trackData);          247   aParticleChange.Initialize(trackData);
215                                                   248 
216   const G4DynamicParticle* aDynamicParticle =  << 249   const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
217                                                   250 
218   G4double gamma =                             << 251   G4double gamma = aDynamicParticle->GetTotalEnergy()/
219     aDynamicParticle->GetTotalEnergy() / (aDyn << 252                    (aDynamicParticle->GetMass()              );
220                                                   253 
221   if(gamma <= 1.0e3)                           << 254   if(gamma <= 1.0e3 ) 
222   {                                               255   {
223     return G4VDiscreteProcess::PostStepDoIt(tr << 256     return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
224   }                                               257   }
225   G4double particleCharge = aDynamicParticle->    258   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
226                                                   259 
227   G4ThreeVector FieldValue;                    << 260   G4ThreeVector  FieldValue;
228   const G4Field* pField = nullptr;             << 261   const G4Field*   pField = 0 ;
229                                                   262 
230   G4FieldManager* fieldMgr = nullptr;          << 263   G4FieldManager* fieldMgr=0;
231   G4bool fieldExertsForce  = false;            << 264   G4bool          fieldExertsForce = false;
232                                                   265 
233   if((particleCharge != 0.0))                  << 266   if( (particleCharge != 0.0) )
234   {                                               267   {
235     fieldMgr = fFieldPropagator->FindAndSetFie << 268     fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
236     if(fieldMgr != nullptr)                    << 269     if ( fieldMgr != 0 ) 
237     {                                             270     {
238       // If the field manager has no field, th    271       // If the field manager has no field, there is no field !
239       fieldExertsForce = (fieldMgr->GetDetecto << 272 
                                                   >> 273       fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
240     }                                             274     }
241   }                                               275   }
242   if(fieldExertsForce)                         << 276   if ( fieldExertsForce )
243   {                                               277   {
244     pField                     = fieldMgr->Get << 278     pField = fieldMgr->GetDetectorField() ;
245     G4ThreeVector globPosition = trackData.Get << 279     G4ThreeVector  globPosition = trackData.GetPosition() ;
246     G4double globPosVec[4], FieldValueVec[6];  << 280     G4double  globPosVec[3], FieldValueVec[3] ;
247     globPosVec[0] = globPosition.x();          << 281     globPosVec[0] = globPosition.x() ;
248     globPosVec[1] = globPosition.y();          << 282     globPosVec[1] = globPosition.y() ;
249     globPosVec[2] = globPosition.z();          << 283     globPosVec[2] = globPosition.z() ;
250     globPosVec[3] = trackData.GetGlobalTime(); << 284 
251                                                << 285     pField->GetFieldValue( globPosVec, FieldValueVec ) ;
252     pField->GetFieldValue(globPosVec, FieldVal << 286     FieldValue = G4ThreeVector( FieldValueVec[0], 
253     FieldValue =                               << 287                                    FieldValueVec[1], 
254       G4ThreeVector(FieldValueVec[0], FieldVal << 288                                    FieldValueVec[2]   );
255                                                << 289        
256     G4ThreeVector unitMomentum = aDynamicParti << 290     G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
257     G4ThreeVector unitMcrossB  = FieldValue.cr << 291     G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
258     G4double perpB             = unitMcrossB.m << 292     G4double perpB = unitMcrossB.mag() ;
259     if(perpB > 0.0)                               293     if(perpB > 0.0)
260     {                                             294     {
261       // M-C of synchrotron photon energy         295       // M-C of synchrotron photon energy
262       G4double energyOfSR = GetRandomEnergySR( << 296 
                                                   >> 297       G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
263                                                   298 
264       if(fVerboseLevel > 0)                       299       if(fVerboseLevel > 0)
265       {                                           300       {
266         G4cout << "SR photon energy = " << ene << 301         G4cout<<"SR photon energy = "<<energyOfSR/keV<<" keV"<<G4endl;
267       }                                           302       }
268                                                << 
269       // check against insufficient energy        303       // check against insufficient energy
270       if(energyOfSR <= 0.0)                    << 
271       {                                        << 
272         return G4VDiscreteProcess::PostStepDoI << 
273       }                                        << 
274       G4double kineticEnergy = aDynamicParticl << 
275       G4ParticleMomentum particleDirection =   << 
276         aDynamicParticle->GetMomentumDirection << 
277                                                   304 
278       // M-C of its direction, simplified dipo << 305       if( energyOfSR <= 0.0 )
279       G4double cosTheta, sinTheta, fcos, beta; << 
280                                                << 
281       do                                       << 
282       {                                           306       {
283         cosTheta = 1. - 2. * G4UniformRand();  << 307         return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
284         fcos     = (1 + cosTheta * cosTheta) * << 
285       }                                           308       }
286       // Loop checking, 07-Aug-2015, Vladimir  << 309       G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
287       while(fcos < G4UniformRand());           << 310       G4ParticleMomentum 
288                                                << 311       particleDirection = aDynamicParticle->GetMomentumDirection();
289       beta = std::sqrt(1. - 1. / (gamma * gamm << 
290                                                << 
291       cosTheta = (cosTheta + beta) / (1. + bet << 
292                                                   312 
293       if(cosTheta > 1.)                        << 313       // M-C of its direction
294         cosTheta = 1.;                         << 314       
295       if(cosTheta < -1.)                       << 315       G4double Teta = G4UniformRand()/gamma ;    // Very roughly
296         cosTheta = -1.;                        << 
297                                                   316 
298       sinTheta = std::sqrt(1. - cosTheta * cos << 317       G4double Phi  = twopi * G4UniformRand() ;
299                                                   318 
300       G4double Phi = twopi * G4UniformRand();  << 319       G4double dirx = sin(Teta)*cos(Phi) , 
                                                   >> 320                diry = sin(Teta)*sin(Phi) , 
                                                   >> 321                dirz = cos(Teta) ;
301                                                   322 
302       G4double dirx = sinTheta * std::cos(Phi) << 323       G4ThreeVector gammaDirection ( dirx, diry, dirz);
303       G4double diry = sinTheta * std::sin(Phi) << 324       gammaDirection.rotateUz(particleDirection);   
304       G4double dirz = cosTheta;                << 325  
                                                   >> 326       // polarization of new gamma
305                                                   327 
306       G4ThreeVector gammaDirection(dirx, diry, << 328       // G4double sx =  cos(Teta)*cos(Phi);
307       gammaDirection.rotateUz(particleDirectio << 329       // G4double sy =  cos(Teta)*sin(Phi);
                                                   >> 330       // G4double sz = -sin(Teta);
308                                                   331 
309       // polarization of new gamma             << 
310       G4ThreeVector gammaPolarization = FieldV    332       G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
311       gammaPolarization               = gammaP << 333       gammaPolarization = gammaPolarization.unit();
                                                   >> 334 
                                                   >> 335       // (sx, sy, sz);
                                                   >> 336       // gammaPolarization.rotateUz(particleDirection);
312                                                   337 
313       // create G4DynamicParticle object for t    338       // create G4DynamicParticle object for the SR photon
314       auto aGamma =                            << 339  
315         new G4DynamicParticle(G4Gamma::Gamma() << 340       G4DynamicParticle* aGamma= new G4DynamicParticle ( G4Gamma::Gamma(),
316       aGamma->SetPolarization(gammaPolarizatio << 341                                                          gammaDirection, 
317                               gammaPolarizatio << 342                                                          energyOfSR      );
                                                   >> 343       aGamma->SetPolarization( gammaPolarization.x(),
                                                   >> 344                                gammaPolarization.y(),
                                                   >> 345                                gammaPolarization.z() );
318                                                   346 
319       aParticleChange.SetNumberOfSecondaries(1 << 
320                                                   347 
321       // Update the incident particle          << 348       aParticleChange.SetNumberOfSecondaries(1);
322       G4double newKinEnergy = kineticEnergy -  << 349       aParticleChange.AddSecondary(aGamma); 
323                                                   350 
324       if(newKinEnergy > 0.)                    << 351       // Update the incident particle 
                                                   >> 352    
                                                   >> 353       G4double newKinEnergy = kineticEnergy - energyOfSR ;
                                                   >> 354       
                                                   >> 355       if (newKinEnergy > 0.)
325       {                                           356       {
326         aParticleChange.ProposeMomentumDirecti << 357         aParticleChange.ProposeMomentumDirection( particleDirection );
327         aParticleChange.ProposeEnergy(newKinEn << 358         aParticleChange.ProposeEnergy( newKinEnergy );
328         aParticleChange.ProposeLocalEnergyDepo << 359         aParticleChange.ProposeLocalEnergyDeposit (0.); 
329       }                                        << 360       } 
330       else                                        361       else
331       {                                        << 362       { 
332         aParticleChange.ProposeEnergy(0.);     << 363         aParticleChange.ProposeEnergy( 0. );
333         aParticleChange.ProposeLocalEnergyDepo << 364         aParticleChange.ProposeLocalEnergyDeposit (0.);
334         G4double charge = aDynamicParticle->Ge << 365         G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();   
335         if(charge < 0.)                        << 366         if (charge<0.)
336         {                                         367         {
337           aParticleChange.ProposeTrackStatus(f << 368            aParticleChange.ProposeTrackStatus(fStopAndKill) ;
338         }                                      << 369   }
339         else                                   << 370         else      
340         {                                         371         {
341           aParticleChange.ProposeTrackStatus(f << 372           aParticleChange.ProposeTrackStatus(fStopButAlive) ;
342         }                                      << 373   }
343       }                                        << 374       } 
344                                                << 375     }       
345       // Create the G4Track                    << 
346       G4Track* aSecondaryTrack = new G4Track(a << 
347       aSecondaryTrack->SetTouchableHandle(step << 
348       aSecondaryTrack->SetParentID(trackData.G << 
349       aSecondaryTrack->SetCreatorModelID(secID << 
350       aParticleChange.AddSecondary(aSecondaryT << 
351     }                                          << 
352     else                                          376     else
353     {                                             377     {
354       return G4VDiscreteProcess::PostStepDoIt( << 378       return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
355     }                                             379     }
356   }                                            << 380   }     
357   return G4VDiscreteProcess::PostStepDoIt(trac << 381   return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
358 }                                                 382 }
359                                                   383 
360 G4double G4SynchrotronRadiationInMat::GetPhoto << 384 
361                                                << 385 G4double 
                                                   >> 386 G4SynchrotronRadiationInMat::GetPhotonEnergy( const G4Track& trackData,
                                                   >> 387                                          const G4Step&  )
                                                   >> 388 
362 {                                                 389 {
363   G4int i;                                     << 390   G4int i ;
364   G4double energyOfSR = -1.0;                  << 391   G4double energyOfSR = -1.0 ;
                                                   >> 392   //G4Material* aMaterial=trackData.GetMaterial() ;
365                                                   393 
366   const G4DynamicParticle* aDynamicParticle =  << 394   const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
367                                                   395 
368   G4double gamma =                             << 396   G4double gamma = aDynamicParticle->GetTotalEnergy()/
369     aDynamicParticle->GetTotalEnergy() / (aDyn << 397                    (aDynamicParticle->GetMass()              ) ;
370                                                   398 
371   G4double particleCharge = aDynamicParticle->    399   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
372                                                   400 
373   G4ThreeVector FieldValue;                    << 401   G4ThreeVector  FieldValue;
374   const G4Field* pField = nullptr;             << 402   const G4Field*   pField = 0 ;
375                                                   403 
376   G4FieldManager* fieldMgr = nullptr;          << 404   G4FieldManager* fieldMgr=0;
377   G4bool fieldExertsForce  = false;            << 405   G4bool          fieldExertsForce = false;
378                                                   406 
379   if((particleCharge != 0.0))                  << 407   if( (particleCharge != 0.0) )
380   {                                               408   {
381     fieldMgr = fFieldPropagator->FindAndSetFie << 409     fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
382     if(fieldMgr != nullptr)                    << 410     if ( fieldMgr != 0 ) 
383     {                                             411     {
384       // If the field manager has no field, th    412       // If the field manager has no field, there is no field !
385       fieldExertsForce = (fieldMgr->GetDetecto << 413 
                                                   >> 414       fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
386     }                                             415     }
387   }                                               416   }
388   if(fieldExertsForce)                         << 417   if ( fieldExertsForce )
389   {                                            << 418   {  
390     pField                     = fieldMgr->Get << 419     pField = fieldMgr->GetDetectorField();
391     G4ThreeVector globPosition = trackData.Get << 420     G4ThreeVector  globPosition = trackData.GetPosition();
392     G4double globPosVec[3], FieldValueVec[3];  << 421     G4double  globPosVec[3], FieldValueVec[3];
393                                                   422 
394     globPosVec[0] = globPosition.x();             423     globPosVec[0] = globPosition.x();
395     globPosVec[1] = globPosition.y();             424     globPosVec[1] = globPosition.y();
396     globPosVec[2] = globPosition.z();             425     globPosVec[2] = globPosition.z();
397                                                   426 
398     pField->GetFieldValue(globPosVec, FieldVal << 427     pField->GetFieldValue( globPosVec, FieldValueVec );
399     FieldValue =                               << 428     FieldValue = G4ThreeVector( FieldValueVec[0], 
400       G4ThreeVector(FieldValueVec[0], FieldVal << 429                                    FieldValueVec[1], 
401                                                << 430                                    FieldValueVec[2]   );
402     G4ThreeVector unitMomentum = aDynamicParti << 431        
403     G4ThreeVector unitMcrossB  = FieldValue.cr << 432     G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
                                                   >> 433     G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum) ;
404     G4double perpB             = unitMcrossB.m    434     G4double perpB             = unitMcrossB.mag();
405     if(perpB > 0.0)                            << 435     if( perpB > 0.0 )
406     {                                             436     {
407       // M-C of synchrotron photon energy         437       // M-C of synchrotron photon energy
408       G4double random = G4UniformRand();       << 438 
409       for(i = 0; i < 200; ++i)                 << 439       G4double random = G4UniformRand() ;
                                                   >> 440       for(i=0;i<200;i++)
410       {                                           441       {
411         if(random >= fIntegralProbabilityOfSR[ << 442         if(random >= fIntegralProbabilityOfSR[i]) break ;
412           break;                               << 
413       }                                           443       }
414       energyOfSR = 0.0001 * i * i * fEnergyCon << 444       energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
415                                                   445 
416       // check against insufficient energy        446       // check against insufficient energy
                                                   >> 447 
417       if(energyOfSR <= 0.0)                       448       if(energyOfSR <= 0.0)
418       {                                           449       {
419         return -1.0;                           << 450         return -1.0 ;
420       }                                           451       }
421     }                                          << 452       //G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
                                                   >> 453       //G4ParticleMomentum 
                                                   >> 454       //particleDirection = aDynamicParticle->GetMomentumDirection();
                                                   >> 455 
                                                   >> 456       // Gamma production cut in this material
                                                   >> 457       //G4double 
                                                   >> 458       //gammaEnergyCut = (G4Gamma::GetCutsInEnergy())[aMaterial->GetIndex()];
                                                   >> 459 
                                                   >> 460       // SR photon has energy more than the current material cut
                                                   >> 461       // M-C of its direction
                                                   >> 462       
                                                   >> 463       //G4double Teta = G4UniformRand()/gamma ;    // Very roughly
                                                   >> 464 
                                                   >> 465       //G4double Phi  = twopi * G4UniformRand() ;
                                                   >> 466     }       
422     else                                          467     else
423     {                                             468     {
424       return -1.0;                             << 469       return -1.0 ;
425     }                                             470     }
426   }                                            << 471   }     
427   return energyOfSR;                           << 472   return energyOfSR ;
428 }                                                 473 }
429                                                   474 
430 //////////////////////////////////////////////    475 /////////////////////////////////////////////////////////////////////////////////
431 G4double G4SynchrotronRadiationInMat::GetRando << 476 //
432                                                << 477 //
                                                   >> 478 
                                                   >> 479 G4double G4SynchrotronRadiationInMat::GetRandomEnergySR(G4double gamma, G4double perpB)
433 {                                                 480 {
434   G4int i;                                     << 481   G4int i, iMax;
435   static constexpr G4int iMax = 200;           << 
436   G4double energySR, random, position;            482   G4double energySR, random, position;
437                                                   483 
                                                   >> 484   iMax = 200;
438   random = G4UniformRand();                       485   random = G4UniformRand();
439                                                   486 
440   for(i = 0; i < iMax; ++i)                    << 487   for( i = 0; i < iMax; i++ )
441   {                                               488   {
442     if(random >= fIntegralProbabilityOfSR[i])  << 489     if( random >= fIntegralProbabilityOfSR[i] ) break;
443       break;                                   << 
444   }                                               490   }
445   if(i <= 0)                                   << 491   if(i <= 0 )        position = G4UniformRand(); // 0.
446     position = G4UniformRand();                << 492   else if( i>= iMax) position = G4double(iMax);
447   else if(i >= iMax)                           << 493   else position = i + G4UniformRand(); // -1
448     position = G4double(iMax);                 << 494   // 
449   else                                         << 495   // it was in initial implementation:
450     position = i + G4UniformRand();            << 496   //       energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
451                                                   497 
452   energySR =                                   << 498   energySR = 0.0001*position*position*fEnergyConst*gamma*gamma*perpB;
453     0.0001 * position * position * fEnergyCons << 
454                                                   499 
455   if(energySR < 0.)                            << 500   if( energySR  < 0. ) energySR = 0.;
456     energySR = 0.;                             << 
457                                                   501 
458   return energySR;                                502   return energySR;
459 }                                                 503 }
460                                                   504 
461 //////////////////////////////////////////////    505 /////////////////////////////////////////////////////////////////////////
462 G4double G4SynchrotronRadiationInMat::GetProbS << 506 //
                                                   >> 507 // return
                                                   >> 508 
                                                   >> 509 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt( G4double t)
463 {                                                 510 {
464   G4double result, hypCos2, hypCos = std::cosh << 511   G4double result, hypCos2, hypCos=std::cosh(t);
465                                                   512 
466   hypCos2 = hypCos * hypCos;                   << 513   hypCos2 = hypCos*hypCos;  
467   result = std::cosh(5. * t / 3.) * std::exp(t << 514   result = std::cosh(5.*t/3.)*std::exp(t-fKsi*hypCos); // fKsi > 0. !
468   result /= hypCos2;                              515   result /= hypCos2;
469   return result;                                  516   return result;
470 }                                                 517 }
471                                                   518 
472 //////////////////////////////////////////////    519 ///////////////////////////////////////////////////////////////////////////
                                                   >> 520 //
473 // return the probability to emit SR photon wi    521 // return the probability to emit SR photon with relative energy
474 // energy/energy_c >= ksi                         522 // energy/energy_c >= ksi
475 // for ksi <= 0. P = 1., however the method wo    523 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
476 G4double G4SynchrotronRadiationInMat::GetIntPr << 524 
                                                   >> 525 G4double G4SynchrotronRadiationInMat::GetIntProbSR( G4double ksi)
477 {                                                 526 {
478   if(ksi <= 0.)                                << 527   if (ksi <= 0.) return 1.0;
479     return 1.0;                                << 528   fKsi = ksi; // should be > 0. !
480   fKsi = ksi;  // should be > 0. !             << 
481   G4int n;                                        529   G4int n;
482   G4double result, a;                             530   G4double result, a;
483                                                   531 
484   a = fAlpha;       // always = 0.             << 532   a = fAlpha;        // always = 0.
485   n = fRootNumber;  // around default = 80     << 533   n = fRootNumber;   // around default = 80
486                                                   534 
487   G4Integrator<G4SynchrotronRadiationInMat,    << 535   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
488                G4double (G4SynchrotronRadiatio << 
489     integral;                                  << 
490                                                   536 
491   result = integral.Laguerre(                  << 537   result = integral.Laguerre(this, 
492     this, &G4SynchrotronRadiationInMat::GetPro << 538            &G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt, a, n);
493                                                   539 
494   result *= 3. / 5. / pi;                      << 540   result *= 3./5./pi;
495                                                   541 
496   return result;                                  542   return result;
497 }                                                 543 }
498                                                   544 
499 //////////////////////////////////////////////    545 /////////////////////////////////////////////////////////////////////////
                                                   >> 546 //
500 // return an auxiliary function for K_5/3 inte    547 // return an auxiliary function for K_5/3 integral representation
501 G4double G4SynchrotronRadiationInMat::GetProbS << 
502 {                                              << 
503   G4double result, hypCos = std::cosh(t);      << 
504                                                   548 
505   result = std::cosh(5. * t / 3.) * std::exp(t << 549 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy( G4double t)
                                                   >> 550 {
                                                   >> 551   G4double result, hypCos=std::cosh(t);
                                                   >> 552   
                                                   >> 553   result = std::cosh(5.*t/3.)*std::exp(t - fKsi*hypCos); // fKsi > 0. !
506   result /= hypCos;                               554   result /= hypCos;
507   return result;                                  555   return result;
508 }                                                 556 }
509                                                   557 
510 //////////////////////////////////////////////    558 ///////////////////////////////////////////////////////////////////////////
                                                   >> 559 //
511 // return the probability to emit SR photon en    560 // return the probability to emit SR photon energy with relative energy
512 // energy/energy_c >= ksi                         561 // energy/energy_c >= ksi
513 // for ksi <= 0. P = 1., however the method wo    562 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
514 G4double G4SynchrotronRadiationInMat::GetEnerg << 563 
                                                   >> 564 G4double G4SynchrotronRadiationInMat::GetEnergyProbSR( G4double ksi)
515 {                                                 565 {
516   if(ksi <= 0.)                                << 566   if (ksi <= 0.) return 1.0;
517     return 1.0;                                << 567   fKsi = ksi; // should be > 0. !
518   fKsi = ksi;  // should be > 0. !             << 
519   G4int n;                                        568   G4int n;
520   G4double result, a;                             569   G4double result, a;
521                                                   570 
522   a = fAlpha;       // always = 0.             << 571   a = fAlpha;        // always = 0.
523   n = fRootNumber;  // around default = 80     << 572   n = fRootNumber;   // around default = 80
524                                                   573 
525   G4Integrator<G4SynchrotronRadiationInMat,    << 574   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
526                G4double (G4SynchrotronRadiatio << 
527     integral;                                  << 
528                                                   575 
529   result = integral.Laguerre(                  << 576   result = integral.Laguerre(this, 
530     this, &G4SynchrotronRadiationInMat::GetPro << 577            &G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy, a, n);
531                                                   578 
532   result *= 9. * std::sqrt(3.) * ksi / 8. / pi << 579   result *= 9.*std::sqrt(3.)*ksi/8./pi;
533                                                   580 
534   return result;                                  581   return result;
535 }                                                 582 }
536                                                   583 
537 //////////////////////////////////////////////    584 /////////////////////////////////////////////////////////////////////////////
538 G4double G4SynchrotronRadiationInMat::GetInteg << 585 //
539 {                                              << 586 // 
540   G4double result, hypCos = std::cosh(t);      << 
541                                                   587 
542   result =                                     << 588 G4double G4SynchrotronRadiationInMat::GetIntegrandForAngleK( G4double t)
543     std::cosh(fOrderAngleK * t) * std::exp(t - << 589 {
                                                   >> 590   G4double result, hypCos=std::cosh(t);
                                                   >> 591   
                                                   >> 592   result  = std::cosh(fOrderAngleK*t)*std::exp(t - fEta*hypCos); // fEta > 0. !
544   result /= hypCos;                               593   result /= hypCos;
545   return result;                                  594   return result;
546 }                                                 595 }
547                                                   596 
548 //////////////////////////////////////////////    597 //////////////////////////////////////////////////////////////////////////
                                                   >> 598 //
549 // Return K 1/3 or 2/3 for angular distributio    599 // Return K 1/3 or 2/3 for angular distribution
550 G4double G4SynchrotronRadiationInMat::GetAngle << 600 
                                                   >> 601 G4double G4SynchrotronRadiationInMat::GetAngleK( G4double eta)
551 {                                                 602 {
552   fEta = eta;  // should be > 0. !             << 603   fEta = eta; // should be > 0. !
553   G4int n;                                        604   G4int n;
554   G4double result, a;                             605   G4double result, a;
555                                                   606 
556   a = fAlpha;       // always = 0.             << 607   a = fAlpha;        // always = 0.
557   n = fRootNumber;  // around default = 80     << 608   n = fRootNumber;   // around default = 80
558                                                   609 
559   G4Integrator<G4SynchrotronRadiationInMat,    << 610   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
560                G4double (G4SynchrotronRadiatio << 
561     integral;                                  << 
562                                                   611 
563   result = integral.Laguerre(                  << 612   result = integral.Laguerre(this, 
564     this, &G4SynchrotronRadiationInMat::GetInt << 613            &G4SynchrotronRadiationInMat::GetIntegrandForAngleK, a, n);
565                                                   614 
566   return result;                                  615   return result;
567 }                                                 616 }
568                                                   617 
569 //////////////////////////////////////////////    618 /////////////////////////////////////////////////////////////////////////
                                                   >> 619 //
570 // Relative angle diff distribution for given     620 // Relative angle diff distribution for given fKsi, which is set externally
571 G4double G4SynchrotronRadiationInMat::GetAngle << 
572 {                                              << 
573   G4double result, funK, funK2, gpsi2 = gpsi * << 
574                                                << 
575   fPsiGamma = gpsi;                            << 
576   fEta      = 0.5 * fKsi * (1. + gpsi2) * std: << 
577                                                   621 
578   fOrderAngleK = 1. / 3.;                      << 622 G4double G4SynchrotronRadiationInMat::GetAngleNumberAtGammaKsi( G4double gpsi)
579   funK         = GetAngleK(fEta);              << 623 {
580   funK2        = funK * funK;                  << 624   G4double result, funK, funK2, gpsi2 = gpsi*gpsi;
581                                                   625 
582   result = gpsi2 * funK2 / (1. + gpsi2);       << 626   fPsiGamma    = gpsi;
                                                   >> 627   fEta         = 0.5*fKsi*(1 + gpsi2)*std::sqrt(1 + gpsi2);
                                                   >> 628  
                                                   >> 629   fOrderAngleK = 1./3.;
                                                   >> 630   funK         = GetAngleK(fEta); 
                                                   >> 631   funK2        = funK*funK;
                                                   >> 632 
                                                   >> 633   result       = gpsi2*funK2/(1 + gpsi2);
                                                   >> 634 
                                                   >> 635   fOrderAngleK = 2./3.;
                                                   >> 636   funK         = GetAngleK(fEta); 
                                                   >> 637   funK2        = funK*funK;
                                                   >> 638 
                                                   >> 639   result      += funK2;
                                                   >> 640   result      *= (1 + gpsi2)*fKsi;
                                                   >> 641      
                                                   >> 642   return result;
                                                   >> 643 }
583                                                   644 
584   fOrderAngleK = 2. / 3.;                      << 
585   funK         = GetAngleK(fEta);              << 
586   funK2        = funK * funK;                  << 
587                                                   645 
588   result += funK2;                             << 646 ///////////////////// end of G4SynchrotronRadiationInMat.cc
589   result *= (1. + gpsi2) * fKsi;               << 
590                                                   647 
591   return result;                               << 
592 }                                              << 
593                                                   648