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 10.1.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 68037 2013-03-13 14:15:08Z gcosmo $
                                                   >>  28 //
 26 // -------------------------------------------     29 // --------------------------------------------------------------
 27 //      GEANT 4 class implementation file          30 //      GEANT 4 class implementation file
                                                   >>  31 //      CERN Geneva Switzerland
 28 //                                                 32 //
 29 //      History: first implementation,             33 //      History: first implementation,
 30 //      21-5-98 V.Grichine                         34 //      21-5-98 V.Grichine
 31 //      28-05-01, V.Ivanchenko minor changes t     35 //      28-05-01, V.Ivanchenko minor changes to provide ANSI -wall compilation
 32 //      04.03.05, V.Grichine: get local field      36 //      04.03.05, V.Grichine: get local field interface
 33 //      19-05-06, V.Ivanchenko rename from G4S     37 //      19-05-06, V.Ivanchenko rename from G4SynchrotronRadiation
 34 //                                                 38 //
                                                   >>  39 //
 35 //////////////////////////////////////////////     40 ///////////////////////////////////////////////////////////////////////////
 36                                                    41 
 37 #include "G4SynchrotronRadiationInMat.hh"          42 #include "G4SynchrotronRadiationInMat.hh"
 38                                                << 
 39 #include "G4EmProcessSubType.hh"               << 
 40 #include "G4Field.hh"                          << 
 41 #include "G4FieldManager.hh"                   << 
 42 #include "G4Integrator.hh"                     << 
 43 #include "G4PhysicalConstants.hh"                  43 #include "G4PhysicalConstants.hh"
 44 #include "G4PropagatorInField.hh"              << 
 45 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 46 #include "G4PhysicsModelCatalog.hh"            <<  45 #include "G4Integrator.hh"
                                                   >>  46 #include "G4EmProcessSubType.hh"
 47                                                    47 
 48 const G4double G4SynchrotronRadiationInMat::fI <<  48 ////////////////////////////////////////////////////////////////////
 49   1.000000e+00, 9.428859e-01, 9.094095e-01, 8. <<  49 //
 50   8.337008e-01, 8.124961e-01, 7.925217e-01, 7. <<  50 // Constant for calculation of mean free path
 51   7.381233e-01, 7.214521e-01, 7.053634e-01, 6. <<  51 //
 52   6.600922e-01, 6.458793e-01, 6.320533e-01, 6. <<  52 
 53   5.926459e-01, 5.801347e-01, 5.679103e-01, 5. <<  53 const G4double
 54   5.328395e-01, 5.216482e-01, 5.106904e-01, 4. <<  54 G4SynchrotronRadiationInMat::fLambdaConst = std::sqrt(3.0)*electron_mass_c2/
 55   4.791351e-01, 4.690316e-01, 4.591249e-01, 4. <<  55                                        (2.5*fine_structure_const*eplus*c_light) ;
 56   4.305320e-01, 4.213608e-01, 4.123623e-01, 4. <<  56 
 57   3.863639e-01, 3.780179e-01, 3.698262e-01, 3. <<  57 /////////////////////////////////////////////////////////////////////
 58   3.461460e-01, 3.385411e-01, 3.310757e-01, 3. <<  58 //
 59   3.094921e-01, 3.025605e-01, 2.957566e-01, 2. <<  59 // Constant for calculation of characterictic energy
 60   2.760907e-01, 2.697773e-01, 2.635817e-01, 2. <<  60 //
 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                                                    61 
 91 ////////////////////////////////////////////// <<  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() ), 
                                                   >> 131   GammaCutInKineticEnergy(0),
                                                   >> 132   ElectronCutInKineticEnergy(0),
                                                   >> 133   PositronCutInKineticEnergy(0),
                                                   >> 134   ParticleCutInKineticEnergy(0),
                                                   >> 135   fAlpha(0.0), fRootNumber(80),
                                                   >> 136   fVerboseLevel( verboseLevel )
