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 10.7.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 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 
                                                   >>  55 //
                                                   >>  56 // Class Description:
                                                   >>  57 //
                                                   >>  58 //
 54 // -------------------------------------------     59 // -------------------------------------------------------------------
 55 //                                                 60 //
 56 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 57 //....oooOO0OOooo........oooOO0OOooo........oo     62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 58                                                    63 
 59 #include "G4MuBremsstrahlungModel.hh"              64 #include "G4MuBremsstrahlungModel.hh"
 60 #include "G4PhysicalConstants.hh"                  65 #include "G4PhysicalConstants.hh"
 61 #include "G4SystemOfUnits.hh"                      66 #include "G4SystemOfUnits.hh"
 62 #include "G4Gamma.hh"                              67 #include "G4Gamma.hh"
 63 #include "G4MuonMinus.hh"                          68 #include "G4MuonMinus.hh"
 64 #include "G4MuonPlus.hh"                           69 #include "G4MuonPlus.hh"
 65 #include "Randomize.hh"                            70 #include "Randomize.hh"
 66 #include "G4Material.hh"                           71 #include "G4Material.hh"
 67 #include "G4Element.hh"                            72 #include "G4Element.hh"
 68 #include "G4ElementVector.hh"                      73 #include "G4ElementVector.hh"
 69 #include "G4ProductionCutsTable.hh"                74 #include "G4ProductionCutsTable.hh"
 70 #include "G4ModifiedMephi.hh"                      75 #include "G4ModifiedMephi.hh"
 71 #include "G4ParticleChangeForLoss.hh"              76 #include "G4ParticleChangeForLoss.hh"
 72 #include "G4Log.hh"                                77 #include "G4Log.hh"
                                                   >>  78 #include "G4Exp.hh"
 73                                                    79 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 75 //....oooOO0OOooo........oooOO0OOooo........oo     81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 76                                                    82 
                                                   >>  83 using namespace std;
                                                   >>  84 
 77 const G4double G4MuBremsstrahlungModel::xgi[]      85 const G4double G4MuBremsstrahlungModel::xgi[] = 
 78   {0.03377,0.16940,0.38069,0.61931,0.83060,0.9     86   {0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
 79 const G4double G4MuBremsstrahlungModel::wgi[]      87 const G4double G4MuBremsstrahlungModel::wgi[] = 
 80   {0.08566,0.18038,0.23396,0.23396,0.18038,0.0     88   {0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
 81 G4double G4MuBremsstrahlungModel::fDN[] = {0.0     89 G4double G4MuBremsstrahlungModel::fDN[] = {0.0};
 82                                                    90 
 83 G4MuBremsstrahlungModel::G4MuBremsstrahlungMod     91 G4MuBremsstrahlungModel::G4MuBremsstrahlungModel(const G4ParticleDefinition* p,
 84                                                    92                                                  const G4String& nam)
 85   : G4VEmModel(nam),                               93   : G4VEmModel(nam),
 86     sqrte(std::sqrt(G4Exp(1.))),               <<  94     particle(nullptr),
 87     lowestKinEnergy(0.1*CLHEP::GeV),           <<  95     sqrte(sqrt(G4Exp(1.))),
 88     minThreshold(0.9*CLHEP::keV)               <<  96     bh(202.4),
                                                   >>  97     bh1(446.),
                                                   >>  98     btf(183.),
                                                   >>  99     btf1(1429.),
                                                   >> 100     fParticleChange(nullptr),
                                                   >> 101     lowestKinEnergy(1.0*GeV),
                                                   >> 102     minThreshold(0.9*keV)
 89 {                                                 103 {
 90   theGamma = G4Gamma::Gamma();                    104   theGamma = G4Gamma::Gamma();
 91   nist = G4NistManager::Instance();            << 105   nist = G4NistManager::Instance();
 92                                                   106 
 93   SetAngularDistribution(new G4ModifiedMephi() << 107   lowestKinEnergy = 1.*GeV;  
 94                                                   108 
 95   if (nullptr != p) { SetParticle(p); }        << 109   mass = rmass = cc = coeff = 1.0;
 96   if (0.0 == fDN[1]) {                         << 110 
 97     for (G4int i=1; i<93; ++i) {               << 111   if(0.0 == fDN[1]) {
                                                   >> 112     for(G4int i=1; i<93; ++i) {
 98       G4double dn = 1.54*nist->GetA27(i);         113       G4double dn = 1.54*nist->GetA27(i);
 99       fDN[i] = dn;                                114       fDN[i] = dn;
100       if(1 < i) {                                 115       if(1 < i) {
101   fDN[i] /= std::pow(dn, 1./G4double(i));      << 116         fDN[i] /= std::pow(dn, 1./G4double(i));
102       }                                           117       }
103     }                                             118     }
104   }                                               119   }
                                                   >> 120   SetAngularDistribution(new G4ModifiedMephi());
                                                   >> 121 
                                                   >> 122   if(p) { SetParticle(p); }
105 }                                                 123 }
106                                                   124 
107 //....oooOO0OOooo........oooOO0OOooo........oo    125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108                                                   126 
109 G4double G4MuBremsstrahlungModel::MinEnergyCut    127 G4double G4MuBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
110                                                   128                                                const G4MaterialCutsCouple*)
111 {                                                 129 {
112   return minThreshold;                            130   return minThreshold;
113 }                                                 131 }
114                                                   132 
115 //....oooOO0OOooo........oooOO0OOooo........oo    133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116                                                   134 
117 G4double G4MuBremsstrahlungModel::MinPrimaryEn    135 G4double G4MuBremsstrahlungModel::MinPrimaryEnergy(const G4Material*,
118                                                   136                                                    const G4ParticleDefinition*,
119                                                   137                                                    G4double cut)
120 {                                                 138 {
121   return std::max(lowestKinEnergy, cut);       << 139   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 }                                                 140 }
136                                                   141 
137 //....oooOO0OOooo........oooOO0OOooo........oo    142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138                                                   143 
139 void G4MuBremsstrahlungModel::Initialise(const    144 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
140                                          const    145                                          const G4DataVector& cuts)
141 {                                                 146 {
142   SetParticle(p);                              << 147   if(p) { SetParticle(p); }
143   if(nullptr == fParticleChange) {             << 148 
144     fParticleChange = GetParticleChangeForLoss << 149   // define pointer to G4ParticleChange
145   }                                            << 150   if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
                                                   >> 151 
146   if(IsMaster() && p == particle && lowestKinE    152   if(IsMaster() && p == particle && lowestKinEnergy < HighEnergyLimit()) { 
147     InitialiseElementSelectors(p, cuts);          153     InitialiseElementSelectors(p, cuts); 
148   }                                               154   }
149 }                                                 155 }
150                                                   156 
151 //....oooOO0OOooo........oooOO0OOooo........oo    157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152                                                   158 
153 void G4MuBremsstrahlungModel::InitialiseLocal(    159 void G4MuBremsstrahlungModel::InitialiseLocal(const G4ParticleDefinition* p,
154                                                   160                                               G4VEmModel* masterModel)
155 {                                                 161 {
156   if(p == particle && lowestKinEnergy < HighEn    162   if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
157     SetElementSelectors(masterModel->GetElemen    163     SetElementSelectors(masterModel->GetElementSelectors());
158   }                                               164   }
159 }                                                 165 }
160                                                   166 
161 //....oooOO0OOooo........oooOO0OOooo........oo    167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162                                                   168 
163 G4double G4MuBremsstrahlungModel::ComputeDEDXP    169 G4double G4MuBremsstrahlungModel::ComputeDEDXPerVolume(
164                                                   170                                               const G4Material* material,
165                                                   171                                               const G4ParticleDefinition*,
166                                                   172                                                     G4double kineticEnergy,
167                                                   173                                                     G4double cutEnergy)
168 {                                                 174 {
169   G4double dedx = 0.0;                            175   G4double dedx = 0.0;
170   if (kineticEnergy <= lowestKinEnergy) { retu    176   if (kineticEnergy <= lowestKinEnergy) { return dedx; }
171                                                   177 
172   G4double cut = std::max(cutEnergy, minThresh << 178   G4double tmax = kineticEnergy;
173   cut = std::min(cut, kineticEnergy);          << 179   G4double cut  = std::min(cutEnergy,tmax);
                                                   >> 180   if(cut < minThreshold) { cut = minThreshold; }
174                                                   181 
175   const G4ElementVector* theElementVector = ma    182   const G4ElementVector* theElementVector = material->GetElementVector();
176   const G4double* theAtomicNumDensityVector =     183   const G4double* theAtomicNumDensityVector =
177     material->GetAtomicNumDensityVector();        184     material->GetAtomicNumDensityVector();
178                                                   185 
179   //  loop for elements in the material           186   //  loop for elements in the material
180   for (size_t i=0; i<material->GetNumberOfElem << 187   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
                                                   >> 188 
181     G4double loss =                               189     G4double loss = 
182       ComputMuBremLoss((*theElementVector)[i]-    190       ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
                                                   >> 191 
183     dedx += loss*theAtomicNumDensityVector[i];    192     dedx += loss*theAtomicNumDensityVector[i];
184   }                                               193   }
185   //  G4cout << "BR e= " << kineticEnergy << "    194   //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
186   dedx = std::max(dedx, 0.);                   << 195   if(dedx < 0.) dedx = 0.;
187   return dedx;                                    196   return dedx;
188 }                                                 197 }
189                                                   198 
190 //....oooOO0OOooo........oooOO0OOooo........oo    199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191                                                   200 
192 G4double G4MuBremsstrahlungModel::ComputMuBrem    201 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
193                                                   202                                                    G4double tkin, G4double cut)
194 {                                                 203 {
195   G4double totalEnergy = mass + tkin;             204   G4double totalEnergy = mass + tkin;
196   static const G4double ak1 = 0.05;               205   static const G4double ak1 = 0.05;
197   static const G4int k2 = 5;                   << 206   static const G4int    k2=5;
198   G4double loss = 0.;                             207   G4double loss = 0.;
199                                                   208 
200   G4double vcut = cut/totalEnergy;                209   G4double vcut = cut/totalEnergy;
201   G4int kkk = (G4int)(vcut/ak1) + k2;          << 210   G4double vmax = tkin/totalEnergy;
202   if (kkk > 8) { kkk = 8; }                    << 211 
203   else if (kkk < 1) { kkk = 1; }               << 212   G4double aaa = 0.;
204   G4double hhh = vcut/(G4double)(kkk);         << 213   G4double bbb = vcut;
205                                                << 214   if(vcut>vmax) { bbb = vmax; }
206   G4double aa = 0.;                            << 215   G4int kkk = (G4int)((bbb-aaa)/ak1)+k2;
207   for(G4int l=0; l<kkk; ++l) {                 << 216   if(kkk < 1) { kkk = 1; }
208     for(G4int i=0; i<6; ++i) {                 << 217 
                                                   >> 218   G4double hhh=(bbb-aaa)/G4double(kkk);
                                                   >> 219 
                                                   >> 220   G4double aa = aaa;
                                                   >> 221   for(G4int l=0; l<kkk; l++)
                                                   >> 222   {
                                                   >> 223     for(G4int i=0; i<6; i++)
                                                   >> 224     {
209       G4double ep = (aa + xgi[i]*hhh)*totalEne    225       G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
210       loss += ep*wgi[i]*ComputeDMicroscopicCro    226       loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
211     }                                             227     }
212     aa += hhh;                                    228     aa += hhh;
213   }                                               229   }
214                                                   230 
215   loss *= hhh*totalEnergy;                     << 231   loss *=hhh*totalEnergy ;
                                                   >> 232 
216   return loss;                                    233   return loss;
217 }                                                 234 }
218                                                   235 
219 //....oooOO0OOooo........oooOO0OOooo........oo    236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
220                                                   237 
221 G4double G4MuBremsstrahlungModel::ComputeMicro    238 G4double G4MuBremsstrahlungModel::ComputeMicroscopicCrossSection(
222                                            G4d    239                                            G4double tkin,
223                                            G4d    240                                            G4double Z,
224                                            G4d    241                                            G4double cut)
225 {                                                 242 {
226   G4double totalEnergy = tkin + mass;             243   G4double totalEnergy = tkin + mass;
227   static const G4double ak1 = 2.3;                244   static const G4double ak1 = 2.3;
228   static const G4int k2 = 4;                   << 245   static const G4int    k2  = 4;
229   G4double cross = 0.;                            246   G4double cross = 0.;
230                                                   247 
231   if(cut >= tkin) return cross;                   248   if(cut >= tkin) return cross;
232                                                   249 
233   G4double vcut = cut/totalEnergy;                250   G4double vcut = cut/totalEnergy;
234   G4double vmax = tkin/totalEnergy;               251   G4double vmax = tkin/totalEnergy;
235                                                   252 
236   G4double aaa = G4Log(vcut);                     253   G4double aaa = G4Log(vcut);
237   G4double bbb = G4Log(vmax);                     254   G4double bbb = G4Log(vmax);
238   G4int kkk = (G4int)((bbb-aaa)/ak1) + k2 ;    << 255   G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
239   if(kkk > 8) { kkk = 8; }                     << 256   if(kkk < 1) { kkk = 1; }
240   else if (kkk < 1) { kkk = 1; }               << 257 
241   G4double hhh = (bbb-aaa)/(G4double)(kkk);    << 258   G4double hhh = (bbb-aaa)/G4double(kkk);
                                                   >> 259 
242   G4double aa = aaa;                              260   G4double aa = aaa;
243                                                   261 
244   for(G4int l=0; l<kkk; ++l) {                 << 262   for(G4int l=0; l<kkk; l++)
245     for(G4int i=0; i<6; ++i) {                 << 263   {
                                                   >> 264     for(G4int i=0; i<6; i++)
                                                   >> 265     {
246       G4double ep = G4Exp(aa + xgi[i]*hhh)*tot    266       G4double ep = G4Exp(aa + xgi[i]*hhh)*totalEnergy;
247       cross += ep*wgi[i]*ComputeDMicroscopicCr    267       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
248     }                                             268     }
249     aa += hhh;                                    269     aa += hhh;
250   }                                               270   }
251                                                   271 
252   cross *= hhh;                                << 272   cross *=hhh;
                                                   >> 273 
253   //G4cout << "BR e= " << tkin<< "  cross= " <    274   //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
                                                   >> 275 
254   return cross;                                   276   return cross;
255 }                                                 277 }
256                                                   278 
257 //....oooOO0OOooo........oooOO0OOooo........oo    279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
258                                                   280 
259 G4double G4MuBremsstrahlungModel::ComputeDMicr    281 G4double G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection(
260                                            G4d    282                                            G4double tkin,
261                                            G4d    283                                            G4double Z,
262                                            G4d    284                                            G4double gammaEnergy)
263 //  differential cross section                    285 //  differential cross section
264 {                                                 286 {
265   G4double dxsection = 0.;                        287   G4double dxsection = 0.;
                                                   >> 288 
266   if(gammaEnergy > tkin) { return dxsection; }    289   if(gammaEnergy > tkin) { return dxsection; }
267                                                   290 
268   G4double E = tkin + mass ;                      291   G4double E = tkin + mass ;
269   G4double v = gammaEnergy/E ;                    292   G4double v = gammaEnergy/E ;
270   G4double delta = 0.5*mass*mass*v/(E-gammaEne    293   G4double delta = 0.5*mass*mass*v/(E-gammaEnergy) ;
271   G4double rab0  = delta*sqrte ;                  294   G4double rab0  = delta*sqrte ;
272                                                   295 
273   G4int iz = G4lrint(Z);                          296   G4int iz = G4lrint(Z);
274   if(iz < 1) { iz = 1; }                          297   if(iz < 1) { iz = 1; }
275   else if(iz > 92) { iz = 92; }                   298   else if(iz > 92) { iz = 92; }
276                                                   299 
277   G4double z13 = 1.0/nist->GetZ13(iz);            300   G4double z13 = 1.0/nist->GetZ13(iz);
278   G4double dnstar = fDN[iz];                      301   G4double dnstar = fDN[iz];
279                                                   302 
280   G4double b,b1;                                  303   G4double b,b1;
                                                   >> 304 
281   if(1 == iz) {                                   305   if(1 == iz) {
282     b  = bh;                                      306     b  = bh;
283     b1 = bh1;                                     307     b1 = bh1;
284   } else {                                        308   } else {
285     b  = btf;                                     309     b  = btf;
286     b1 = btf1;                                    310     b1 = btf1;
287   }                                               311   }
288                                                   312 
289   // nucleus contribution logarithm               313   // nucleus contribution logarithm
290   G4double rab1 = b*z13;                       << 314   G4double rab1=b*z13;
291   G4double fn = G4Log(rab1/(dnstar*(CLHEP::ele << 315   G4double fn=G4Log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
292   (mass + delta*(dnstar*sqrte-2.)));           << 316               (mass+delta*(dnstar*sqrte-2.))) ;
293   fn = std::max(fn, 0.);                       << 317   if(fn <0.) { fn = 0.; }
294   // electron contribution logarithm              318   // electron contribution logarithm
295   G4double epmax1 = E/(1.+0.5*mass*rmass/E);   << 319   G4double epmax1=E/(1.+0.5*mass*rmass/E) ;
296   G4double fe = 0.;                            << 320   G4double fe=0.;
297   if(gammaEnergy < epmax1) {                   << 321   if(gammaEnergy<epmax1)
298     G4double rab2 = b1*z13*z13;                << 322   {
299     fe = G4Log(rab2*mass/((1.+delta*rmass/(CLH << 323     G4double rab2=b1*z13*z13 ;
300   (CLHEP::electron_mass_c2+rab0*rab2)));       << 324     fe=G4Log(rab2*mass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
301     fe = std::max(fe, 0.);                     << 325                               (electron_mass_c2+rab0*rab2))) ;
                                                   >> 326     if(fe<0.) { fe=0.; }
302   }                                               327   }
303                                                   328 
304   dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn    329   dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn*Z + fe)/gammaEnergy;
305   dxsection = std::max(dxsection, 0.0);        << 330 
306   return dxsection;                               331   return dxsection;
307 }                                                 332 }
308                                                   333 
309 //....oooOO0OOooo........oooOO0OOooo........oo    334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310                                                   335 
311 G4double G4MuBremsstrahlungModel::ComputeCross    336 G4double G4MuBremsstrahlungModel::ComputeCrossSectionPerAtom(
312                                            con    337                                            const G4ParticleDefinition*,
313                                                   338                                                  G4double kineticEnergy,
314                                                   339                                                  G4double Z, G4double,
315                                                   340                                                  G4double cutEnergy,
316                                                   341                                                  G4double maxEnergy)
317 {                                                 342 {
318   G4double cross = 0.0;                           343   G4double cross = 0.0;
319   if (kineticEnergy <= lowestKinEnergy) return    344   if (kineticEnergy <= lowestKinEnergy) return cross;
320   G4double tmax = std::min(maxEnergy, kineticE    345   G4double tmax = std::min(maxEnergy, kineticEnergy);
321   G4double cut  = std::min(cutEnergy, kineticE    346   G4double cut  = std::min(cutEnergy, kineticEnergy);
322   if (cut < minThreshold) cut = minThreshold;  << 347   if(cut < minThreshold) cut = minThreshold;
323   if (cut >= tmax) return cross;                  348   if (cut >= tmax) return cross;
324                                                   349 
325   cross = ComputeMicroscopicCrossSection (kine    350   cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
326   if(tmax < kineticEnergy) {                      351   if(tmax < kineticEnergy) {
327     cross -= ComputeMicroscopicCrossSection(ki    352     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
328   }                                               353   }
329   return cross;                                   354   return cross;
330 }                                                 355 }
331                                                   356 
332 //....oooOO0OOooo........oooOO0OOooo........oo    357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333                                                   358 
334 void G4MuBremsstrahlungModel::SampleSecondarie    359 void G4MuBremsstrahlungModel::SampleSecondaries(
335                               std::vector<G4Dy    360                               std::vector<G4DynamicParticle*>* vdp,
336                               const G4Material    361                               const G4MaterialCutsCouple* couple,
337                               const G4DynamicP    362                               const G4DynamicParticle* dp,
338                               G4double minEner    363                               G4double minEnergy,
339                               G4double maxEner    364                               G4double maxEnergy)
340 {                                                 365 {
341   G4double kineticEnergy = dp->GetKineticEnerg    366   G4double kineticEnergy = dp->GetKineticEnergy();
342   // check against insufficient energy            367   // check against insufficient energy
343   G4double tmax = std::min(kineticEnergy, maxE    368   G4double tmax = std::min(kineticEnergy, maxEnergy);
344   G4double tmin = std::min(kineticEnergy, minE    369   G4double tmin = std::min(kineticEnergy, minEnergy);
345   tmin = std::max(tmin, minThreshold);         << 370   if(tmin < minThreshold) tmin = minThreshold;
346   if(tmin >= tmax) return;                        371   if(tmin >= tmax) return;
347                                                   372 
348   // ===== sampling of energy transfer ======     373   // ===== sampling of energy transfer ======
349                                                   374 
350   G4ParticleMomentum partDirection = dp->GetMo    375   G4ParticleMomentum partDirection = dp->GetMomentumDirection();
351                                                   376 
352   // select randomly one element constituing t    377   // select randomly one element constituing the material
353   const G4Element* anElement = SelectRandomAto    378   const G4Element* anElement = SelectRandomAtom(couple,particle,kineticEnergy);
354   G4double Z = anElement->GetZ();                 379   G4double Z = anElement->GetZ();
                                                   >> 380 
355   G4double func1 = tmin*                          381   G4double func1 = tmin*
356     ComputeDMicroscopicCrossSection(kineticEne << 382     ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
357                                                   383 
358   G4double gEnergy;                            << 384   G4double lnepksi, epksi;
359   G4double func2;                                 385   G4double func2;
360                                                   386 
361   G4double xmin = G4Log(tmin/minThreshold);    << 387   G4double xmin = G4Log(tmin/CLHEP::MeV);
362   G4double xmax = G4Log(tmax/tmin);            << 388   G4double xmax = G4Log(kineticEnergy/tmin);
363                                                   389 
364   do {                                            390   do {
365     gEnergy = minThreshold*G4Exp(xmin + G4Unif << 391     lnepksi = xmin + G4UniformRand()*xmax;
366     func2 = gEnergy*ComputeDMicroscopicCrossSe << 392     epksi   = CLHEP::MeV*G4Exp(lnepksi);
367                                                << 393     func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
                                                   >> 394 
368     // Loop checking, 03-Aug-2015, Vladimir Iv    395     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
369   } while(func2 < func1*G4UniformRand());         396   } while(func2 < func1*G4UniformRand());
370                                                   397 
                                                   >> 398   G4double gEnergy = epksi;
                                                   >> 399 
                                                   >> 400   // ===== sample angle =====
                                                   >> 401 
                                                   >> 402   //
371   // angles of the emitted gamma using general    403   // angles of the emitted gamma using general interface
                                                   >> 404 
372   G4ThreeVector gamDir =                          405   G4ThreeVector gamDir = 
373     GetAngularDistribution()->SampleDirection(    406     GetAngularDistribution()->SampleDirection(dp, gEnergy, Z, 
374                                                   407                                               couple->GetMaterial());
375   // create G4DynamicParticle object for the G    408   // create G4DynamicParticle object for the Gamma
376   G4DynamicParticle* gamma = new G4DynamicPart    409   G4DynamicParticle* gamma = new G4DynamicParticle(theGamma, gamDir, gEnergy);
377   vdp->push_back(gamma);                          410   vdp->push_back(gamma);
378                                                << 
379   // compute post-interaction kinematics of pr    411   // compute post-interaction kinematics of primary e-/e+ based on
380   // energy-momentum conservation                 412   // energy-momentum conservation
381   const G4double totMomentum =                    413   const G4double totMomentum = 
382     std::sqrt(kineticEnergy*(kineticEnergy + 2    414     std::sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
383   G4ThreeVector dir =                             415   G4ThreeVector dir =
384     (totMomentum*dp->GetMomentumDirection() -  << 416     (totMomentum*dp->GetMomentumDirection()-gEnergy*gamDir).unit();
385   const G4double finalE = kineticEnergy - gEne    417   const G4double finalE = kineticEnergy - gEnergy;
386                                                << 
387   // if secondary gamma energy is higher than     418   // if secondary gamma energy is higher than threshold(very high by default)
388   // then stop tracking the primary particle a    419   // then stop tracking the primary particle and create new secondary e-/e+
389   // instead of the primary one                   420   // instead of the primary one
390   if (gEnergy > SecondaryThreshold()) {           421   if (gEnergy > SecondaryThreshold()) {
391     fParticleChange->ProposeTrackStatus(fStopA    422     fParticleChange->ProposeTrackStatus(fStopAndKill);
392     fParticleChange->SetProposedKineticEnergy(    423     fParticleChange->SetProposedKineticEnergy(0.0);
393     G4DynamicParticle* newdp = new G4DynamicPa    424     G4DynamicParticle* newdp = new G4DynamicParticle(particle, dir, finalE);
394     vdp->push_back(newdp);                        425     vdp->push_back(newdp);
395   } else {                                     << 426   } else { // continue tracking the primary e-/e+ otherwise
396     // continue tracking the primary e-/e+ oth << 
397     fParticleChange->SetProposedMomentumDirect    427     fParticleChange->SetProposedMomentumDirection(dir);
398     fParticleChange->SetProposedKineticEnergy(    428     fParticleChange->SetProposedKineticEnergy(finalE);
399   }                                               429   }
400 }                                                 430 }
401                                                   431 
402 //....oooOO0OOooo........oooOO0OOooo........oo    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
403                                                   433