Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.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/muons/src/G4MuBremsstrahlungModel.cc (Version 11.3.0) and /processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc (Version 11.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 //                                                 26 //
 27 // -------------------------------------------     27 // -------------------------------------------------------------------
 28 //                                                 28 //
 29 // GEANT4 Class file                               29 // GEANT4 Class file
 30 //                                                 30 //
 31 //                                                 31 //
 32 // File name:     G4MuBremsstrahlungModel          32 // File name:     G4MuBremsstrahlungModel
 33 //                                                 33 //
 34 // Author:        Vladimir Ivanchenko on base      34 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //                                                 35 //
 36 // Creation date: 24.06.2002                       36 // Creation date: 24.06.2002
 37 //                                                 37 //
 38 // Modifications:                                  38 // Modifications:
 39 //                                                 39 //
 40 // 04-12-02 Change G4DynamicParticle construct     40 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
 41 // 23-12-02 Change interface in order to move      41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 42 // 24-01-03 Fix for compounds (V.Ivanchenko)       42 // 24-01-03 Fix for compounds (V.Ivanchenko)
 43 // 27-01-03 Make models region aware (V.Ivanch     43 // 27-01-03 Make models region aware (V.Ivanchenko)
 44 // 13-02-03 Add name (V.Ivanchenko)                44 // 13-02-03 Add name (V.Ivanchenko)
 45 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)     45 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
 46 // 08-04-05 Major optimisation of internal int     46 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
 47 // 03-08-05 Angular correlations according to      47 // 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
 48 // 13-02-06 add ComputeCrossSectionPerAtom (mm     48 // 13-02-06 add ComputeCrossSectionPerAtom (mma)
 49 // 21-03-06 Fix problem of initialisation in c     49 // 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
 50 // 07-11-07 Improve sampling of final state (A     50 // 07-11-07 Improve sampling of final state (A.Bogdanov)
 51 // 28-02-08 Use precomputed Z^1/3 and Log(A) (     51 // 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
 52 // 31-05-13 Use element selectors instead of l     52 // 31-05-13 Use element selectors instead of local data structure (V.Ivanchenko)
 53 //                                                 53 //
 54 // -------------------------------------------     54 // -------------------------------------------------------------------
 55 //                                                 55 //
 56 //....oooOO0OOooo........oooOO0OOooo........oo     56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 57 //....oooOO0OOooo........oooOO0OOooo........oo     57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    58 
 59 #include "G4MuBremsstrahlungModel.hh"              59 #include "G4MuBremsstrahlungModel.hh"
 60 #include "G4PhysicalConstants.hh"                  60 #include "G4PhysicalConstants.hh"
 61 #include "G4SystemOfUnits.hh"                      61 #include "G4SystemOfUnits.hh"
 62 #include "G4Gamma.hh"                              62 #include "G4Gamma.hh"
 63 #include "G4MuonMinus.hh"                          63 #include "G4MuonMinus.hh"
 64 #include "G4MuonPlus.hh"                           64 #include "G4MuonPlus.hh"
 65 #include "Randomize.hh"                            65 #include "Randomize.hh"
 66 #include "G4Material.hh"                           66 #include "G4Material.hh"
 67 #include "G4Element.hh"                            67 #include "G4Element.hh"
 68 #include "G4ElementVector.hh"                      68 #include "G4ElementVector.hh"
 69 #include "G4ProductionCutsTable.hh"                69 #include "G4ProductionCutsTable.hh"
 70 #include "G4ModifiedMephi.hh"                      70 #include "G4ModifiedMephi.hh"
 71 #include "G4ParticleChangeForLoss.hh"              71 #include "G4ParticleChangeForLoss.hh"
 72 #include "G4Log.hh"                                72 #include "G4Log.hh"
                                                   >>  73 #include "G4Exp.hh"
                                                   >>  74 #include "G4NistManager.hh"
 73                                                    75 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 75 //....oooOO0OOooo........oooOO0OOooo........oo     77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 76                                                    78 
                                                   >>  79 using namespace std;
                                                   >>  80 
 77 const G4double G4MuBremsstrahlungModel::xgi[]      81 const G4double G4MuBremsstrahlungModel::xgi[] = 
 78   {0.03377,0.16940,0.38069,0.61931,0.83060,0.9     82   {0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
 79 const G4double G4MuBremsstrahlungModel::wgi[]      83 const G4double G4MuBremsstrahlungModel::wgi[] = 
 80   {0.08566,0.18038,0.23396,0.23396,0.18038,0.0     84   {0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
 81 G4double G4MuBremsstrahlungModel::fDN[] = {0.0     85 G4double G4MuBremsstrahlungModel::fDN[] = {0.0};
                                                   >>  86 #ifdef G4MULTITHREADED
                                                   >>  87   G4Mutex G4MuBremsstrahlungModel::theMuBremMutex = G4MUTEX_INITIALIZER;
                                                   >>  88 #endif
 82                                                    89 
 83 G4MuBremsstrahlungModel::G4MuBremsstrahlungMod     90 G4MuBremsstrahlungModel::G4MuBremsstrahlungModel(const G4ParticleDefinition* p,
 84                                                    91                                                  const G4String& nam)
 85   : G4VEmModel(nam),                               92   : G4VEmModel(nam),
 86     sqrte(std::sqrt(G4Exp(1.))),               <<  93     sqrte(sqrt(G4Exp(1.))),
 87     lowestKinEnergy(0.1*CLHEP::GeV),           <<  94     lowestKinEnergy(1.0*CLHEP::GeV),
 88     minThreshold(0.9*CLHEP::keV)                   95     minThreshold(0.9*CLHEP::keV)
 89 {                                                  96 {
 90   theGamma = G4Gamma::Gamma();                     97   theGamma = G4Gamma::Gamma();
 91   nist = G4NistManager::Instance();                98   nist = G4NistManager::Instance();  
 92                                                    99 
 93   SetAngularDistribution(new G4ModifiedMephi() << 100   if(0.0 == fDN[1]) {
 94                                                << 101 #ifdef G4MULTITHREADED
 95   if (nullptr != p) { SetParticle(p); }        << 102     G4MUTEXLOCK(&theMuBremMutex);
 96   if (0.0 == fDN[1]) {                         << 103     if(0.0 == fDN[1]) {
 97     for (G4int i=1; i<93; ++i) {               << 104 #endif
 98       G4double dn = 1.54*nist->GetA27(i);      << 105       for(G4int i=1; i<93; ++i) {
 99       fDN[i] = dn;                             << 106         G4double dn = 1.54*nist->GetA27(i);
100       if(1 < i) {                              << 107         fDN[i] = dn;
101   fDN[i] /= std::pow(dn, 1./G4double(i));      << 108         if(1 < i) {
                                                   >> 109           fDN[i] /= std::pow(dn, 1./G4double(i));
                                                   >> 110         }
102       }                                           111       }
                                                   >> 112 #ifdef G4MULTITHREADED
103     }                                             113     }
                                                   >> 114     G4MUTEXUNLOCK(&theMuBremMutex);
                                                   >> 115 #endif
104   }                                               116   }
                                                   >> 117   SetAngularDistribution(new G4ModifiedMephi());
                                                   >> 118 
                                                   >> 119   if(nullptr != p) { SetParticle(p); }
105 }                                                 120 }
106                                                   121 
107 //....oooOO0OOooo........oooOO0OOooo........oo    122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108                                                   123 
109 G4double G4MuBremsstrahlungModel::MinEnergyCut    124 G4double G4MuBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
110                                                   125                                                const G4MaterialCutsCouple*)
111 {                                                 126 {
112   return minThreshold;                            127   return minThreshold;
113 }                                                 128 }
114                                                   129 
115 //....oooOO0OOooo........oooOO0OOooo........oo    130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116                                                   131 
117 G4double G4MuBremsstrahlungModel::MinPrimaryEn    132 G4double G4MuBremsstrahlungModel::MinPrimaryEnergy(const G4Material*,
118                                                   133                                                    const G4ParticleDefinition*,
119                                                   134                                                    G4double cut)
120 {                                                 135 {
121   return std::max(lowestKinEnergy, cut);       << 136   return std::max(lowestKinEnergy,cut);
122 }                                              << 
123                                                << 
124 //....oooOO0OOooo........oooOO0OOooo........oo << 
125                                                << 
126 void G4MuBremsstrahlungModel::SetParticle(cons << 
127 {                                              << 
128   if(nullptr == particle) {                    << 
129     particle = p;                              << 
130     mass = particle->GetPDGMass();             << 
131     rmass = mass/CLHEP::electron_mass_c2 ;     << 
132     cc = CLHEP::classic_electr_radius/rmass ;  << 
133     coeff = 16.*CLHEP::fine_structure_const*cc << 
134   }                                            << 
135 }                                                 137 }
136                                                   138 
137 //....oooOO0OOooo........oooOO0OOooo........oo    139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138                                                   140 
139 void G4MuBremsstrahlungModel::Initialise(const    141 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
140                                          const    142                                          const G4DataVector& cuts)
141 {                                                 143 {
142   SetParticle(p);                                 144   SetParticle(p);
                                                   >> 145 
143   if(nullptr == fParticleChange) {                146   if(nullptr == fParticleChange) { 
144     fParticleChange = GetParticleChangeForLoss    147     fParticleChange = GetParticleChangeForLoss(); 
145   }                                               148   }
                                                   >> 149 
146   if(IsMaster() && p == particle && lowestKinE    150   if(IsMaster() && p == particle && lowestKinEnergy < HighEnergyLimit()) { 
147     InitialiseElementSelectors(p, cuts);          151     InitialiseElementSelectors(p, cuts); 
148   }                                               152   }
149 }                                                 153 }
150                                                   154 
151 //....oooOO0OOooo........oooOO0OOooo........oo    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152                                                   156 
153 void G4MuBremsstrahlungModel::InitialiseLocal(    157 void G4MuBremsstrahlungModel::InitialiseLocal(const G4ParticleDefinition* p,
154                                                   158                                               G4VEmModel* masterModel)
155 {                                                 159 {
156   if(p == particle && lowestKinEnergy < HighEn    160   if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
157     SetElementSelectors(masterModel->GetElemen    161     SetElementSelectors(masterModel->GetElementSelectors());
158   }                                               162   }
159 }                                                 163 }
160                                                   164 
161 //....oooOO0OOooo........oooOO0OOooo........oo    165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162                                                   166 
163 G4double G4MuBremsstrahlungModel::ComputeDEDXP    167 G4double G4MuBremsstrahlungModel::ComputeDEDXPerVolume(
164                                                   168                                               const G4Material* material,
165                                                   169                                               const G4ParticleDefinition*,
166                                                   170                                                     G4double kineticEnergy,
167                                                   171                                                     G4double cutEnergy)
168 {                                                 172 {
169   G4double dedx = 0.0;                            173   G4double dedx = 0.0;
170   if (kineticEnergy <= lowestKinEnergy) { retu    174   if (kineticEnergy <= lowestKinEnergy) { return dedx; }
171                                                   175 
172   G4double cut = std::max(cutEnergy, minThresh << 176   G4double tmax = kineticEnergy;
173   cut = std::min(cut, kineticEnergy);          << 177   G4double cut  = std::min(cutEnergy,tmax);
                                                   >> 178   if(cut < minThreshold) { cut = minThreshold; }
174                                                   179 
175   const G4ElementVector* theElementVector = ma    180   const G4ElementVector* theElementVector = material->GetElementVector();
176   const G4double* theAtomicNumDensityVector =     181   const G4double* theAtomicNumDensityVector =
177     material->GetAtomicNumDensityVector();        182     material->GetAtomicNumDensityVector();
178                                                   183 
179   //  loop for elements in the material           184   //  loop for elements in the material
180   for (size_t i=0; i<material->GetNumberOfElem    185   for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
                                                   >> 186 
181     G4double loss =                               187     G4double loss = 
182       ComputMuBremLoss((*theElementVector)[i]-    188       ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
                                                   >> 189 
183     dedx += loss*theAtomicNumDensityVector[i];    190     dedx += loss*theAtomicNumDensityVector[i];
184   }                                               191   }
185   //  G4cout << "BR e= " << kineticEnergy << "    192   //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
186   dedx = std::max(dedx, 0.);                   << 193   if(dedx < 0.) dedx = 0.;
187   return dedx;                                    194   return dedx;
188 }                                                 195 }
189                                                   196 
190 //....oooOO0OOooo........oooOO0OOooo........oo    197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191                                                   198 
192 G4double G4MuBremsstrahlungModel::ComputMuBrem    199 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
193                                                   200                                                    G4double tkin, G4double cut)
194 {                                                 201 {
195   G4double totalEnergy = mass + tkin;             202   G4double totalEnergy = mass + tkin;
196   static const G4double ak1 = 0.05;               203   static const G4double ak1 = 0.05;
197   static const G4int k2 = 5;                   << 204   static const G4int    k2=5;
198   G4double loss = 0.;                             205   G4double loss = 0.;
199                                                   206 
200   G4double vcut = cut/totalEnergy;                207   G4double vcut = cut/totalEnergy;
201   G4int kkk = (G4int)(vcut/ak1) + k2;          << 208   G4double vmax = tkin/totalEnergy;
202   if (kkk > 8) { kkk = 8; }                    << 209 
203   else if (kkk < 1) { kkk = 1; }               << 210   G4double aaa = 0.;
204   G4double hhh = vcut/(G4double)(kkk);         << 211   G4double bbb = vcut;
205                                                << 212   if(vcut>vmax) { bbb = vmax; }
206   G4double aa = 0.;                            << 213   G4int kkk = (G4int)((bbb-aaa)/ak1)+k2;
207   for(G4int l=0; l<kkk; ++l) {                 << 214   if(kkk < 1) { kkk = 1; }
208     for(G4int i=0; i<6; ++i) {                 << 215 
                                                   >> 216   G4double hhh=(bbb-aaa)/G4double(kkk);
                                                   >> 217 
                                                   >> 218   G4double aa = aaa;
                                                   >> 219   for(G4int l=0; l<kkk; l++)
                                                   >> 220   {
                                                   >> 221     for(G4int i=0; i<6; i++)
                                                   >> 222     {
209       G4double ep = (aa + xgi[i]*hhh)*totalEne    223       G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
210       loss += ep*wgi[i]*ComputeDMicroscopicCro    224       loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
211     }                                             225     }
212     aa += hhh;                                    226     aa += hhh;
213   }                                               227   }
214                                                   228 
215   loss *= hhh*totalEnergy;                     << 229   loss *=hhh*totalEnergy ;
                                                   >> 230 
216   return loss;                                    231   return loss;
217 }                                                 232 }
218                                                   233 
219 //....oooOO0OOooo........oooOO0OOooo........oo    234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
220                                                   235 
221 G4double G4MuBremsstrahlungModel::ComputeMicro    236 G4double G4MuBremsstrahlungModel::ComputeMicroscopicCrossSection(
222                                            G4d    237                                            G4double tkin,
223                                            G4d    238                                            G4double Z,
224                                            G4d    239                                            G4double cut)
225 {                                                 240 {
226   G4double totalEnergy = tkin + mass;             241   G4double totalEnergy = tkin + mass;
227   static const G4double ak1 = 2.3;                242   static const G4double ak1 = 2.3;
228   static const G4int k2 = 4;                   << 243   static const G4int    k2  = 4;
229   G4double cross = 0.;                            244   G4double cross = 0.;
230                                                   245 
231   if(cut >= tkin) return cross;                   246   if(cut >= tkin) return cross;
232                                                   247 
233   G4double vcut = cut/totalEnergy;                248   G4double vcut = cut/totalEnergy;
234   G4double vmax = tkin/totalEnergy;               249   G4double vmax = tkin/totalEnergy;
235                                                   250 
236   G4double aaa = G4Log(vcut);                     251   G4double aaa = G4Log(vcut);
237   G4double bbb = G4Log(vmax);                     252   G4double bbb = G4Log(vmax);
238   G4int kkk = (G4int)((bbb-aaa)/ak1) + k2 ;    << 253   G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
239   if(kkk > 8) { kkk = 8; }                     << 254   if(kkk < 1) { kkk = 1; }
240   else if (kkk < 1) { kkk = 1; }               << 255 
241   G4double hhh = (bbb-aaa)/(G4double)(kkk);    << 256   G4double hhh = (bbb-aaa)/G4double(kkk);
                                                   >> 257 
242   G4double aa = aaa;                              258   G4double aa = aaa;
243                                                   259 
244   for(G4int l=0; l<kkk; ++l) {                 << 260   for(G4int l=0; l<kkk; ++l)
245     for(G4int i=0; i<6; ++i) {                 << 261   {
                                                   >> 262     for(G4int i=0; i<6; ++i)
                                                   >> 263     {
246       G4double ep = G4Exp(aa + xgi[i]*hhh)*tot    264       G4double ep = G4Exp(aa + xgi[i]*hhh)*totalEnergy;
247       cross += ep*wgi[i]*ComputeDMicroscopicCr    265       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
248     }                                             266     }
249     aa += hhh;                                    267     aa += hhh;
250   }                                               268   }
251                                                   269 
252   cross *= hhh;                                << 270   cross *=hhh;
                                                   >> 271 
253   //G4cout << "BR e= " << tkin<< "  cross= " <    272   //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
                                                   >> 273 
254   return cross;                                   274   return cross;
255 }                                                 275 }
256                                                   276 
257 //....oooOO0OOooo........oooOO0OOooo........oo    277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
258                                                   278 
259 G4double G4MuBremsstrahlungModel::ComputeDMicr    279 G4double G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection(
260                                            G4d    280                                            G4double tkin,
261                                            G4d    281                                            G4double Z,
262                                            G4d    282                                            G4double gammaEnergy)
263 //  differential cross section                    283 //  differential cross section
264 {                                                 284 {
265   G4double dxsection = 0.;                        285   G4double dxsection = 0.;
                                                   >> 286 
266   if(gammaEnergy > tkin) { return dxsection; }    287   if(gammaEnergy > tkin) { return dxsection; }
267                                                   288 
268   G4double E = tkin + mass ;                      289   G4double E = tkin + mass ;
269   G4double v = gammaEnergy/E ;                    290   G4double v = gammaEnergy/E ;
270   G4double delta = 0.5*mass*mass*v/(E-gammaEne    291   G4double delta = 0.5*mass*mass*v/(E-gammaEnergy) ;
271   G4double rab0  = delta*sqrte ;                  292   G4double rab0  = delta*sqrte ;
272                                                   293 
273   G4int iz = G4lrint(Z);                          294   G4int iz = G4lrint(Z);
274   if(iz < 1) { iz = 1; }                          295   if(iz < 1) { iz = 1; }
275   else if(iz > 92) { iz = 92; }                   296   else if(iz > 92) { iz = 92; }
276                                                   297 
277   G4double z13 = 1.0/nist->GetZ13(iz);            298   G4double z13 = 1.0/nist->GetZ13(iz);
278   G4double dnstar = fDN[iz];                      299   G4double dnstar = fDN[iz];
279                                                   300 
280   G4double b,b1;                                  301   G4double b,b1;
                                                   >> 302 
281   if(1 == iz) {                                   303   if(1 == iz) {
282     b  = bh;                                      304     b  = bh;
283     b1 = bh1;                                     305     b1 = bh1;
284   } else {                                        306   } else {
285     b  = btf;                                     307     b  = btf;
286     b1 = btf1;                                    308     b1 = btf1;
287   }                                               309   }
288                                                   310 
289   // nucleus contribution logarithm               311   // nucleus contribution logarithm
290   G4double rab1 = b*z13;                       << 312   G4double rab1=b*z13;
291   G4double fn = G4Log(rab1/(dnstar*(CLHEP::ele << 313   G4double fn=G4Log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
292   (mass + delta*(dnstar*sqrte-2.)));           << 314               (mass+delta*(dnstar*sqrte-2.))) ;
293   fn = std::max(fn, 0.);                       << 315   if(fn <0.) { fn = 0.; }
294   // electron contribution logarithm              316   // electron contribution logarithm
295   G4double epmax1 = E/(1.+0.5*mass*rmass/E);   << 317   G4double epmax1=E/(1.+0.5*mass*rmass/E) ;
296   G4double fe = 0.;                            << 318   G4double fe=0.;
297   if(gammaEnergy < epmax1) {                   << 319   if(gammaEnergy<epmax1)
298     G4double rab2 = b1*z13*z13;                << 320   {
299     fe = G4Log(rab2*mass/((1.+delta*rmass/(CLH << 321     G4double rab2=b1*z13*z13 ;
300   (CLHEP::electron_mass_c2+rab0*rab2)));       << 322     fe=G4Log(rab2*mass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
301     fe = std::max(fe, 0.);                     << 323                               (electron_mass_c2+rab0*rab2))) ;
                                                   >> 324     if(fe<0.) { fe=0.; }
302   }                                               325   }
303                                                   326 
304   dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn    327   dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn*Z + fe)/gammaEnergy;
305   dxsection = std::max(dxsection, 0.0);        << 328 
306   return dxsection;                               329   return dxsection;
307 }                                                 330 }
308                                                   331 
309 //....oooOO0OOooo........oooOO0OOooo........oo    332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310                                                   333 
311 G4double G4MuBremsstrahlungModel::ComputeCross    334 G4double G4MuBremsstrahlungModel::ComputeCrossSectionPerAtom(
312                                            con    335                                            const G4ParticleDefinition*,
313                                                   336                                                  G4double kineticEnergy,
314                                                   337                                                  G4double Z, G4double,
315                                                   338                                                  G4double cutEnergy,
316                                                   339                                                  G4double maxEnergy)
317 {                                                 340 {
318   G4double cross = 0.0;                           341   G4double cross = 0.0;
319   if (kineticEnergy <= lowestKinEnergy) return    342   if (kineticEnergy <= lowestKinEnergy) return cross;
320   G4double tmax = std::min(maxEnergy, kineticE    343   G4double tmax = std::min(maxEnergy, kineticEnergy);
321   G4double cut  = std::min(cutEnergy, kineticE    344   G4double cut  = std::min(cutEnergy, kineticEnergy);
322   if (cut < minThreshold) cut = minThreshold;  << 345   if(cut < minThreshold) cut = minThreshold;
323   if (cut >= tmax) return cross;                  346   if (cut >= tmax) return cross;
324                                                   347 
325   cross = ComputeMicroscopicCrossSection (kine    348   cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
326   if(tmax < kineticEnergy) {                      349   if(tmax < kineticEnergy) {
327     cross -= ComputeMicroscopicCrossSection(ki    350     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
328   }                                               351   }
329   return cross;                                   352   return cross;
330 }                                                 353 }
331                                                   354 
332 //....oooOO0OOooo........oooOO0OOooo........oo    355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333                                                   356 
334 void G4MuBremsstrahlungModel::SampleSecondarie    357 void G4MuBremsstrahlungModel::SampleSecondaries(
335                               std::vector<G4Dy    358                               std::vector<G4DynamicParticle*>* vdp,
336                               const G4Material    359                               const G4MaterialCutsCouple* couple,
337                               const G4DynamicP    360                               const G4DynamicParticle* dp,
338                               G4double minEner    361                               G4double minEnergy,
339                               G4double maxEner    362                               G4double maxEnergy)
340 {                                                 363 {
341   G4double kineticEnergy = dp->GetKineticEnerg    364   G4double kineticEnergy = dp->GetKineticEnergy();
342   // check against insufficient energy            365   // check against insufficient energy
343   G4double tmax = std::min(kineticEnergy, maxE    366   G4double tmax = std::min(kineticEnergy, maxEnergy);
344   G4double tmin = std::min(kineticEnergy, minE    367   G4double tmin = std::min(kineticEnergy, minEnergy);
345   tmin = std::max(tmin, minThreshold);         << 368   if(tmin < minThreshold) tmin = minThreshold;
346   if(tmin >= tmax) return;                        369   if(tmin >= tmax) return;
347                                                   370 
348   // ===== sampling of energy transfer ======     371   // ===== sampling of energy transfer ======
349                                                   372 
350   G4ParticleMomentum partDirection = dp->GetMo    373   G4ParticleMomentum partDirection = dp->GetMomentumDirection();
351                                                   374 
352   // select randomly one element constituing t    375   // select randomly one element constituing the material
353   const G4Element* anElement = SelectRandomAto    376   const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
354   G4double Z = anElement->GetZ();                 377   G4double Z = anElement->GetZ();
                                                   >> 378 
355   G4double func1 = tmin*                          379   G4double func1 = tmin*
356     ComputeDMicroscopicCrossSection(kineticEne << 380     ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
357                                                   381 
358   G4double gEnergy;                            << 382   G4double lnepksi, epksi;
359   G4double func2;                                 383   G4double func2;
360                                                   384 
361   G4double xmin = G4Log(tmin/minThreshold);    << 385   G4double xmin = G4Log(tmin/CLHEP::MeV);
362   G4double xmax = G4Log(tmax/tmin);            << 386   G4double xmax = G4Log(kineticEnergy/tmin);
363                                                   387 
364   do {                                            388   do {
365     gEnergy = minThreshold*G4Exp(xmin + G4Unif << 389     lnepksi = xmin + G4UniformRand()*xmax;
366     func2 = gEnergy*ComputeDMicroscopicCrossSe << 390     epksi   = CLHEP::MeV*G4Exp(lnepksi);
367                                                << 391     func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
                                                   >> 392 
368     // Loop checking, 03-Aug-2015, Vladimir Iv    393     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
369   } while(func2 < func1*G4UniformRand());         394   } while(func2 < func1*G4UniformRand());
370                                                   395 
                                                   >> 396   G4double gEnergy = epksi;
                                                   >> 397 
                                                   >> 398   // ===== sample angle =====
                                                   >> 399 
                                                   >> 400   //
371   // angles of the emitted gamma using general    401   // angles of the emitted gamma using general interface
                                                   >> 402 
372   G4ThreeVector gamDir =                          403   G4ThreeVector gamDir = 
373     GetAngularDistribution()->SampleDirection(    404     GetAngularDistribution()->SampleDirection(dp, gEnergy, Z, 
374                                                   405                                               couple->GetMaterial());
375   // create G4DynamicParticle object for the G    406   // create G4DynamicParticle object for the Gamma
376   G4DynamicParticle* gamma = new G4DynamicPart    407   G4DynamicParticle* gamma = new G4DynamicParticle(theGamma, gamDir, gEnergy);
377   vdp->push_back(gamma);                          408   vdp->push_back(gamma);
378                                                << 
379   // compute post-interaction kinematics of pr    409   // compute post-interaction kinematics of primary e-/e+ based on
380   // energy-momentum conservation                 410   // energy-momentum conservation
381   const G4double totMomentum =                    411   const G4double totMomentum = 
382     std::sqrt(kineticEnergy*(kineticEnergy + 2    412     std::sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
383   G4ThreeVector dir =                             413   G4ThreeVector dir =
384     (totMomentum*dp->GetMomentumDirection() -  << 414     (totMomentum*dp->GetMomentumDirection()-gEnergy*gamDir).unit();
385   const G4double finalE = kineticEnergy - gEne    415   const G4double finalE = kineticEnergy - gEnergy;
386                                                << 
387   // if secondary gamma energy is higher than     416   // if secondary gamma energy is higher than threshold(very high by default)
388   // then stop tracking the primary particle a    417   // then stop tracking the primary particle and create new secondary e-/e+
389   // instead of the primary one                   418   // instead of the primary one
390   if (gEnergy > SecondaryThreshold()) {           419   if (gEnergy > SecondaryThreshold()) {
391     fParticleChange->ProposeTrackStatus(fStopA    420     fParticleChange->ProposeTrackStatus(fStopAndKill);
392     fParticleChange->SetProposedKineticEnergy(    421     fParticleChange->SetProposedKineticEnergy(0.0);
393     G4DynamicParticle* newdp = new G4DynamicPa    422     G4DynamicParticle* newdp = new G4DynamicParticle(particle, dir, finalE);
394     vdp->push_back(newdp);                        423     vdp->push_back(newdp);
395   } else {                                     << 424   } else { // continue tracking the primary e-/e+ otherwise
396     // continue tracking the primary e-/e+ oth << 
397     fParticleChange->SetProposedMomentumDirect    425     fParticleChange->SetProposedMomentumDirection(dir);
398     fParticleChange->SetProposedKineticEnergy(    426     fParticleChange->SetProposedKineticEnergy(finalE);
399   }                                               427   }
400 }                                                 428 }
401                                                   429 
402 //....oooOO0OOooo........oooOO0OOooo........oo    430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
403                                                   431