103 {                                                 137 {
104   G4TransportationManager* transportMgr =      << 138   G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
105     G4TransportationManager::GetTransportation << 
106                                                   139 
107   fFieldPropagator = transportMgr->GetPropagat    140   fFieldPropagator = transportMgr->GetPropagatorInField();
108   secID = G4PhysicsModelCatalog::GetModelID("m << 
109   SetProcessSubType(fSynchrotronRadiation);       141   SetProcessSubType(fSynchrotronRadiation);
110   CutInRange = GammaCutInKineticEnergyNow = El << 142   CutInRange = GammaCutInKineticEnergyNow = ElectronCutInKineticEnergyNow = 
111     PositronCutInKineticEnergyNow = ParticleCu << 143     PositronCutInKineticEnergyNow =  ParticleCutInKineticEnergyNow = fKsi = 
112       fPsiGamma = fEta = fOrderAngleK = 0.0;   << 144     fPsiGamma = fEta = fOrderAngleK = 0.0;
113 }                                                 145 }
114                                                << 146  
115 //////////////////////////////////////////////    147 /////////////////////////////////////////////////////////////////////////
                                                   >> 148 //
116 // Destructor                                     149 // Destructor
117 G4SynchrotronRadiationInMat::~G4SynchrotronRad << 150 //
                                                   >> 151  
                                                   >> 152 G4SynchrotronRadiationInMat::~G4SynchrotronRadiationInMat()
                                                   >> 153 {}
                                                   >> 154  
118                                                   155 
119 G4bool G4SynchrotronRadiationInMat::IsApplicab << 156 G4bool
120   const G4ParticleDefinition& particle)        << 157 G4SynchrotronRadiationInMat::IsApplicable( const G4ParticleDefinition& particle )
121 {                                                 158 {
122   return ((&particle == (const G4ParticleDefin << 
123           (&particle == (const G4ParticleDefin << 
124 }                                              << 
125                                                   159 
126 G4double G4SynchrotronRadiationInMat::GetLambd << 160   return ( ( &particle == (const G4ParticleDefinition *)theElectron ) ||
127                                                << 161      ( &particle == (const G4ParticleDefinition *)thePositron )    );
128 G4double G4SynchrotronRadiationInMat::GetEnerg << 
129                                                   162 
                                                   >> 163 }
                                                   >> 164  
                                                   >> 165 /////////////////////////////// METHODS /////////////////////////////////
                                                   >> 166 //
                                                   >> 167 //
130 // Production of synchrotron X-ray photon         168 // Production of synchrotron X-ray photon
131 // Geant4 internal units.                      << 169 // GEANT4 internal units.
132 G4double G4SynchrotronRadiationInMat::GetMeanF << 170 //
133   const G4Track& trackData, G4double, G4ForceC << 171 
                                                   >> 172 
                                                   >> 173 G4double 
                                                   >> 174 G4SynchrotronRadiationInMat::GetMeanFreePath( const G4Track& trackData,
                                                   >> 175                                                G4double,
                                                   >> 176                                                G4ForceCondition* condition)
134 {                                                 177 {
135   // gives the MeanFreePath in GEANT4 internal    178   // gives the MeanFreePath in GEANT4 internal units
136   G4double MeanFreePath;                          179   G4double MeanFreePath;
137                                                   180 
138   const G4DynamicParticle* aDynamicParticle =     181   const G4DynamicParticle* aDynamicParticle = trackData.GetDynamicParticle();
                                                   >> 182   // G4Material* aMaterial = trackData.GetMaterial();
139                                                   183 
140   *condition = NotForced;                      << 184   //G4bool isOutRange ;
                                                   >> 185  
                                                   >> 186   *condition = NotForced ;
141                                                   187 
142   G4double gamma =                             << 188   G4double gamma = aDynamicParticle->GetTotalEnergy()/
143     aDynamicParticle->GetTotalEnergy() / aDyna << 189                    aDynamicParticle->GetMass();
144                                                   190 
145   G4double particleCharge = aDynamicParticle->    191   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
146                                                   192 
147   G4double KineticEnergy = aDynamicParticle->G    193   G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
148                                                   194 
149   if(KineticEnergy < LowestKineticEnergy || ga << 195   if ( KineticEnergy <  LowestKineticEnergy || gamma < 1.0e3 )  MeanFreePath = DBL_MAX; 
150     MeanFreePath = DBL_MAX;                    << 
151   else                                            196   else
152   {                                               197   {
153     G4ThreeVector FieldValue;                  << 
154     const G4Field* pField = nullptr;           << 
155                                                   198 
156     G4FieldManager* fieldMgr = nullptr;        << 199     G4ThreeVector  FieldValue;
157     G4bool fieldExertsForce  = false;          << 200     const G4Field*   pField = 0;
                                                   >> 201 
                                                   >> 202     G4FieldManager* fieldMgr=0;
                                                   >> 203     G4bool          fieldExertsForce = false;
158                                                   204 
159     if((particleCharge != 0.0))                << 205     if( (particleCharge != 0.0) )
160     {                                             206     {
161       fieldMgr =                               << 207       fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
162         fFieldPropagator->FindAndSetFieldManag << 
163                                                   208 
164       if(fieldMgr != nullptr)                  << 209       if ( fieldMgr != 0 ) 
165       {                                           210       {
166         // If the field manager has no field,     211         // If the field manager has no field, there is no field !
167         fieldExertsForce = (fieldMgr->GetDetec << 212 
                                                   >> 213         fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
168       }                                           214       }
169     }                                             215     }
170     if(fieldExertsForce)                       << 216     if ( fieldExertsForce )
171     {                                          << 217     {  
172       pField                     = fieldMgr->G << 218       pField = fieldMgr->GetDetectorField() ;
173       G4ThreeVector globPosition = trackData.G << 219       G4ThreeVector  globPosition = trackData.GetPosition();
174                                                   220 
175       G4double globPosVec[4], FieldValueVec[6] << 221       G4double  globPosVec[4], FieldValueVec[6];
176                                                   222 
177       globPosVec[0] = globPosition.x();           223       globPosVec[0] = globPosition.x();
178       globPosVec[1] = globPosition.y();           224       globPosVec[1] = globPosition.y();
179       globPosVec[2] = globPosition.z();           225       globPosVec[2] = globPosition.z();
180       globPosVec[3] = trackData.GetGlobalTime(    226       globPosVec[3] = trackData.GetGlobalTime();
181                                                   227 
182       pField->GetFieldValue(globPosVec, FieldV << 228       pField->GetFieldValue( globPosVec, FieldValueVec );
183                                                   229 
184       FieldValue =                             << 230       FieldValue = G4ThreeVector( FieldValueVec[0], 
185         G4ThreeVector(FieldValueVec[0], FieldV << 231                                    FieldValueVec[1], 
                                                   >> 232                                    FieldValueVec[2]   );
                                                   >> 233 
                                                   >> 234        
                                                   >> 235 
                                                   >> 236       G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
                                                   >> 237       G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum) ;
                                                   >> 238       G4double perpB             = unitMcrossB.mag() ;
                                                   >> 239       G4double beta              = aDynamicParticle->GetTotalMomentum()/
                                                   >> 240                                    (aDynamicParticle->GetTotalEnergy()    );
186                                                   241 
187       G4ThreeVector unitMomentum = aDynamicPar << 242       if( perpB > 0.0 ) MeanFreePath = fLambdaConst*beta/perpB;
188       G4ThreeVector unitMcrossB  = FieldValue. << 243       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     }                                             244     }
198     else                                       << 245     else  MeanFreePath = DBL_MAX;    
199       MeanFreePath = DBL_MAX;                  << 
200   }                                               246   }
201   if(fVerboseLevel > 0)                           247   if(fVerboseLevel > 0)
202   {                                               248   {
203     G4cout << "G4SynchrotronRadiationInMat::Me << 249     G4cout<<"G4SynchrotronRadiationInMat::MeanFreePath = "<<MeanFreePath/m<<" m"<<G4endl;
204            << " m" << G4endl;                  << 
205   }                                               250   }
206   return MeanFreePath;                         << 251   return MeanFreePath; 
207 }                                              << 252 } 
208                                                   253 
209 //////////////////////////////////////////////    254 ////////////////////////////////////////////////////////////////////////////////
210 G4VParticleChange* G4SynchrotronRadiationInMat << 255 //
211   const G4Track& trackData, const G4Step& step << 256 //
                                                   >> 257 
                                                   >> 258 G4VParticleChange*
                                                   >> 259 G4SynchrotronRadiationInMat::PostStepDoIt(const G4Track& trackData,
                                                   >> 260                                      const G4Step&  stepData      )
212                                                   261 
213 {                                                 262 {
214   aParticleChange.Initialize(trackData);          263   aParticleChange.Initialize(trackData);
215                                                   264 
216   const G4DynamicParticle* aDynamicParticle =  << 265   const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
217                                                   266 
218   G4double gamma =                             << 267   G4double gamma = aDynamicParticle->GetTotalEnergy()/
219     aDynamicParticle->GetTotalEnergy() / (aDyn << 268                    (aDynamicParticle->GetMass()              );
220                                                   269 
221   if(gamma <= 1.0e3)                           << 270   if(gamma <= 1.0e3 ) 
222   {                                               271   {
223     return G4VDiscreteProcess::PostStepDoIt(tr << 272     return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
224   }                                               273   }
225   G4double particleCharge = aDynamicParticle->    274   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
226                                                   275 
227   G4ThreeVector FieldValue;                    << 276   G4ThreeVector  FieldValue;
228   const G4Field* pField = nullptr;             << 277   const G4Field*   pField = 0 ;
229                                                   278 
230   G4FieldManager* fieldMgr = nullptr;          << 279   G4FieldManager* fieldMgr=0;
231   G4bool fieldExertsForce  = false;            << 280   G4bool          fieldExertsForce = false;
232                                                   281 
233   if((particleCharge != 0.0))                  << 282   if( (particleCharge != 0.0) )
234   {                                               283   {
235     fieldMgr = fFieldPropagator->FindAndSetFie << 284     fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
236     if(fieldMgr != nullptr)                    << 285     if ( fieldMgr != 0 ) 
237     {                                             286     {
238       // If the field manager has no field, th    287       // If the field manager has no field, there is no field !
239       fieldExertsForce = (fieldMgr->GetDetecto << 288 
                                                   >> 289       fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
240     }                                             290     }
241   }                                               291   }
242   if(fieldExertsForce)                         << 292   if ( fieldExertsForce )
243   {                                               293   {
244     pField                     = fieldMgr->Get << 294     pField = fieldMgr->GetDetectorField() ;
245     G4ThreeVector globPosition = trackData.Get << 295     G4ThreeVector  globPosition = trackData.GetPosition() ;
246     G4double globPosVec[4], FieldValueVec[6];  << 296     G4double  globPosVec[4], FieldValueVec[6] ;
247     globPosVec[0] = globPosition.x();          << 297     globPosVec[0] = globPosition.x() ;
248     globPosVec[1] = globPosition.y();          << 298     globPosVec[1] = globPosition.y() ;
249     globPosVec[2] = globPosition.z();          << 299     globPosVec[2] = globPosition.z() ;
250     globPosVec[3] = trackData.GetGlobalTime();    300     globPosVec[3] = trackData.GetGlobalTime();
251                                                   301 
252     pField->GetFieldValue(globPosVec, FieldVal << 302     pField->GetFieldValue( globPosVec, FieldValueVec ) ;
253     FieldValue =                               << 303     FieldValue = G4ThreeVector( FieldValueVec[0], 
254       G4ThreeVector(FieldValueVec[0], FieldVal << 304                                    FieldValueVec[1], 
255                                                << 305                                    FieldValueVec[2]   );
256     G4ThreeVector unitMomentum = aDynamicParti << 306        
257     G4ThreeVector unitMcrossB  = FieldValue.cr << 307     G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
258     G4double perpB             = unitMcrossB.m << 308     G4ThreeVector unitMcrossB = FieldValue.cross(unitMomentum);
                                                   >> 309     G4double perpB = unitMcrossB.mag() ;
259     if(perpB > 0.0)                               310     if(perpB > 0.0)
260     {                                             311     {
261       // M-C of synchrotron photon energy         312       // M-C of synchrotron photon energy
262       G4double energyOfSR = GetRandomEnergySR( << 313 
                                                   >> 314       G4double energyOfSR = GetRandomEnergySR(gamma,perpB);
263                                                   315 
264       if(fVerboseLevel > 0)                       316       if(fVerboseLevel > 0)
265       {                                           317       {
266         G4cout << "SR photon energy = " << ene << 318         G4cout<<"SR photon energy = "<<energyOfSR/keV<<" keV"<<G4endl;
267       }                                           319       }
268                                                << 
269       // check against insufficient energy        320       // check against insufficient energy
270       if(energyOfSR <= 0.0)                    << 321 
                                                   >> 322       if( energyOfSR <= 0.0 )
271       {                                           323       {
272         return G4VDiscreteProcess::PostStepDoI << 324         return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
273       }                                           325       }
274       G4double kineticEnergy = aDynamicParticl    326       G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
275       G4ParticleMomentum particleDirection =   << 327       G4ParticleMomentum 
276         aDynamicParticle->GetMomentumDirection << 328       particleDirection = aDynamicParticle->GetMomentumDirection();
277                                                   329 
278       // M-C of its direction, simplified dipo    330       // M-C of its direction, simplified dipole busted approach
279       G4double cosTheta, sinTheta, fcos, beta; << 331       
                                                   >> 332       // G4double Teta = G4UniformRand()/gamma ;    // Very roughly
280                                                   333 
281       do                                       << 334       G4double cosTheta, sinTheta, fcos, beta;
282       {                                        << 
283         cosTheta = 1. - 2. * G4UniformRand();  << 
284         fcos     = (1 + cosTheta * cosTheta) * << 
285       }                                        << 
286       // Loop checking, 07-Aug-2015, Vladimir  << 
287       while(fcos < G4UniformRand());           << 
288                                                   335 
289       beta = std::sqrt(1. - 1. / (gamma * gamm << 336   do
                                                   >> 337   { 
                                                   >> 338     cosTheta = 1. - 2.*G4UniformRand();
                                                   >> 339     fcos     = (1 + cosTheta*cosTheta)*0.5;
                                                   >> 340   }
                                                   >> 341   while( fcos < G4UniformRand() );
290                                                   342 
291       cosTheta = (cosTheta + beta) / (1. + bet << 343   beta = std::sqrt(1. - 1./(gamma*gamma));
292                                                   344 
293       if(cosTheta > 1.)                        << 345   cosTheta = (cosTheta + beta)/(1. + beta*cosTheta);
294         cosTheta = 1.;                         << 
295       if(cosTheta < -1.)                       << 
296         cosTheta = -1.;                        << 
297                                                   346 
298       sinTheta = std::sqrt(1. - cosTheta * cos << 347   if( cosTheta >  1. ) cosTheta =  1.;
                                                   >> 348   if( cosTheta < -1. ) cosTheta = -1.;
299                                                   349 
300       G4double Phi = twopi * G4UniformRand();  << 350   sinTheta = std::sqrt(1. - cosTheta*cosTheta );
301                                                   351 
302       G4double dirx = sinTheta * std::cos(Phi) << 352       G4double Phi  = twopi * G4UniformRand() ;
303       G4double diry = sinTheta * std::sin(Phi) << 
304       G4double dirz = cosTheta;                << 
305                                                   353 
306       G4ThreeVector gammaDirection(dirx, diry, << 354       G4double dirx = sinTheta*std::cos(Phi) , 
307       gammaDirection.rotateUz(particleDirectio << 355                diry = sinTheta*std::sin(Phi) , 
                                                   >> 356                dirz = cosTheta;
308                                                   357 
                                                   >> 358       G4ThreeVector gammaDirection ( dirx, diry, dirz);
                                                   >> 359       gammaDirection.rotateUz(particleDirection);   
                                                   >> 360  
309       // polarization of new gamma                361       // polarization of new gamma
                                                   >> 362 
                                                   >> 363       // G4double sx =  std::cos(Teta)*std::cos(Phi);
                                                   >> 364       // G4double sy =  std::cos(Teta)*std::sin(Phi);
                                                   >> 365       // G4double sz = -std::sin(Teta);
                                                   >> 366 
310       G4ThreeVector gammaPolarization = FieldV    367       G4ThreeVector gammaPolarization = FieldValue.cross(gammaDirection);
311       gammaPolarization               = gammaP << 368       gammaPolarization = gammaPolarization.unit();
                                                   >> 369 
                                                   >> 370       // (sx, sy, sz);
                                                   >> 371       // gammaPolarization.rotateUz(particleDirection);
312                                                   372 
313       // create G4DynamicParticle object for t    373       // create G4DynamicParticle object for the SR photon
314       auto aGamma =                            << 374  
315         new G4DynamicParticle(G4Gamma::Gamma() << 375       G4DynamicParticle* aGamma= new G4DynamicParticle ( G4Gamma::Gamma(),
316       aGamma->SetPolarization(gammaPolarizatio << 376                                                          gammaDirection, 
317                               gammaPolarizatio << 377                                                          energyOfSR      );
                                                   >> 378       aGamma->SetPolarization( gammaPolarization.x(),
                                                   >> 379                                gammaPolarization.y(),
                                                   >> 380                                gammaPolarization.z() );
318                                                   381 
319       aParticleChange.SetNumberOfSecondaries(1 << 
320                                                   382 
321       // Update the incident particle          << 383       aParticleChange.SetNumberOfSecondaries(1);
322       G4double newKinEnergy = kineticEnergy -  << 384       aParticleChange.AddSecondary(aGamma); 
323                                                   385 
324       if(newKinEnergy > 0.)                    << 386       // Update the incident particle 
                                                   >> 387    
                                                   >> 388       G4double newKinEnergy = kineticEnergy - energyOfSR ;
                                                   >> 389       
                                                   >> 390       if (newKinEnergy > 0.)
325       {                                           391       {
326         aParticleChange.ProposeMomentumDirecti << 392         aParticleChange.ProposeMomentumDirection( particleDirection );
327         aParticleChange.ProposeEnergy(newKinEn << 393         aParticleChange.ProposeEnergy( newKinEnergy );
328         aParticleChange.ProposeLocalEnergyDepo << 394         aParticleChange.ProposeLocalEnergyDeposit (0.); 
329       }                                        << 395       } 
330       else                                        396       else
331       {                                        << 397       { 
332         aParticleChange.ProposeEnergy(0.);     << 398         aParticleChange.ProposeEnergy( 0. );
333         aParticleChange.ProposeLocalEnergyDepo << 399         aParticleChange.ProposeLocalEnergyDeposit (0.);
334         G4double charge = aDynamicParticle->Ge << 400         G4double charge = aDynamicParticle->GetDefinition()->GetPDGCharge();   
335         if(charge < 0.)                        << 401         if (charge<0.)
336         {                                         402         {
337           aParticleChange.ProposeTrackStatus(f << 403            aParticleChange.ProposeTrackStatus(fStopAndKill) ;
338         }                                      << 404   }
339         else                                   << 405         else      
340         {                                         406         {
341           aParticleChange.ProposeTrackStatus(f << 407           aParticleChange.ProposeTrackStatus(fStopButAlive) ;
342         }                                      << 408   }
343       }                                        << 409       } 
344                                                << 410     }       
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                                          411     else
353     {                                             412     {
354       return G4VDiscreteProcess::PostStepDoIt( << 413       return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
355     }                                             414     }
356   }                                            << 415   }     
357   return G4VDiscreteProcess::PostStepDoIt(trac << 416   return G4VDiscreteProcess::PostStepDoIt(trackData,stepData);
358 }                                                 417 }
359                                                   418 
360 G4double G4SynchrotronRadiationInMat::GetPhoto << 419 
361                                                << 420 G4double 
                                                   >> 421 G4SynchrotronRadiationInMat::GetPhotonEnergy( const G4Track& trackData,
                                                   >> 422                                          const G4Step&  )
                                                   >> 423 
362 {                                                 424 {
363   G4int i;                                     << 425   G4int i ;
364   G4double energyOfSR = -1.0;                  << 426   G4double energyOfSR = -1.0 ;
                                                   >> 427   //G4Material* aMaterial=trackData.GetMaterial() ;
365                                                   428 
366   const G4DynamicParticle* aDynamicParticle =  << 429   const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
367                                                   430 
368   G4double gamma =                             << 431   G4double gamma = aDynamicParticle->GetTotalEnergy()/
369     aDynamicParticle->GetTotalEnergy() / (aDyn << 432                    (aDynamicParticle->GetMass()              ) ;
370                                                   433 
371   G4double particleCharge = aDynamicParticle->    434   G4double particleCharge = aDynamicParticle->GetDefinition()->GetPDGCharge();
372                                                   435 
373   G4ThreeVector FieldValue;                    << 436   G4ThreeVector  FieldValue;
374   const G4Field* pField = nullptr;             << 437   const G4Field*   pField = 0 ;
375                                                   438 
376   G4FieldManager* fieldMgr = nullptr;          << 439   G4FieldManager* fieldMgr=0;
377   G4bool fieldExertsForce  = false;            << 440   G4bool          fieldExertsForce = false;
378                                                   441 
379   if((particleCharge != 0.0))                  << 442   if( (particleCharge != 0.0) )
380   {                                               443   {
381     fieldMgr = fFieldPropagator->FindAndSetFie << 444     fieldMgr = fFieldPropagator->FindAndSetFieldManager( trackData.GetVolume() );
382     if(fieldMgr != nullptr)                    << 445     if ( fieldMgr != 0 ) 
383     {                                             446     {
384       // If the field manager has no field, th    447       // If the field manager has no field, there is no field !
385       fieldExertsForce = (fieldMgr->GetDetecto << 448 
                                                   >> 449       fieldExertsForce = ( fieldMgr->GetDetectorField() != 0 );
386     }                                             450     }
387   }                                               451   }
388   if(fieldExertsForce)                         << 452   if ( fieldExertsForce )
389   {                                            << 453   {  
390     pField                     = fieldMgr->Get << 454     pField = fieldMgr->GetDetectorField();
391     G4ThreeVector globPosition = trackData.Get << 455     G4ThreeVector  globPosition = trackData.GetPosition();
392     G4double globPosVec[3], FieldValueVec[3];  << 456     G4double  globPosVec[3], FieldValueVec[3];
393                                                   457 
394     globPosVec[0] = globPosition.x();             458     globPosVec[0] = globPosition.x();
395     globPosVec[1] = globPosition.y();             459     globPosVec[1] = globPosition.y();
396     globPosVec[2] = globPosition.z();             460     globPosVec[2] = globPosition.z();
397                                                   461 
398     pField->GetFieldValue(globPosVec, FieldVal << 462     pField->GetFieldValue( globPosVec, FieldValueVec );
399     FieldValue =                               << 463     FieldValue = G4ThreeVector( FieldValueVec[0], 
400       G4ThreeVector(FieldValueVec[0], FieldVal << 464                                    FieldValueVec[1], 
401                                                << 465                                    FieldValueVec[2]   );
402     G4ThreeVector unitMomentum = aDynamicParti << 466        
403     G4ThreeVector unitMcrossB  = FieldValue.cr << 467     G4ThreeVector unitMomentum = aDynamicParticle->GetMomentumDirection(); 
                                                   >> 468     G4ThreeVector unitMcrossB  = FieldValue.cross(unitMomentum) ;
404     G4double perpB             = unitMcrossB.m    469     G4double perpB             = unitMcrossB.mag();
405     if(perpB > 0.0)                            << 470     if( perpB > 0.0 )
406     {                                             471     {
407       // M-C of synchrotron photon energy         472       // M-C of synchrotron photon energy
408       G4double random = G4UniformRand();       << 473 
409       for(i = 0; i < 200; ++i)                 << 474       G4double random = G4UniformRand() ;
                                                   >> 475       for(i=0;i<200;i++)
410       {                                           476       {
411         if(random >= fIntegralProbabilityOfSR[ << 477         if(random >= fIntegralProbabilityOfSR[i]) break ;
412           break;                               << 
413       }                                           478       }
414       energyOfSR = 0.0001 * i * i * fEnergyCon << 479       energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
415                                                   480 
416       // check against insufficient energy        481       // check against insufficient energy
                                                   >> 482 
417       if(energyOfSR <= 0.0)                       483       if(energyOfSR <= 0.0)
418       {                                           484       {
419         return -1.0;                           << 485         return -1.0 ;
420       }                                           486       }
421     }                                          << 487       //G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
                                                   >> 488       //G4ParticleMomentum 
                                                   >> 489       //particleDirection = aDynamicParticle->GetMomentumDirection();
                                                   >> 490 
                                                   >> 491       // Gamma production cut in this material
                                                   >> 492       //G4double 
                                                   >> 493       //gammaEnergyCut = (G4Gamma::GetCutsInEnergy())[aMaterial->GetIndex()];
                                                   >> 494 
                                                   >> 495       // SR photon has energy more than the current material cut
                                                   >> 496       // M-C of its direction
                                                   >> 497       
                                                   >> 498       //G4double Teta = G4UniformRand()/gamma ;    // Very roughly
                                                   >> 499 
                                                   >> 500       //G4double Phi  = twopi * G4UniformRand() ;
                                                   >> 501     }       
422     else                                          502     else
423     {                                             503     {
424       return -1.0;                             << 504       return -1.0 ;
425     }                                             505     }
426   }                                            << 506   }     
427   return energyOfSR;                           << 507   return energyOfSR ;
428 }                                                 508 }
429                                                   509 
430 //////////////////////////////////////////////    510 /////////////////////////////////////////////////////////////////////////////////
431 G4double G4SynchrotronRadiationInMat::GetRando << 511 //
432                                                << 512 //
                                                   >> 513 
                                                   >> 514 G4double G4SynchrotronRadiationInMat::GetRandomEnergySR(G4double gamma, G4double perpB)
433 {                                                 515 {
434   G4int i;                                     << 516   G4int i, iMax;
435   static constexpr G4int iMax = 200;           << 
436   G4double energySR, random, position;            517   G4double energySR, random, position;
437                                                   518 
                                                   >> 519   iMax = 200;
438   random = G4UniformRand();                       520   random = G4UniformRand();
439                                                   521 
440   for(i = 0; i < iMax; ++i)                    << 522   for( i = 0; i < iMax; i++ )
441   {                                               523   {
442     if(random >= fIntegralProbabilityOfSR[i])  << 524     if( random >= fIntegralProbabilityOfSR[i] ) break;
443       break;                                   << 
444   }                                               525   }
445   if(i <= 0)                                   << 526   if(i <= 0 )        position = G4UniformRand(); // 0.
446     position = G4UniformRand();                << 527   else if( i>= iMax) position = G4double(iMax);
447   else if(i >= iMax)                           << 528   else position = i + G4UniformRand(); // -1
448     position = G4double(iMax);                 << 529   // 
449   else                                         << 530   // it was in initial implementation:
450     position = i + G4UniformRand();            << 531   //       energyOfSR = 0.0001*i*i*fEnergyConst*gamma*gamma*perpB ;
451                                                   532 
452   energySR =                                   << 533   energySR = 0.0001*position*position*fEnergyConst*gamma*gamma*perpB;
453     0.0001 * position * position * fEnergyCons << 
454                                                   534 
455   if(energySR < 0.)                            << 535   if( energySR  < 0. ) energySR = 0.;
456     energySR = 0.;                             << 
457                                                   536 
458   return energySR;                                537   return energySR;
459 }                                                 538 }
460                                                   539 
461 //////////////////////////////////////////////    540 /////////////////////////////////////////////////////////////////////////
462 G4double G4SynchrotronRadiationInMat::GetProbS << 541 //
                                                   >> 542 // return
                                                   >> 543 
                                                   >> 544 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt( G4double t)
463 {                                                 545 {
464   G4double result, hypCos2, hypCos = std::cosh << 546   G4double result, hypCos2, hypCos=std::cosh(t);
465                                                   547 
466   hypCos2 = hypCos * hypCos;                   << 548   hypCos2 = hypCos*hypCos;  
467   result = std::cosh(5. * t / 3.) * std::exp(t << 549   result = std::cosh(5.*t/3.)*std::exp(t-fKsi*hypCos); // fKsi > 0. !
468   result /= hypCos2;                              550   result /= hypCos2;
469   return result;                                  551   return result;
470 }                                                 552 }
471                                                   553 
472 //////////////////////////////////////////////    554 ///////////////////////////////////////////////////////////////////////////
                                                   >> 555 //
473 // return the probability to emit SR photon wi    556 // return the probability to emit SR photon with relative energy
474 // energy/energy_c >= ksi                         557 // energy/energy_c >= ksi
475 // for ksi <= 0. P = 1., however the method wo    558 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
476 G4double G4SynchrotronRadiationInMat::GetIntPr << 559 
                                                   >> 560 G4double G4SynchrotronRadiationInMat::GetIntProbSR( G4double ksi)
477 {                                                 561 {
478   if(ksi <= 0.)                                << 562   if (ksi <= 0.) return 1.0;
479     return 1.0;                                << 563   fKsi = ksi; // should be > 0. !
480   fKsi = ksi;  // should be > 0. !             << 
481   G4int n;                                        564   G4int n;
482   G4double result, a;                             565   G4double result, a;
483                                                   566 
484   a = fAlpha;       // always = 0.             << 567   a = fAlpha;        // always = 0.
485   n = fRootNumber;  // around default = 80     << 568   n = fRootNumber;   // around default = 80
486                                                   569 
487   G4Integrator<G4SynchrotronRadiationInMat,    << 570   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
488                G4double (G4SynchrotronRadiatio << 
489     integral;                                  << 
490                                                   571 
491   result = integral.Laguerre(                  << 572   result = integral.Laguerre(this, 
492     this, &G4SynchrotronRadiationInMat::GetPro << 573            &G4SynchrotronRadiationInMat::GetProbSpectrumSRforInt, a, n);
493                                                   574 
494   result *= 3. / 5. / pi;                      << 575   result *= 3./5./pi;
495                                                   576 
496   return result;                                  577   return result;
497 }                                                 578 }
498                                                   579 
499 //////////////////////////////////////////////    580 /////////////////////////////////////////////////////////////////////////
                                                   >> 581 //
500 // return an auxiliary function for K_5/3 inte    582 // return an auxiliary function for K_5/3 integral representation
501 G4double G4SynchrotronRadiationInMat::GetProbS << 
502 {                                              << 
503   G4double result, hypCos = std::cosh(t);      << 
504                                                   583 
505   result = std::cosh(5. * t / 3.) * std::exp(t << 584 G4double G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy( G4double t)
                                                   >> 585 {
                                                   >> 586   G4double result, hypCos=std::cosh(t);
                                                   >> 587   
                                                   >> 588   result = std::cosh(5.*t/3.)*std::exp(t - fKsi*hypCos); // fKsi > 0. !
506   result /= hypCos;                               589   result /= hypCos;
507   return result;                                  590   return result;
508 }                                                 591 }
509                                                   592 
510 //////////////////////////////////////////////    593 ///////////////////////////////////////////////////////////////////////////
                                                   >> 594 //
511 // return the probability to emit SR photon en    595 // return the probability to emit SR photon energy with relative energy
512 // energy/energy_c >= ksi                         596 // energy/energy_c >= ksi
513 // for ksi <= 0. P = 1., however the method wo    597 // for ksi <= 0. P = 1., however the method works for ksi > 0 only!
514 G4double G4SynchrotronRadiationInMat::GetEnerg << 598 
                                                   >> 599 G4double G4SynchrotronRadiationInMat::GetEnergyProbSR( G4double ksi)
515 {                                                 600 {
516   if(ksi <= 0.)                                << 601   if (ksi <= 0.) return 1.0;
517     return 1.0;                                << 602   fKsi = ksi; // should be > 0. !
518   fKsi = ksi;  // should be > 0. !             << 
519   G4int n;                                        603   G4int n;
520   G4double result, a;                             604   G4double result, a;
521                                                   605 
522   a = fAlpha;       // always = 0.             << 606   a = fAlpha;        // always = 0.
523   n = fRootNumber;  // around default = 80     << 607   n = fRootNumber;   // around default = 80
524                                                   608 
525   G4Integrator<G4SynchrotronRadiationInMat,    << 609   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
526                G4double (G4SynchrotronRadiatio << 
527     integral;                                  << 
528                                                   610 
529   result = integral.Laguerre(                  << 611   result = integral.Laguerre(this, 
530     this, &G4SynchrotronRadiationInMat::GetPro << 612            &G4SynchrotronRadiationInMat::GetProbSpectrumSRforEnergy, a, n);
531                                                   613 
532   result *= 9. * std::sqrt(3.) * ksi / 8. / pi << 614   result *= 9.*std::sqrt(3.)*ksi/8./pi;
533                                                   615 
534   return result;                                  616   return result;
535 }                                                 617 }
536                                                   618 
537 //////////////////////////////////////////////    619 /////////////////////////////////////////////////////////////////////////////
538 G4double G4SynchrotronRadiationInMat::GetInteg << 620 //
539 {                                              << 621 // 
540   G4double result, hypCos = std::cosh(t);      << 
541                                                   622 
542   result =                                     << 623 G4double G4SynchrotronRadiationInMat::GetIntegrandForAngleK( G4double t)
543     std::cosh(fOrderAngleK * t) * std::exp(t - << 624 {
                                                   >> 625   G4double result, hypCos=std::cosh(t);
                                                   >> 626   
                                                   >> 627   result  = std::cosh(fOrderAngleK*t)*std::exp(t - fEta*hypCos); // fEta > 0. !
544   result /= hypCos;                               628   result /= hypCos;
545   return result;                                  629   return result;
546 }                                                 630 }
547                                                   631 
548 //////////////////////////////////////////////    632 //////////////////////////////////////////////////////////////////////////
                                                   >> 633 //
549 // Return K 1/3 or 2/3 for angular distributio    634 // Return K 1/3 or 2/3 for angular distribution
550 G4double G4SynchrotronRadiationInMat::GetAngle << 635 
                                                   >> 636 G4double G4SynchrotronRadiationInMat::GetAngleK( G4double eta)
551 {                                                 637 {
552   fEta = eta;  // should be > 0. !             << 638   fEta = eta; // should be > 0. !
553   G4int n;                                        639   G4int n;
554   G4double result, a;                             640   G4double result, a;
555                                                   641 
556   a = fAlpha;       // always = 0.             << 642   a = fAlpha;        // always = 0.
557   n = fRootNumber;  // around default = 80     << 643   n = fRootNumber;   // around default = 80
558                                                   644 
559   G4Integrator<G4SynchrotronRadiationInMat,    << 645   G4Integrator<G4SynchrotronRadiationInMat, G4double(G4SynchrotronRadiationInMat::*)(G4double)> integral;
560                G4double (G4SynchrotronRadiatio << 
561     integral;                                  << 
562                                                   646 
563   result = integral.Laguerre(                  << 647   result = integral.Laguerre(this, 
564     this, &G4SynchrotronRadiationInMat::GetInt << 648            &G4SynchrotronRadiationInMat::GetIntegrandForAngleK, a, n);
565                                                   649 
566   return result;                                  650   return result;
567 }                                                 651 }
568                                                   652 
569 //////////////////////////////////////////////    653 /////////////////////////////////////////////////////////////////////////
                                                   >> 654 //
570 // Relative angle diff distribution for given     655 // 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                                                   656 
578   fOrderAngleK = 1. / 3.;                      << 657 G4double G4SynchrotronRadiationInMat::GetAngleNumberAtGammaKsi( G4double gpsi)
579   funK         = GetAngleK(fEta);              << 658 {
580   funK2        = funK * funK;                  << 659   G4double result, funK, funK2, gpsi2 = gpsi*gpsi;
581                                                   660 
582   result = gpsi2 * funK2 / (1. + gpsi2);       << 661   fPsiGamma    = gpsi;
                                                   >> 662   fEta         = 0.5*fKsi*(1 + gpsi2)*std::sqrt(1 + gpsi2);
                                                   >> 663  
                                                   >> 664   fOrderAngleK = 1./3.;
                                                   >> 665   funK         = GetAngleK(fEta); 
                                                   >> 666   funK2        = funK*funK;
                                                   >> 667 
                                                   >> 668   result       = gpsi2*funK2/(1 + gpsi2);
                                                   >> 669 
                                                   >> 670   fOrderAngleK = 2./3.;
                                                   >> 671   funK         = GetAngleK(fEta); 
                                                   >> 672   funK2        = funK*funK;
                                                   >> 673 
                                                   >> 674   result      += funK2;
                                                   >> 675   result      *= (1 + gpsi2)*fKsi;
                                                   >> 676      
                                                   >> 677   return result;
                                                   >> 678 }
583                                                   679 
584   fOrderAngleK = 2. / 3.;                      << 
585   funK         = GetAngleK(fEta);              << 
586   funK2        = funK * funK;                  << 
587                                                   680 
588   result += funK2;                             << 681 ///////////////////// end of G4SynchrotronRadiationInMat.cc
589   result *= (1. + gpsi2) * fKsi;               << 
590                                                   682 
591   return result;                               << 
592 }                                              << 
593                                                   683