Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/highenergy/src/G4mplIonisationWithDeltaModel.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/highenergy/src/G4mplIonisationWithDeltaModel.cc (Version 11.3.0) and /processes/electromagnetic/highenergy/src/G4mplIonisationWithDeltaModel.cc (Version 10.3)


  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 // $Id: G4mplIonisationWithDeltaModel.cc 97391 2016-06-02 10:08:45Z gcosmo $
 26 //                                                 27 //
 27 // -------------------------------------------     28 // -------------------------------------------------------------------
 28 //                                                 29 //
 29 // GEANT4 Class header file                        30 // GEANT4 Class header file
 30 //                                                 31 //
 31 //                                                 32 //
 32 // File name:     G4mplIonisationWithDeltaMode     33 // File name:     G4mplIonisationWithDeltaModel
 33 //                                                 34 //
 34 // Author:        Vladimir Ivanchenko              35 // Author:        Vladimir Ivanchenko 
 35 //                                                 36 //
 36 // Creation date: 06.09.2005                       37 // Creation date: 06.09.2005
 37 //                                                 38 //
 38 // Modifications:                                  39 // Modifications:
 39 // 12.08.2007 Changing low energy approximatio     40 // 12.08.2007 Changing low energy approximation and extrapolation. 
 40 //            Small bug fixing and refactoring     41 //            Small bug fixing and refactoring (M. Vladymyrov)
 41 // 13.11.2007 Use low-energy asymptotic from [     42 // 13.11.2007 Use low-energy asymptotic from [3] (V.Ivanchenko) 
 42 //                                                 43 //
 43 //                                                 44 //
 44 // -------------------------------------------     45 // -------------------------------------------------------------------
 45 // References                                      46 // References
 46 // [1] Steven P. Ahlen: Energy loss of relativ     47 // [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles, 
 47 //     S.P. Ahlen, Rev. Mod. Phys 52(1980), p1     48 //     S.P. Ahlen, Rev. Mod. Phys 52(1980), p121
 48 // [2] K.A. Milton arXiv:hep-ex/0602040            49 // [2] K.A. Milton arXiv:hep-ex/0602040
 49 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev.     50 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347
 50                                                    51 
 51                                                    52 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    55 
 55 #include "G4mplIonisationWithDeltaModel.hh"        56 #include "G4mplIonisationWithDeltaModel.hh"
 56 #include "Randomize.hh"                            57 #include "Randomize.hh"
 57 #include "G4PhysicalConstants.hh"                  58 #include "G4PhysicalConstants.hh"
 58 #include "G4SystemOfUnits.hh"                      59 #include "G4SystemOfUnits.hh"
 59 #include "G4ParticleChangeForLoss.hh"              60 #include "G4ParticleChangeForLoss.hh"
 60 #include "G4Electron.hh"                           61 #include "G4Electron.hh"
 61 #include "G4DynamicParticle.hh"                    62 #include "G4DynamicParticle.hh"
 62 #include "G4ProductionCutsTable.hh"                63 #include "G4ProductionCutsTable.hh"
 63 #include "G4MaterialCutsCouple.hh"                 64 #include "G4MaterialCutsCouple.hh"
 64 #include "G4Log.hh"                                65 #include "G4Log.hh"
 65 #include "G4Pow.hh"                            << 
 66                                                    66 
 67 //....oooOO0OOooo........oooOO0OOooo........oo     67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 68                                                    68 
 69 using namespace std;                               69 using namespace std;
 70                                                    70 
 71 std::vector<G4double>* G4mplIonisationWithDelt     71 std::vector<G4double>* G4mplIonisationWithDeltaModel::dedx0 = nullptr;
 72                                                    72 
 73 G4mplIonisationWithDeltaModel::G4mplIonisation     73 G4mplIonisationWithDeltaModel::G4mplIonisationWithDeltaModel(G4double mCharge,
 74                                                <<  74                    const G4String& nam)
 75   : G4VEmModel(nam),G4VEmFluctuationModel(nam)     75   : G4VEmModel(nam),G4VEmFluctuationModel(nam),
 76   magCharge(mCharge),                              76   magCharge(mCharge),
 77   twoln10(std::log(100.0)),                    <<  77   twoln10(log(100.0)),
 78   betalow(0.01),                                   78   betalow(0.01),
 79   betalim(0.1),                                    79   betalim(0.1),
 80   beta2lim(betalim*betalim),                       80   beta2lim(betalim*betalim),
 81   bg2lim(beta2lim*(1.0 + beta2lim))                81   bg2lim(beta2lim*(1.0 + beta2lim))
 82 {                                                  82 {
 83   nmpl = G4lrint(std::abs(magCharge) * 2 * fin <<  83   nmpl = G4lrint(std::fabs(magCharge) * 2 * fine_structure_const);
 84   if(nmpl > 6)      { nmpl = 6; }                  84   if(nmpl > 6)      { nmpl = 6; }
 85   else if(nmpl < 1) { nmpl = 1; }                  85   else if(nmpl < 1) { nmpl = 1; }
 86   pi_hbarc2_over_mc2 = pi * hbarc * hbarc / el     86   pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
 87   chargeSquare = magCharge * magCharge;            87   chargeSquare = magCharge * magCharge;
 88   dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;               88   dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
 89   fParticleChange = nullptr;                       89   fParticleChange = nullptr;
 90   theElectron = G4Electron::Electron();            90   theElectron = G4Electron::Electron();
 91   G4cout << "### Monopole ionisation model wit     91   G4cout << "### Monopole ionisation model with d-electron production, Gmag= " 
 92          << magCharge/eplus << G4endl;         <<  92    << magCharge/eplus << G4endl;
 93   monopole = nullptr;                              93   monopole = nullptr;
 94   mass = 0.0;                                      94   mass = 0.0;
 95 }                                                  95 }
 96                                                    96 
 97 //....oooOO0OOooo........oooOO0OOooo........oo     97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98                                                    98 
 99 G4mplIonisationWithDeltaModel::~G4mplIonisatio     99 G4mplIonisationWithDeltaModel::~G4mplIonisationWithDeltaModel()
100 {                                                 100 {
101   if(IsMaster()) { delete dedx0; }                101   if(IsMaster()) { delete dedx0; }
102 }                                                 102 }
103                                                   103 
104 //....oooOO0OOooo........oooOO0OOooo........oo    104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105                                                   105 
106 void G4mplIonisationWithDeltaModel::SetParticl    106 void G4mplIonisationWithDeltaModel::SetParticle(const G4ParticleDefinition* p)
107 {                                                 107 {
108   monopole = p;                                   108   monopole = p;
109   mass     = monopole->GetPDGMass();              109   mass     = monopole->GetPDGMass();
110   G4double emin =                                 110   G4double emin = 
111     std::min(LowEnergyLimit(),0.1*mass*(1./sqr    111     std::min(LowEnergyLimit(),0.1*mass*(1./sqrt(1. - betalow*betalow) - 1.)); 
112   G4double emax =                                 112   G4double emax = 
113     std::max(HighEnergyLimit(),10*mass*(1./sqr    113     std::max(HighEnergyLimit(),10*mass*(1./sqrt(1. - beta2lim) - 1.)); 
114   SetLowEnergyLimit(emin);                        114   SetLowEnergyLimit(emin);
115   SetHighEnergyLimit(emax);                       115   SetHighEnergyLimit(emax);
116 }                                                 116 }
117                                                   117 
118 //....oooOO0OOooo........oooOO0OOooo........oo    118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119                                                   119 
120 void                                              120 void 
121 G4mplIonisationWithDeltaModel::Initialise(cons    121 G4mplIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
122                                           cons << 122             const G4DataVector&)
123 {                                                 123 {
124   if(!monopole) { SetParticle(p); }               124   if(!monopole) { SetParticle(p); }
125   if(!fParticleChange) { fParticleChange = Get    125   if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
126   if(IsMaster()) {                                126   if(IsMaster()) {
127     if(!dedx0) { dedx0 = new std::vector<G4dou    127     if(!dedx0) { dedx0 = new std::vector<G4double>; }
128     G4ProductionCutsTable* theCoupleTable=        128     G4ProductionCutsTable* theCoupleTable=
129       G4ProductionCutsTable::GetProductionCuts    129       G4ProductionCutsTable::GetProductionCutsTable();
130     G4int numOfCouples = (G4int)theCoupleTable << 130     G4int numOfCouples = theCoupleTable->GetTableSize();
131     G4int n = (G4int)dedx0->size();            << 131     G4int n = dedx0->size();
132     if(n < numOfCouples) { dedx0->resize(numOf    132     if(n < numOfCouples) { dedx0->resize(numOfCouples); }
133     G4Pow* g4calc = G4Pow::GetInstance();      << 
134                                                   133 
135     // initialise vector assuming low conducti << 134     // initialise vector
136     for(G4int i=0; i<numOfCouples; ++i) {         135     for(G4int i=0; i<numOfCouples; ++i) {
137                                                   136 
138       const G4Material* material =                137       const G4Material* material = 
139         theCoupleTable->GetMaterialCutsCouple( << 138   theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
140       G4double eDensity = material->GetElectro    139       G4double eDensity = material->GetElectronDensity();
141       G4double vF2 = 2*electron_Compton_length << 140       G4double vF = electron_Compton_length*pow(3.*pi*pi*eDensity,0.3333333333);
142       (*dedx0)[i] = pi_hbarc2_over_mc2*eDensit    141       (*dedx0)[i] = pi_hbarc2_over_mc2*eDensity*nmpl*nmpl*
143         (G4Log(vF2/fine_structure_const) - 0.5 << 142   (G4Log(2*vF/fine_structure_const) - 0.5)/vF;
144     }                                             143     }
145   }                                               144   }
146 }                                                 145 }
147                                                   146 
148 //....oooOO0OOooo........oooOO0OOooo........oo << 
149                                                << 
150 G4double                                       << 
151 G4mplIonisationWithDeltaModel::MinEnergyCut(co << 
152                                             co << 
153 {                                              << 
154   return couple->GetMaterial()->GetIonisation( << 
155 }                                              << 
156                                                << 
157 //....oooOO0OOooo........oooOO0OOooo........oo    147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158                                                   148 
159 G4double                                          149 G4double 
160 G4mplIonisationWithDeltaModel::ComputeDEDXPerV    150 G4mplIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
161                                                << 151                 const G4ParticleDefinition* p,
162                                                << 152                 G4double kineticEnergy,
163                                                << 153                 G4double maxEnergy)
164 {                                                 154 {
165   if(!monopole) { SetParticle(p); }               155   if(!monopole) { SetParticle(p); }
166   G4double tmax = MaxSecondaryEnergy(p,kinetic    156   G4double tmax = MaxSecondaryEnergy(p,kineticEnergy);
167   G4double cutEnergy = std::min(tmax, maxEnerg    157   G4double cutEnergy = std::min(tmax, maxEnergy);
168   cutEnergy = std::max(LowEnergyLimit(), cutEn    158   cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
169   G4double tau   = kineticEnergy / mass;          159   G4double tau   = kineticEnergy / mass;
170   G4double gam   = tau + 1.0;                     160   G4double gam   = tau + 1.0;
171   G4double bg2   = tau * (tau + 2.0);             161   G4double bg2   = tau * (tau + 2.0);
172   G4double beta2 = bg2 / (gam * gam);             162   G4double beta2 = bg2 / (gam * gam);
173   G4double beta  = sqrt(beta2);                   163   G4double beta  = sqrt(beta2);
174                                                   164 
175   // low-energy asymptotic formula                165   // low-energy asymptotic formula
                                                   >> 166   //G4double dedx  = dedxlim*beta*material->GetDensity();
176   G4double dedx = (*dedx0)[CurrentCouple()->Ge    167   G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()]*beta;
177                                                   168 
178   // above asymptotic                             169   // above asymptotic
179   if(beta > betalow) {                            170   if(beta > betalow) {
180                                                   171 
181     // high energy                                172     // high energy
182     if(beta >= betalim) {                         173     if(beta >= betalim) {
183       dedx = ComputeDEDXAhlen(material, bg2, c    174       dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
184                                                   175 
185     } else {                                      176     } else {
                                                   >> 177 
                                                   >> 178       //G4double dedx1 = dedxlim*betalow*material->GetDensity();
186       G4double dedx1 = (*dedx0)[CurrentCouple(    179       G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()]*betalow;
187       G4double dedx2 = ComputeDEDXAhlen(materi    180       G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
188                                                   181 
189       // extrapolation between two formula        182       // extrapolation between two formula 
190       G4double kapa2 = beta - betalow;            183       G4double kapa2 = beta - betalow;
191       G4double kapa1 = betalim - beta;            184       G4double kapa1 = betalim - beta;
192       dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa    185       dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
193     }                                             186     }
194   }                                               187   }
195   return dedx;                                    188   return dedx;
196 }                                                 189 }
197                                                   190 
198 //....oooOO0OOooo........oooOO0OOooo........oo    191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199                                                   192 
200 G4double                                          193 G4double 
201 G4mplIonisationWithDeltaModel::ComputeDEDXAhle    194 G4mplIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material, 
202                                                << 195             G4double bg2, 
203                                                << 196             G4double cutEnergy)
204 {                                                 197 {
205   G4double eDensity = material->GetElectronDen    198   G4double eDensity = material->GetElectronDensity();
206   G4double eexc  = material->GetIonisation()->    199   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
207                                                   200 
208   // Ahlen's formula for nonconductors, [1]p15    201   // Ahlen's formula for nonconductors, [1]p157, f(5.7)
209   G4double dedx =                                 202   G4double dedx = 
210     0.5*(G4Log(2.0*electron_mass_c2*bg2*cutEne << 203     0.5*(log(2.0 * electron_mass_c2 * bg2*cutEnergy / (eexc*eexc)) - 1.0);
211                                                   204 
212   // Kazama et al. cross-section correction       205   // Kazama et al. cross-section correction
213   G4double  k = 0.406;                            206   G4double  k = 0.406;
214   if(nmpl > 1) { k = 0.346; }                     207   if(nmpl > 1) { k = 0.346; }
215                                                   208 
216   // Bloch correction                             209   // Bloch correction
217   const G4double B[7] = { 0.0, 0.248, 0.672, 1    210   const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685}; 
218                                                   211 
219   dedx += 0.5 * k - B[nmpl];                      212   dedx += 0.5 * k - B[nmpl];
220                                                   213 
221   // density effect correction                    214   // density effect correction
222   G4double x = G4Log(bg2)/twoln10;                215   G4double x = G4Log(bg2)/twoln10;
223   dedx -= material->GetIonisation()->DensityCo    216   dedx -= material->GetIonisation()->DensityCorrection(x);
224                                                   217 
225   // now compute the total ionization loss        218   // now compute the total ionization loss
226   dedx *=  pi_hbarc2_over_mc2 * eDensity * nmp    219   dedx *=  pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
227                                                   220 
228   dedx = std::max(dedx, 0.0);                  << 221   if (dedx < 0.0) { dedx = 0.; }
229   return dedx;                                    222   return dedx;
230 }                                                 223 }
231                                                   224 
232 //....oooOO0OOooo........oooOO0OOooo........oo    225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233                                                   226 
234 G4double                                          227 G4double
235 G4mplIonisationWithDeltaModel::ComputeCrossSec    228 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerElectron(
236                                            con    229                                            const G4ParticleDefinition* p,
237                                            G4d << 230              G4double kineticEnergy,
238                                            G4d << 231              G4double cut,
239                                            G4d << 232              G4double maxKinEnergy)
240 {                                                 233 {
241   if(!monopole) { SetParticle(p); }               234   if(!monopole) { SetParticle(p); }
                                                   >> 235   G4double cross = 0.0;
242   G4double tmax = MaxSecondaryEnergy(p, kineti    236   G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
243   G4double maxEnergy = std::min(tmax, maxKinEn << 237   G4double maxEnergy = std::min(tmax,maxKinEnergy);
244   G4double cutEnergy = std::max(LowEnergyLimit    238   G4double cutEnergy = std::max(LowEnergyLimit(), cut);
245   G4double cross = (cutEnergy < maxEnergy)     << 239   if(cutEnergy < maxEnergy) {
246     ? (0.5/cutEnergy - 0.5/maxEnergy)*pi_hbarc << 240     cross = (0.5/cutEnergy - 0.5/maxEnergy)*pi_hbarc2_over_mc2 * nmpl * nmpl;
                                                   >> 241   }
247   return cross;                                   242   return cross;
248 }                                                 243 }
249                                                   244 
250 //....oooOO0OOooo........oooOO0OOooo........oo    245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
251                                                   246 
252 G4double                                          247 G4double 
253 G4mplIonisationWithDeltaModel::ComputeCrossSec    248 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerAtom(
254                                           cons << 249             const G4ParticleDefinition* p,
255                                           G4do << 250             G4double kineticEnergy,
256                                           G4do << 251             G4double Z, G4double,
257                                           G4do << 252             G4double cutEnergy,
258                                           G4do << 253             G4double maxEnergy)
259 {                                                 254 {
260   G4double cross =                                255   G4double cross = 
261     Z*ComputeCrossSectionPerElectron(p,kinetic    256     Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
262   return cross;                                   257   return cross;
263 }                                                 258 }
264                                                   259 
265 //....oooOO0OOooo........oooOO0OOooo........oo    260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266                                                   261 
267 void                                              262 void 
268 G4mplIonisationWithDeltaModel::SampleSecondari    263 G4mplIonisationWithDeltaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
269                                                << 264              const G4MaterialCutsCouple*,
270                                                << 265              const G4DynamicParticle* dp,
271                                                << 266              G4double minKinEnergy,
272                                                << 267              G4double maxEnergy)
273 {                                                 268 {
274   G4double kineticEnergy = dp->GetKineticEnerg    269   G4double kineticEnergy = dp->GetKineticEnergy();
275   G4double tmax = MaxSecondaryEnergy(dp->GetDe    270   G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
276                                                   271 
277   G4double maxKinEnergy = std::min(maxEnergy,t    272   G4double maxKinEnergy = std::min(maxEnergy,tmax);
278   if(minKinEnergy >= maxKinEnergy) { return; }    273   if(minKinEnergy >= maxKinEnergy) { return; }
279                                                   274 
280   //G4cout << "G4mplIonisationWithDeltaModel::    275   //G4cout << "G4mplIonisationWithDeltaModel::SampleSecondaries: E(GeV)= "
281   //         << kineticEnergy/GeV << " M(GeV)= << 276   //   << kineticEnergy/GeV << " M(GeV)= " << mass/GeV
282   //         << " tmin(MeV)= " << minKinEnergy << 277   //   << " tmin(MeV)= " << minKinEnergy/MeV << G4endl;
283                                                   278 
284   G4double totEnergy     = kineticEnergy + mas    279   G4double totEnergy     = kineticEnergy + mass;
285   G4double etot2         = totEnergy*totEnergy    280   G4double etot2         = totEnergy*totEnergy;
286   G4double beta2         = kineticEnergy*(kine    281   G4double beta2         = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
287                                                   282   
288   // sampling without nuclear size effect         283   // sampling without nuclear size effect
289   G4double q = G4UniformRand();                   284   G4double q = G4UniformRand();
290   G4double deltaKinEnergy = minKinEnergy*maxKi    285   G4double deltaKinEnergy = minKinEnergy*maxKinEnergy
291     /(minKinEnergy*(1.0 - q) + maxKinEnergy*q)    286     /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
292                                                   287 
293   // delta-electron is produced                   288   // delta-electron is produced
294   G4double totMomentum = totEnergy*sqrt(beta2)    289   G4double totMomentum = totEnergy*sqrt(beta2);
295   G4double deltaMomentum =                        290   G4double deltaMomentum =
296            sqrt(deltaKinEnergy * (deltaKinEner    291            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
297   G4double cost = deltaKinEnergy * (totEnergy     292   G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
298                                    (deltaMomen    293                                    (deltaMomentum * totMomentum);
299   cost = std::min(cost, 1.0);                  << 294   if(cost > 1.0) { cost = 1.0; }
300                                                   295 
301   G4double sint = sqrt((1.0 - cost)*(1.0 + cos    296   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
302                                                   297 
303   G4double phi = twopi * G4UniformRand() ;        298   G4double phi = twopi * G4UniformRand() ;
304                                                   299 
305   G4ThreeVector deltaDirection(sint*cos(phi),s    300   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
306   G4ThreeVector direction = dp->GetMomentumDir    301   G4ThreeVector direction = dp->GetMomentumDirection();
307   deltaDirection.rotateUz(direction);             302   deltaDirection.rotateUz(direction);
308                                                   303 
309   // create G4DynamicParticle object for delta    304   // create G4DynamicParticle object for delta ray
310   G4DynamicParticle* delta =                      305   G4DynamicParticle* delta = 
311     new G4DynamicParticle(theElectron,deltaDir    306     new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
312                                                   307 
313   vdp->push_back(delta);                          308   vdp->push_back(delta);
314                                                   309 
315   // Change kinematics of primary particle        310   // Change kinematics of primary particle
316   kineticEnergy       -= deltaKinEnergy;          311   kineticEnergy       -= deltaKinEnergy;
317   G4ThreeVector finalP = direction*totMomentum    312   G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
318   finalP               = finalP.unit();           313   finalP               = finalP.unit();
319                                                   314 
320   fParticleChange->SetProposedKineticEnergy(ki    315   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
321   fParticleChange->SetProposedMomentumDirectio    316   fParticleChange->SetProposedMomentumDirection(finalP);
322 }                                                 317 }
323                                                   318 
324 //....oooOO0OOooo........oooOO0OOooo........oo    319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325                                                   320 
326 G4double G4mplIonisationWithDeltaModel::Sample    321 G4double G4mplIonisationWithDeltaModel::SampleFluctuations(
327                                        const G << 322                const G4MaterialCutsCouple* couple,
328                                        const G << 323                const G4DynamicParticle* dp,
329                                        const G << 324                G4double tmax,
330                                        const G << 325                G4double length,
331                                        const G << 326                G4double meanLoss)
332                                        const G << 
333 {                                                 327 {
334   G4double siga = Dispersion(couple->GetMateri << 328   G4double siga = Dispersion(couple->GetMaterial(),dp,tmax,length);
335   G4double loss = meanLoss;                       329   G4double loss = meanLoss;
336   siga = std::sqrt(siga);                      << 330   siga = sqrt(siga);
337   G4double twomeanLoss = meanLoss + meanLoss;     331   G4double twomeanLoss = meanLoss + meanLoss;
338                                                   332 
339   if(twomeanLoss < siga) {                        333   if(twomeanLoss < siga) {
340     G4double x;                                   334     G4double x;
341     do {                                          335     do {
342       loss = twomeanLoss*G4UniformRand();         336       loss = twomeanLoss*G4UniformRand();
343       x = (loss - meanLoss)/siga;                 337       x = (loss - meanLoss)/siga;
344       // Loop checking, 07-Aug-2015, Vladimir     338       // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
345     } while (1.0 - 0.5*x*x < G4UniformRand());    339     } while (1.0 - 0.5*x*x < G4UniformRand());
346   } else {                                        340   } else {
347     do {                                          341     do {
348       loss = G4RandGauss::shoot(meanLoss,siga)    342       loss = G4RandGauss::shoot(meanLoss,siga);
349       // Loop checking, 07-Aug-2015, Vladimir     343       // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
350     } while (0.0 > loss || loss > twomeanLoss)    344     } while (0.0 > loss || loss > twomeanLoss);
351   }                                               345   }
352   return loss;                                    346   return loss;
353 }                                                 347 }
354                                                   348 
355 //....oooOO0OOooo........oooOO0OOooo........oo    349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
356                                                   350 
357 G4double                                          351 G4double 
358 G4mplIonisationWithDeltaModel::Dispersion(cons    352 G4mplIonisationWithDeltaModel::Dispersion(const G4Material* material,
359                                           cons << 353             const G4DynamicParticle* dp,
360             const G4double tcut,               << 354             G4double tmax,
361             const G4double tmax,               << 355             G4double length)
362             const G4double length)             << 
363 {                                                 356 {
364   G4double siga = 0.0;                            357   G4double siga = 0.0;
365   G4double tau   = dp->GetKineticEnergy()/mass    358   G4double tau   = dp->GetKineticEnergy()/mass;
366   if(tau > 0.0) {                                 359   if(tau > 0.0) { 
367     const G4double beta = dp->GetBeta();       << 360     G4double electronDensity = material->GetElectronDensity();
368     siga  = (tmax/(beta*beta) - 0.5*tcut) * tw << 361     G4double gam   = tau + 1.0;
369       * material->GetElectronDensity() * charg << 362     G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
                                                   >> 363     siga  = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
                                                   >> 364       * electronDensity * chargeSquare;
370   }                                               365   }
371   return siga;                                    366   return siga;
372 }                                                 367 }
373                                                   368 
374 //....oooOO0OOooo........oooOO0OOooo........oo    369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
375                                                   370 
376 G4double                                          371 G4double 
377 G4mplIonisationWithDeltaModel::MaxSecondaryEne    372 G4mplIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
378                                                << 373               G4double kinEnergy)
379 {                                                 374 {
380   G4double tau = kinEnergy/mass;                  375   G4double tau = kinEnergy/mass;
381   return 2.0*electron_mass_c2*tau*(tau + 2.);     376   return 2.0*electron_mass_c2*tau*(tau + 2.);
382 }                                                 377 }
383                                                   378 
384 //....oooOO0OOooo........oooOO0OOooo........oo    379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
385                                                   380