Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4BraggModel.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/standard/src/G4BraggModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4BraggModel.cc (Version 9.2.p4)


  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: G4BraggModel.cc,v 1.20 2008/10/22 16:01:46 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // -------------------------------------------------------------------
 28 //                                                 30 //
 29 // GEANT4 Class file                               31 // GEANT4 Class file
 30 //                                                 32 //
 31 //                                                 33 //
 32 // File name:   G4BraggModel                       34 // File name:   G4BraggModel
 33 //                                                 35 //
 34 // Author:        Vladimir Ivanchenko              36 // Author:        Vladimir Ivanchenko
 35 //                                                 37 //
 36 // Creation date: 03.01.2002                       38 // Creation date: 03.01.2002
 37 //                                                 39 //
 38 // Modifications:                                  40 // Modifications: 
 39 //                                                 41 //
 40 // 04-12-02 Fix problem of G4DynamicParticle c     42 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
 41 // 23-12-02 Change interface in order to move      43 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 42 // 27-01-03 Make models region aware (V.Ivanch     44 // 27-01-03 Make models region aware (V.Ivanchenko)
 43 // 13-02-03 Add name (V.Ivanchenko)                45 // 13-02-03 Add name (V.Ivanchenko)
 44 // 04-06-03 Fix compilation warnings (V.Ivanch     46 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
 45 // 12-09-04 Add lowestKinEnergy and change ord     47 // 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI)
 46 // 11-04-05 Major optimisation of internal int     48 // 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
 47 // 16-06-05 Fix problem of chemical formula (V     49 // 16-06-05 Fix problem of chemical formula (V.Ivantchenko)
 48 // 15-02-06 ComputeCrossSectionPerElectron, Co     50 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 49 // 25-04-06 Add stopping data from PSTAR (V.Iv     51 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
 50 // 12-08-08 Added methods GetParticleCharge, G     52 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 
 51 //          CorrectionsAlongStep needed for io     53 //          CorrectionsAlongStep needed for ions(V.Ivanchenko)
 52                                                    54 
 53 // Class Description:                              55 // Class Description:
 54 //                                                 56 //
 55 // Implementation of energy loss and delta-ele     57 // Implementation of energy loss and delta-electron production by
 56 // slow charged heavy particles                    58 // slow charged heavy particles
 57                                                    59 
 58 // -------------------------------------------     60 // -------------------------------------------------------------------
 59 //                                                 61 //
 60                                                    62 
 61                                                    63 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 63 //....oooOO0OOooo........oooOO0OOooo........oo     65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64                                                    66 
 65 #include "G4BraggModel.hh"                         67 #include "G4BraggModel.hh"
 66 #include "G4PhysicalConstants.hh"              << 
 67 #include "G4SystemOfUnits.hh"                  << 
 68 #include "Randomize.hh"                            68 #include "Randomize.hh"
 69 #include "G4Electron.hh"                           69 #include "G4Electron.hh"
 70 #include "G4ParticleChangeForLoss.hh"              70 #include "G4ParticleChangeForLoss.hh"
 71 #include "G4LossTableManager.hh"                   71 #include "G4LossTableManager.hh"
 72 #include "G4EmCorrections.hh"                      72 #include "G4EmCorrections.hh"
 73 #include "G4EmParameters.hh"                   << 
 74 #include "G4DeltaAngle.hh"                     << 
 75 #include "G4ICRU90StoppingData.hh"             << 
 76 #include "G4NistManager.hh"                    << 
 77 #include "G4Log.hh"                            << 
 78 #include "G4Exp.hh"                            << 
 79 #include "G4AutoLock.hh"                       << 
 80                                                    73 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 82                                                    75 
 83 G4ICRU90StoppingData* G4BraggModel::fICRU90 =  <<  76 using namespace std;
 84 G4PSTARStopping* G4BraggModel::fPSTAR = nullpt << 
 85                                                << 
 86 namespace                                      << 
 87 {                                              << 
 88   G4Mutex ionMutex = G4MUTEX_INITIALIZER;      << 
 89 }                                              << 
 90                                                    77 
 91 G4BraggModel::G4BraggModel(const G4ParticleDef     78 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
 92   : G4VEmModel(nam)                            <<  79   : G4VEmModel(nam),
                                                   >>  80     particle(0),
                                                   >>  81     protonMassAMU(1.007276),
                                                   >>  82     iMolecula(0),
                                                   >>  83     isIon(false),
                                                   >>  84     isInitialised(false)
 93 {                                                  85 {
 94   SetHighEnergyLimit(2.0*CLHEP::MeV);          <<  86   if(p) SetParticle(p);
                                                   >>  87   SetHighEnergyLimit(2.0*MeV);
 95                                                    88 
 96   lowestKinEnergy  = 0.25*CLHEP::keV;          <<  89   lowestKinEnergy  = 1.0*keV;
 97   theZieglerFactor = CLHEP::eV*CLHEP::cm2*1.0e <<  90   theZieglerFactor = eV*cm2*1.0e-15;
 98   theElectron = G4Electron::Electron();            91   theElectron = G4Electron::Electron();
 99   expStopPower125 = 0.0;                       << 
100                                                << 
101   corr = G4LossTableManager::Instance()->EmCor << 
102   if(nullptr != p) { SetParticle(p); }         << 
103   else { SetParticle(theElectron); }           << 
104 }                                                  92 }
105                                                    93 
106 //....oooOO0OOooo........oooOO0OOooo........oo     94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107                                                    95 
108 G4BraggModel::~G4BraggModel()                      96 G4BraggModel::~G4BraggModel()
                                                   >>  97 {}
                                                   >>  98 
                                                   >>  99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 100 
                                                   >> 101 G4double G4BraggModel::MinEnergyCut(const G4ParticleDefinition*,
                                                   >> 102                                     const G4MaterialCutsCouple* couple)
109 {                                                 103 {
110   if(isFirst) {                                << 104   return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
111     delete fPSTAR;                             << 
112     fPSTAR = nullptr;                          << 
113   }                                            << 
114 }                                                 105 }
115                                                   106 
116 //....oooOO0OOooo........oooOO0OOooo........oo    107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117                                                   108 
118 void G4BraggModel::Initialise(const G4Particle    109 void G4BraggModel::Initialise(const G4ParticleDefinition* p,
119                               const G4DataVect    110                               const G4DataVector&)
120 {                                                 111 {
121   if(p != particle) { SetParticle(p); }        << 112   if(p != particle) SetParticle(p);
122                                                << 
123   // always false before the run               << 
124   SetDeexcitationFlag(false);                  << 
125                                                   113 
126   // initialise data only once                 << 114   if(!isInitialised) {
127   if(nullptr == fPSTAR) {                      << 115     isInitialised = true;
128     G4AutoLock l(&ionMutex);                   << 
129     if(nullptr == fPSTAR) {                    << 
130       isFirst = true;                          << 
131       fPSTAR = new G4PSTARStopping();          << 
132       if(G4EmParameters::Instance()->UseICRU90 << 
133   fICRU90 = G4NistManager::Instance()->GetICRU << 
134       }                                        << 
135     }                                          << 
136     l.unlock();                                << 
137   }                                            << 
138   if(isFirst) {                                << 
139     if(nullptr != fICRU90) { fICRU90->Initiali << 
140     fPSTAR->Initialise();                      << 
141   }                                            << 
142                                                << 
143   if(nullptr == fParticleChange) {             << 
144                                                   116 
145     if(UseAngularGeneratorFlag() && !GetAngula << 
146       SetAngularDistribution(new G4DeltaAngle( << 
147     }                                          << 
148     G4String pname = particle->GetParticleName    117     G4String pname = particle->GetParticleName();
149     if(particle->GetParticleType() == "nucleus    118     if(particle->GetParticleType() == "nucleus" && 
150        pname != "deuteron" && pname != "triton << 119        pname != "deuteron" && pname != "triton") isIon = true;
151        pname != "alpha+"   && pname != "helium << 
152        pname != "hydrogen") { isIon = true; }  << 
153                                                   120 
154     fParticleChange = GetParticleChangeForLoss << 121     corr = G4LossTableManager::Instance()->EmCorrections();
155   }                                            << 
156 }                                              << 
157                                                << 
158 //....oooOO0OOooo........oooOO0OOooo........oo << 
159                                                   122 
160 void G4BraggModel::SetParticle(const G4Particl << 123     if(pParticleChange) {
161 {                                              << 124       fParticleChange = 
162   particle = p;                                << 125   reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
163   mass = particle->GetPDGMass();               << 126     } else {
164   spin = particle->GetPDGSpin();               << 127       fParticleChange = new G4ParticleChangeForLoss();
165   G4double q = particle->GetPDGCharge()/CLHEP: << 128     }
166   chargeSquare = q*q;                          << 129   }
167   massRate = mass/CLHEP::proton_mass_c2;       << 
168   ratio = CLHEP::electron_mass_c2/mass;        << 
169 }                                                 130 }
170                                                   131 
171 //....oooOO0OOooo........oooOO0OOooo........oo    132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172                                                   133 
173 G4double G4BraggModel::GetChargeSquareRatio(co    134 G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
174                                             co << 135               const G4Material* mat,
175                                             G4 << 136               G4double kineticEnergy)
176 {                                                 137 {
177   // this method is called only for ions          138   // this method is called only for ions
178   chargeSquare = corr->EffectiveChargeSquareRa << 139   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
179   return chargeSquare;                         << 140   GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
180 }                                              << 141   return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
181                                                << 
182 //....oooOO0OOooo........oooOO0OOooo........oo << 
183                                                << 
184 G4double G4BraggModel::MinEnergyCut(const G4Pa << 
185                                     const G4Ma << 
186 {                                              << 
187   return couple->GetMaterial()->GetIonisation( << 
188 }                                                 142 }
189                                                   143 
190 //....oooOO0OOooo........oooOO0OOooo........oo    144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191                                                   145 
192 G4double G4BraggModel::GetParticleCharge(const    146 G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
193                                          const << 147            const G4Material* mat,
194                                          G4dou << 148            G4double kineticEnergy)
195 {                                                 149 {
196   // this method is called only for ions, so n << 150   // this method is called only for ions
197   return corr->GetParticleCharge(p,mat,kinetic    151   return corr->GetParticleCharge(p,mat,kineticEnergy);
198 }                                                 152 }
199                                                   153 
200 //....oooOO0OOooo........oooOO0OOooo........oo    154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
201                                                   155 
202 G4double G4BraggModel::ComputeCrossSectionPerE    156 G4double G4BraggModel::ComputeCrossSectionPerElectron(
203                                            con    157                                            const G4ParticleDefinition* p,
204                                                   158                                                  G4double kineticEnergy,
205                                                << 159                                                  G4double cutEnergy,
206                                                   160                                                  G4double maxKinEnergy)
207 {                                                 161 {
208   G4double cross = 0.0;                        << 162   G4double cross     = 0.0;
209   const G4double tmax      = MaxSecondaryEnerg << 163   G4double tmax      = MaxSecondaryEnergy(p, kineticEnergy);
210   const G4double maxEnergy = std::min(tmax, ma << 164   G4double maxEnergy = std::min(tmax,maxKinEnergy);
211   const G4double cutEnergy = std::max(cut, low << 165   if(cutEnergy < tmax) {
212   if(cutEnergy < maxEnergy) {                  << 166 
213                                                << 167     G4double energy  = kineticEnergy + mass;
214     const G4double energy  = kineticEnergy + m << 168     G4double energy2 = energy*energy;
215     const G4double energy2 = energy*energy;    << 169     G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
216     const G4double beta2   = kineticEnergy*(ki << 170     cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
217     cross = (maxEnergy - cutEnergy)/(cutEnergy << 
218       - beta2*G4Log(maxEnergy/cutEnergy)/tmax; << 
219                                                << 
220     if( 0.0 < spin ) { cross += 0.5*(maxEnergy << 
221                                                   171 
222     cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 172     cross *= twopi_mc2_rcl2*chargeSquare/beta2;
223   }                                               173   }
224   //   G4cout << "BR: e= " << kineticEnergy << << 174  //   G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy 
225   //          << " tmax= " << tmax << " cross= << 175  //          << " tmax= " << tmax << " cross= " << cross << G4endl;
                                                   >> 176  
226   return cross;                                   177   return cross;
227 }                                                 178 }
228                                                   179 
229 //....oooOO0OOooo........oooOO0OOooo........oo    180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230                                                   181 
231 G4double                                       << 182 G4double G4BraggModel::ComputeCrossSectionPerAtom(
232 G4BraggModel::ComputeCrossSectionPerAtom(const << 183                                            const G4ParticleDefinition* p,
233                                          G4dou << 184                                                  G4double kineticEnergy,
234                                          G4dou << 185              G4double Z, G4double,
235                                          G4dou << 186                                                  G4double cutEnergy,
236                                          G4dou << 187                                                  G4double maxEnergy)
237 {                                                 188 {
238   return                                       << 189   G4double cross = Z*ComputeCrossSectionPerElectron
239     Z*ComputeCrossSectionPerElectron(p,kinetic << 190                                          (p,kineticEnergy,cutEnergy,maxEnergy);
                                                   >> 191   return cross;
240 }                                                 192 }
241                                                   193 
242 //....oooOO0OOooo........oooOO0OOooo........oo    194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243                                                   195 
244 G4double G4BraggModel::CrossSectionPerVolume(c << 196 G4double G4BraggModel::CrossSectionPerVolume(
245                                              c << 197              const G4Material* material,
246                                              G << 198                                            const G4ParticleDefinition* p,
247                                              G << 199                                                  G4double kineticEnergy,
248                                              G << 200                                                  G4double cutEnergy,
                                                   >> 201                                                  G4double maxEnergy)
249 {                                                 202 {
250   return material->GetElectronDensity()        << 203   G4double eDensity = material->GetElectronDensity();
251     *ComputeCrossSectionPerElectron(p,kineticE << 204   G4double cross = eDensity*ComputeCrossSectionPerElectron
                                                   >> 205                                          (p,kineticEnergy,cutEnergy,maxEnergy);
                                                   >> 206   return cross;
252 }                                                 207 }
253                                                   208 
254 //....oooOO0OOooo........oooOO0OOooo........oo    209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255                                                   210 
256 G4double G4BraggModel::ComputeDEDXPerVolume(co    211 G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
257                                             co << 212               const G4ParticleDefinition* p,
258                                             G4 << 213               G4double kineticEnergy,
259                                             G4 << 214               G4double cutEnergy)
260 {                                              << 215 {
261   const G4double tmax = MaxSecondaryEnergy(p,  << 216   G4double tmax  = MaxSecondaryEnergy(p, kineticEnergy);
262   const G4double tkin = kinEnergy/massRate;    << 217   G4double tkin  = kineticEnergy/massRate;
263   const G4double cutEnergy = std::max(cut, low << 
264   G4double dedx  = 0.0;                           218   G4double dedx  = 0.0;
                                                   >> 219   if(tkin > lowestKinEnergy) dedx = DEDX(material, tkin);
                                                   >> 220   else      dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
265                                                   221 
266   // tkin is the scaled energy to proton       << 222   if (cutEnergy < tmax) {
267   if (tkin < lowestKinEnergy) {                << 
268     dedx = DEDX(material, lowestKinEnergy)*std << 
269   } else {                                     << 
270     dedx = DEDX(material, tkin);               << 
271                                                   223 
272     // correction for low cut value using Beth << 224     G4double tau   = kineticEnergy/mass;
273     if (cutEnergy < tmax) {                    << 225     G4double gam   = tau + 1.0;
274       const G4double tau = kinEnergy/mass;     << 226     G4double bg2   = tau * (tau+2.0);
275       const G4double x = cutEnergy/tmax;       << 227     G4double beta2 = bg2/(gam*gam);
                                                   >> 228     G4double x     = cutEnergy/tmax;
276                                                   229 
277       dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/ << 230     dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
278   CLHEP::twopi_mc2_rcl2 * material->GetElectro << 231           * (material->GetElectronDensity())/beta2;
279     }                                          << 
280   }                                               232   }
281   dedx = std::max(dedx, 0.0) * chargeSquare;   << 
282                                                   233 
283   //G4cout << "E(MeV)= " << tkin/MeV << " dedx << 234   // now compute the total ionization loss
284   //         << "  " << material->GetName() << << 235 
                                                   >> 236   if (dedx < 0.0) dedx = 0.0 ;
                                                   >> 237 
                                                   >> 238   dedx *= chargeSquare;
                                                   >> 239 
285   return dedx;                                    240   return dedx;
286 }                                                 241 }
287                                                   242 
288 //....oooOO0OOooo........oooOO0OOooo........oo    243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289                                                   244 
290 void G4BraggModel::SampleSecondaries(std::vect << 245 void G4BraggModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
291                                      const G4M << 246           const G4DynamicParticle* dp,
292                                      const G4D << 247           G4double& eloss,
293                                      G4double  << 248           G4double&,
294                                      G4double  << 249           G4double length)
295 {                                              << 250 {
296   const G4double tmax = MaxSecondaryKinEnergy( << 251   if(nuclearStopping) {
297   const G4double xmax = std::min(tmax, maxEner << 252 
298   const G4double xmin = std::max(lowestKinEner << 253     G4double preKinEnergy = dp->GetKineticEnergy();
299   if(xmin >= xmax) { return; }                 << 254     G4double e = preKinEnergy - eloss*0.5;
                                                   >> 255     if(e < 0.0) e = preKinEnergy*0.5;
                                                   >> 256     G4double nloss = length*corr->NuclearDEDX(dp->GetDefinition(),
                                                   >> 257                 couple->GetMaterial(),
                                                   >> 258                 e,false);
                                                   >> 259 
                                                   >> 260     // too big energy loss
                                                   >> 261     if(eloss + nloss > preKinEnergy) {
                                                   >> 262       nloss *= (preKinEnergy/(eloss + nloss));
                                                   >> 263       eloss = preKinEnergy;
                                                   >> 264     } else {
                                                   >> 265       eloss += nloss;
                                                   >> 266     }
                                                   >> 267     /*
                                                   >> 268     G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
                                                   >> 269          << " de= " << eloss << " NIEL= " << nloss 
                                                   >> 270      << " dynQ= " << dp->GetCharge()/eplus << G4endl;
                                                   >> 271     */
                                                   >> 272     fParticleChange->ProposeNonIonizingEnergyDeposit(nloss);
                                                   >> 273   }
                                                   >> 274 }
                                                   >> 275 
                                                   >> 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 277 
                                                   >> 278 void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
                                                   >> 279              const G4MaterialCutsCouple*,
                                                   >> 280              const G4DynamicParticle* dp,
                                                   >> 281              G4double xmin,
                                                   >> 282              G4double maxEnergy)
                                                   >> 283 {
                                                   >> 284   G4double tmax = MaxSecondaryKinEnergy(dp);
                                                   >> 285   G4double xmax = std::min(tmax, maxEnergy);
                                                   >> 286   if(xmin >= xmax) return;
300                                                   287 
301   G4double kineticEnergy = dp->GetKineticEnerg    288   G4double kineticEnergy = dp->GetKineticEnergy();
302   const G4double energy  = kineticEnergy + mas << 289   G4double energy  = kineticEnergy + mass;
303   const G4double energy2 = energy*energy;      << 290   G4double energy2 = energy*energy;
304   const G4double beta2   = kineticEnergy*(kine << 291   G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
305   const G4double grej = 1.0;                   << 292   G4double grej    = 1.0;
306   G4double deltaKinEnergy, f;                     293   G4double deltaKinEnergy, f;
307                                                   294 
308   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 295   G4ThreeVector direction = dp->GetMomentumDirection();
309   G4double rndm[2];                            << 
310                                                   296 
311   // sampling follows ...                         297   // sampling follows ...
312   do {                                            298   do {
313     rndmEngineMod->flatArray(2, rndm);         << 299     G4double q = G4UniformRand();
314     deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rn << 300     deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
315                                                   301 
316     f = 1.0 - beta2*deltaKinEnergy/tmax;          302     f = 1.0 - beta2*deltaKinEnergy/tmax;
317                                                   303 
318     if(f > grej) {                                304     if(f > grej) {
319         G4cout << "G4BraggModel::SampleSeconda    305         G4cout << "G4BraggModel::SampleSecondary Warning! "
320                << "Majorant " << grej << " < "    306                << "Majorant " << grej << " < "
321                << f << " for e= " << deltaKinE    307                << f << " for e= " << deltaKinEnergy
322                << G4endl;                         308                << G4endl;
323     }                                             309     }
324                                                   310 
325     // Loop checking, 03-Aug-2015, Vladimir Iv << 311   } while( grej*G4UniformRand() >= f );
326   } while( grej*rndm[1] >= f );                << 
327                                                << 
328   G4ThreeVector deltaDirection;                << 
329                                                << 
330   if(UseAngularGeneratorFlag()) {              << 
331     const G4Material* mat =  couple->GetMateri << 
332     G4int Z = SelectRandomAtomNumber(mat);     << 
333                                                   312 
334     deltaDirection =                           << 313   G4double deltaMomentum =
335       GetAngularDistribution()->SampleDirectio << 314            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
                                                   >> 315   G4double totMomentum = energy*sqrt(beta2);
                                                   >> 316   G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
                                                   >> 317                                    (deltaMomentum * totMomentum);
                                                   >> 318   if(cost > 1.0) cost = 1.0;
                                                   >> 319   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
336                                                   320 
337   } else {                                     << 321   G4double phi = twopi * G4UniformRand() ;
338                                                << 
339     G4double deltaMomentum =                   << 
340       std::sqrt(deltaKinEnergy * (deltaKinEner << 
341     G4double cost = deltaKinEnergy * (energy + << 
342       (deltaMomentum * dp->GetTotalMomentum()) << 
343     if(cost > 1.0) { cost = 1.0; }             << 
344     G4double sint = std::sqrt((1.0 - cost)*(1. << 
345                                                << 
346     G4double phi = twopi*rndmEngineMod->flat() << 
347                                                << 
348     deltaDirection.set(sint*std::cos(phi),sint << 
349     deltaDirection.rotateUz(dp->GetMomentumDir << 
350   }                                            << 
351                                                   322 
352   // create G4DynamicParticle object for delta << 323   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
353   auto delta = new G4DynamicParticle(theElectr << 324   deltaDirection.rotateUz(direction);
354                                                   325 
355   // Change kinematics of primary particle        326   // Change kinematics of primary particle
356   kineticEnergy -= deltaKinEnergy;             << 327   kineticEnergy       -= deltaKinEnergy;
357   G4ThreeVector finalP = dp->GetMomentum() - d << 328   G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
358   finalP = finalP.unit();                      << 329   finalP               = finalP.unit();
359                                                   330   
360   fParticleChange->SetProposedKineticEnergy(ki    331   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
361   fParticleChange->SetProposedMomentumDirectio    332   fParticleChange->SetProposedMomentumDirection(finalP);
362                                                   333 
363   vdp->push_back(delta);                       << 334   // create G4DynamicParticle object for delta ray
364 }                                              << 335   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
365                                                << 336                deltaKinEnergy);
366 //....oooOO0OOooo........oooOO0OOooo........oo << 
367                                                   337 
368 G4double G4BraggModel::MaxSecondaryEnergy(cons << 338   vdp->push_back(delta);
369                                           G4do << 
370 {                                              << 
371   if(pd != particle) { SetParticle(pd); }      << 
372   G4double tau  = kinEnergy/mass;              << 
373   G4double tmax = 2.0*electron_mass_c2*tau*(ta << 
374                   (1. + 2.0*(tau + 1.)*ratio + << 
375   return tmax;                                 << 
376 }                                                 339 }
377                                                   340 
378 //....oooOO0OOooo........oooOO0OOooo........oo    341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379                                                   342 
380 void G4BraggModel::HasMaterial(const G4Materia << 343 G4bool G4BraggModel::HasMaterial(const G4Material* material)
381 {                                                 344 {
382   const G4String& chFormula = mat->GetChemical << 345   const size_t numberOfMolecula = 11 ;
383   if(chFormula.empty()) { return; }            << 346   SetMoleculaNumber(numberOfMolecula) ;
                                                   >> 347   G4String chFormula = material->GetChemicalFormula() ;
384                                                   348 
385   // ICRU Report N49, 1993. Power's model for  << 349   // ICRU Report N49, 1993. Power's model for He.
386   static const G4int numberOfMolecula = 11;    << 350   static G4String molName[numberOfMolecula] = {
387   static const G4String molName[numberOfMolecu << 
388     "Al_2O_3",                 "CO_2",            351     "Al_2O_3",                 "CO_2",                      "CH_4",
389     "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol    352     "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene",  "(C_8H_8)_N",
390     "C_3H_8",                  "SiO_2",           353     "C_3H_8",                  "SiO_2",                     "H_2O",
391     "H_2O-Gas",                "Graphite" } ;     354     "H_2O-Gas",                "Graphite" } ;
392                                                   355 
                                                   >> 356   // Special treatment for water in gas state
                                                   >> 357   const G4State theState = material->GetState() ;
                                                   >> 358 
                                                   >> 359   if( theState == kStateGas && "H_2O" == chFormula) {
                                                   >> 360     chFormula = G4String("H_2O-Gas");
                                                   >> 361   }
                                                   >> 362 
393   // Search for the material in the table         363   // Search for the material in the table
394   for (G4int i=0; i<numberOfMolecula; ++i) {   << 364   for (size_t i=0; i<numberOfMolecula; i++) {
395     if (chFormula == molName[i]) {             << 365       if (chFormula == molName[i]) {
396       iMolecula = i;                           << 366         SetMoleculaNumber(i) ;
397       return;                                  << 367   return true ;
398     }                                          << 368       }
399   }                                               369   }
400   return;                                      << 370   return false ;
401 }                                                 371 }
402                                                   372 
403 //....oooOO0OOooo........oooOO0OOooo........oo    373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404                                                   374 
405 G4double G4BraggModel::StoppingPower(const G4M    375 G4double G4BraggModel::StoppingPower(const G4Material* material,
406                                      G4double  << 376                                            G4double kineticEnergy) 
407 {                                                 377 {
408   G4double ionloss = 0.0 ;                        378   G4double ionloss = 0.0 ;
409                                                   379 
410   if (iMolecula >= 0) {                        << 380   if (iMolecula < 11) {
411                                                   381   
412     // The data and the fit from:                 382     // The data and the fit from: 
413     // ICRU Report N49, 1993. Ziegler's model     383     // ICRU Report N49, 1993. Ziegler's model for protons.
414     // Proton kinetic energy for parametrisati    384     // Proton kinetic energy for parametrisation (keV/amu)  
415                                                   385 
416     G4double T = kineticEnergy/(keV*protonMass    386     G4double T = kineticEnergy/(keV*protonMassAMU) ; 
417                                                   387 
418     static const G4float a[11][5] = {          << 388      static G4double a[11][5] = {
419    {1.187E+1f, 1.343E+1f, 1.069E+4f, 7.723E+2f << 389    {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
420    {7.802E+0f, 8.814E+0f, 8.303E+3f, 7.446E+2f << 390    {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3}, 
421    {7.294E+0f, 8.284E+0f, 5.010E+3f, 4.544E+2f << 391    {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3}, 
422    {8.646E+0f, 9.800E+0f, 7.066E+3f, 4.581E+2f << 392    {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3}, 
423    {1.286E+1f, 1.462E+1f, 5.625E+3f, 2.621E+3f << 393    {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2}, 
424    {3.229E+1f, 3.696E+1f, 8.918E+3f, 3.244E+3f << 394    {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1}, 
425    {1.604E+1f, 1.825E+1f, 6.967E+3f, 2.307E+3f << 395    {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2}, 
426    {8.049E+0f, 9.099E+0f, 9.257E+3f, 3.846E+2f << 396    {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
427    {4.015E+0f, 4.542E+0f, 3.955E+3f, 4.847E+2f << 397    {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3}, 
428    {4.571E+0f, 5.173E+0f, 4.346E+3f, 4.779E+2f << 398    {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
429    {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f << 399    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
430                                                << 400 
431     static const G4float atomicWeight[11] = {  << 401      static G4double atomicWeight[11] = {
432     101.96128f, 44.0098f, 16.0426f, 28.0536f,  << 402     101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
433     104.1512f, 44.665f, 60.0843f, 18.0152f, 18 << 403     104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};       
434                                                   404 
435     if ( T < 10.0 ) {                             405     if ( T < 10.0 ) {
436       ionloss = ((G4double)(a[iMolecula][0]))  << 406       ionloss = a[iMolecula][0] * sqrt(T) ;
437                                                   407     
438     } else if ( T < 10000.0 ) {                   408     } else if ( T < 10000.0 ) {
439       G4double x1 = (G4double)(a[iMolecula][1] << 409       G4double slow  = a[iMolecula][1] * pow(T, 0.45) ;
440       G4double x2 = (G4double)(a[iMolecula][2] << 410       G4double shigh = log( 1.0 + a[iMolecula][3]/T  
441       G4double x3 = (G4double)(a[iMolecula][3] << 411                      + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
442       G4double x4 = (G4double)(a[iMolecula][4] << 
443       G4double slow  = x1 * G4Exp(G4Log(T)* 0. << 
444       G4double shigh = G4Log( 1.0 + x3/T  + x4 << 
445       ionloss = slow*shigh / (slow + shigh) ;     412       ionloss = slow*shigh / (slow + shigh) ;     
446     }                                             413     } 
447                                                   414 
448     ionloss = std::max(ionloss, 0.0);          << 415     if ( ionloss < 0.0) ionloss = 0.0 ;
449     if ( 10 == iMolecula ) {                      416     if ( 10 == iMolecula ) { 
450       static const G4double invLog10 = 1.0/G4L << 
451                                                << 
452       if (T < 100.0) {                            417       if (T < 100.0) {
453         ionloss *= (1.0+0.023+0.0066*G4Log(T)* << 418   ionloss *= (1.0+0.023+0.0066*log10(T));  
454       }                                           419       }
455       else if (T < 700.0) {                       420       else if (T < 700.0) {   
456         ionloss *=(1.0+0.089-0.0248*G4Log(T-99 << 421   ionloss *=(1.0+0.089-0.0248*log10(T-99.));
457       }                                           422       } 
458       else if (T < 10000.0) {                     423       else if (T < 10000.0) {    
459         ionloss *=(1.0+0.089-0.0248*G4Log(700. << 424   ionloss *=(1.0+0.089-0.0248*log10(700.-99.));
460       }                                           425       }
461     }                                             426     }
462     ionloss /= (G4double)atomicWeight[iMolecul << 427     ionloss /= atomicWeight[iMolecula];
463                                                   428 
464   // pure material (normally not the case for     429   // pure material (normally not the case for this function)
465   } else if(1 == (material->GetNumberOfElement    430   } else if(1 == (material->GetNumberOfElements())) {
466     G4double z = material->GetZ() ;               431     G4double z = material->GetZ() ;
467     ionloss = ElectronicStoppingPower( z, kine    432     ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;  
468   }                                               433   }
469                                                   434   
470   return ionloss;                                 435   return ionloss;
471 }                                                 436 }
472                                                   437 
473 //....oooOO0OOooo........oooOO0OOooo........oo    438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474                                                   439 
475 G4double G4BraggModel::ElectronicStoppingPower    440 G4double G4BraggModel::ElectronicStoppingPower(G4double z,
476                                                   441                                                G4double kineticEnergy) const
477 {                                                 442 {
478   G4double ionloss ;                              443   G4double ionloss ;
479   G4int i = std::min(std::max(G4lrint(z)-1,0), << 444   G4int i = G4int(z)-1 ;  // index of atom
480                                                << 445   if(i < 0)  i = 0 ;
                                                   >> 446   if(i > 91) i = 91 ;
                                                   >> 447   
481   // The data and the fit from:                   448   // The data and the fit from: 
482   // ICRU Report 49, 1993. Ziegler's type of p    449   // ICRU Report 49, 1993. Ziegler's type of parametrisations.
483   // Proton kinetic energy for parametrisation    450   // Proton kinetic energy for parametrisation (keV/amu)  
484                                                   451 
485   G4double T = kineticEnergy/(keV*protonMassAM    452   G4double T = kineticEnergy/(keV*protonMassAMU) ; 
486                                                   453   
487   static const G4float a[92][5] = {            << 454   static G4double a[92][5] = {
488    {1.254E+0f, 1.440E+0f, 2.426E+2f, 1.200E+4f << 455    {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
489    {1.229E+0f, 1.397E+0f, 4.845E+2f, 5.873E+3f << 456    {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
490    {1.411E+0f, 1.600E+0f, 7.256E+2f, 3.013E+3f << 457    {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
491    {2.248E+0f, 2.590E+0f, 9.660E+2f, 1.538E+2f << 458    {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
492    {2.474E+0f, 2.815E+0f, 1.206E+3f, 1.060E+3f << 459    {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
493    {2.631E+0f, 2.601E+0f, 1.701E+3f, 1.279E+3f << 460    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
494    {2.954E+0f, 3.350E+0f, 1.683E+3f, 1.900E+3f << 461    {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
495    {2.652E+0f, 3.000E+0f, 1.920E+3f, 2.000E+3f << 462    {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
496    {2.085E+0f, 2.352E+0f, 2.157E+3f, 2.634E+3f << 463    {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
497    {1.951E+0f, 2.199E+0f, 2.393E+3f, 2.699E+3f << 464    {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
498        // Z= 11-20                                465        // Z= 11-20
499    {2.542E+0f, 2.869E+0f, 2.628E+3f, 1.854E+3f << 466    {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
500    {3.791E+0f, 4.293E+0f, 2.862E+3f, 1.009E+3f << 467    {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
501    {4.154E+0f, 4.739E+0f, 2.766E+3f, 1.645E+2f << 468    {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
502    {4.914E+0f, 5.598E+0f, 3.193E+3f, 2.327E+2f << 469    {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
503    {3.232E+0f, 3.647E+0f, 3.561E+3f, 1.560E+3f << 470    {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
504    {3.447E+0f, 3.891E+0f, 3.792E+3f, 1.219E+3f << 471    {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
505    {5.301E+0f, 6.008E+0f, 3.969E+3f, 6.451E+2f << 472    {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
506    {5.731E+0f, 6.500E+0f, 4.253E+3f, 5.300E+2f << 473    {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
507    {5.152E+0f, 5.833E+0f, 4.482E+3f, 5.457E+2f << 474    {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
508    {5.521E+0f, 6.252E+0f, 4.710E+3f, 5.533E+2f << 475    {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
509        // Z= 21-30                                476        // Z= 21-30
510    {5.201E+0f, 5.884E+0f, 4.938E+3f, 5.609E+2f << 477    {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
511    {4.858E+0f, 5.489E+0f, 5.260E+3f, 6.511E+2f << 478    {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
512    {4.479E+0f, 5.055E+0f, 5.391E+3f, 9.523E+2f << 479    {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
513    {3.983E+0f, 4.489E+0f, 5.616E+3f, 1.336E+3f << 480    {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
514    {3.469E+0f, 3.907E+0f, 5.725E+3f, 1.461E+3f << 481    {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
515    {3.519E+0f, 3.963E+0f, 6.065E+3f, 1.243E+3f << 482    {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
516    {3.140E+0f, 3.535E+0f, 6.288E+3f, 1.372E+3f << 483    {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
517    {3.553E+0f, 4.004E+0f, 6.205E+3f, 5.551E+2f << 484    {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
518    {3.696E+0f, 4.194E+0f, 4.649E+3f, 8.113E+1f << 485    {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
519    {4.210E+0f, 4.750E+0f, 6.953E+3f, 2.952E+2f << 486    {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
520        // Z= 31-40                                487        // Z= 31-40
521    {5.041E+0f, 5.697E+0f, 7.173E+3f, 2.026E+2f << 488    {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
522    {5.554E+0f, 6.300E+0f, 6.496E+3f, 1.100E+2f << 489    {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
523    {5.323E+0f, 6.012E+0f, 7.611E+3f, 2.925E+2f << 490    {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
524    {5.874E+0f, 6.656E+0f, 7.395E+3f, 1.175E+2f << 491    {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
525    {6.658E+0f, 7.536E+0f, 7.694E+3f, 2.223E+2f << 492    {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
526    {6.413E+0f, 7.240E+0f, 1.185E+4f, 1.537E+2f << 493    {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
527    {5.694E+0f, 6.429E+0f, 8.478E+3f, 2.929E+2f << 494    {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
528    {6.339E+0f, 7.159E+0f, 8.693E+3f, 3.303E+2f << 495    {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
529    {6.407E+0f, 7.234E+0f, 8.907E+3f, 3.678E+2f << 496    {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
530    {6.734E+0f, 7.603E+0f, 9.120E+3f, 4.052E+2f << 497    {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
531        // Z= 41-50                                498        // Z= 41-50
532    {6.901E+0f, 7.791E+0f, 9.333E+3f, 4.427E+2f << 499    {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
533    {6.424E+0f, 7.248E+0f, 9.545E+3f, 4.802E+2f << 500    {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
534    {6.799E+0f, 7.671E+0f, 9.756E+3f, 5.176E+2f << 501    {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
535    {6.109E+0f, 6.887E+0f, 9.966E+3f, 5.551E+2f << 502    {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
536    {5.924E+0f, 6.677E+0f, 1.018E+4f, 5.925E+2f << 503    {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
537    {5.238E+0f, 5.900E+0f, 1.038E+4f, 6.300E+2f << 504    {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
538    // {5.623f,    6.354f,    7160.0f,   337.6f << 505    // {5.623,    6.354,    7160.0,   337.6,    0.013940}, // Ag Ziegler77
539    {5.345E+0f, 6.038E+0f, 6.790E+3f, 3.978E+2f << 506    {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2}, // Ag ICRU49
540    {5.814E+0f, 6.554E+0f, 1.080E+4f, 3.555E+2f << 507    {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
541    {6.229E+0f, 7.024E+0f, 1.101E+4f, 3.709E+2f << 508    {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
542    {6.409E+0f, 7.227E+0f, 1.121E+4f, 3.864E+2f << 509    {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
543        // Z= 51-60                                510        // Z= 51-60
544    {7.500E+0f, 8.480E+0f, 8.608E+3f, 3.480E+2f << 511    {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
545    {6.979E+0f, 7.871E+0f, 1.162E+4f, 3.924E+2f << 512    {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
546    {7.725E+0f, 8.716E+0f, 1.183E+4f, 3.948E+2f << 513    {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
547    {8.337E+0f, 9.425E+0f, 1.051E+4f, 2.696E+2f << 514    {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
548    {7.287E+0f, 8.218E+0f, 1.223E+4f, 3.997E+2f << 515    {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
549    {7.899E+0f, 8.911E+0f, 1.243E+4f, 4.021E+2f << 516    {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
550    {8.041E+0f, 9.071E+0f, 1.263E+4f, 4.045E+2f << 517    {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
551    {7.488E+0f, 8.444E+0f, 1.283E+4f, 4.069E+2f << 518    {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
552    {7.291E+0f, 8.219E+0f, 1.303E+4f, 4.093E+2f << 519    {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
553    {7.098E+0f, 8.000E+0f, 1.323E+4f, 4.118E+2f << 520    {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
554        // Z= 61-70                                521        // Z= 61-70
555    {6.909E+0f, 7.786E+0f, 1.343E+4f, 4.142E+2f << 522    {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
556    {6.728E+0f, 7.580E+0f, 1.362E+4f, 4.166E+2f << 523    {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
557    {6.551E+0f, 7.380E+0f, 1.382E+4f, 4.190E+2f << 524    {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
558    {6.739E+0f, 7.592E+0f, 1.402E+4f, 4.214E+2f << 525    {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
559    {6.212E+0f, 6.996E+0f, 1.421E+4f, 4.239E+2f << 526    {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
560    {5.517E+0f, 6.210E+0f, 1.440E+4f, 4.263E+2f << 527    {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
561    {5.220E+0f, 5.874E+0f, 1.460E+4f, 4.287E+2f << 528    {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
562    {5.071E+0f, 5.706E+0f, 1.479E+4f, 4.330E+2f << 529    {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
563    {4.926E+0f, 5.542E+0f, 1.498E+4f, 4.335E+2f << 530    {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
564    {4.788E+0f, 5.386E+0f, 1.517E+4f, 4.359E+2f << 531    {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
565        // Z= 71-80                                532        // Z= 71-80
566    {4.893E+0f, 5.505E+0f, 1.536E+4f, 4.384E+2f << 533    {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
567    {5.028E+0f, 5.657E+0f, 1.555E+4f, 4.408E+2f << 534    {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
568    {4.738E+0f, 5.329E+0f, 1.574E+4f, 4.432E+2f << 535    {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
569    {4.587E+0f, 5.160E+0f, 1.541E+4f, 4.153E+2f << 536    {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
570    {5.201E+0f, 5.851E+0f, 1.612E+4f, 4.416E+2f << 537    {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
571    {5.071E+0f, 5.704E+0f, 1.630E+4f, 4.409E+2f << 538    {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
572    {4.946E+0f, 5.563E+0f, 1.649E+4f, 4.401E+2f << 539    {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
573    {4.477E+0f, 5.034E+0f, 1.667E+4f, 4.393E+2f << 540    {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
574    //  {4.856f,    5.460f,    18320.0f,  438.5 << 541    //  {4.856,    5.460,    18320.0,  438.5,    0.002542}, //Ziegler77
575    {4.844E+0f, 5.458E+0f, 7.852E+3f, 9.758E+2f << 542    {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2}, //ICRU49
576    {4.307E+0f, 4.843E+0f, 1.704E+4f, 4.878E+2f << 543    {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
577        // Z= 81-90                                544        // Z= 81-90
578    {4.723E+0f, 5.311E+0f, 1.722E+4f, 5.370E+2f << 545    {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
579    {5.319E+0f, 5.982E+0f, 1.740E+4f, 5.863E+2f << 546    {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
580    {5.956E+0f, 6.700E+0f, 1.780E+4f, 6.770E+2f << 547    {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
581    {6.158E+0f, 6.928E+0f, 1.777E+4f, 5.863E+2f << 548    {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
582    {6.203E+0f, 6.979E+0f, 1.795E+4f, 5.863E+2f << 549    {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
583    {6.181E+0f, 6.954E+0f, 1.812E+4f, 5.863E+2f << 550    {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
584    {6.949E+0f, 7.820E+0f, 1.830E+4f, 5.863E+2f << 551    {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
585    {7.506E+0f, 8.448E+0f, 1.848E+4f, 5.863E+2f << 552    {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
586    {7.648E+0f, 8.609E+0f, 1.866E+4f, 5.863E+2f << 553    {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
587    {7.711E+0f, 8.679E+0f, 1.883E+4f, 5.863E+2f << 554    {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
588        // Z= 91-92                                555        // Z= 91-92
589    {7.407E+0f, 8.336E+0f, 1.901E+4f, 5.863E+2f << 556    {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
590    {7.290E+0f, 8.204E+0f, 1.918E+4f, 5.863E+2f << 557    {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
591   };                                              558   };
592                                                   559 
593   G4double fac = 1.0 ;                            560   G4double fac = 1.0 ;
594                                                   561 
595     // Carbon specific case for E < 40 keV        562     // Carbon specific case for E < 40 keV
596   if ( T < 40.0 && 5 == i) {                      563   if ( T < 40.0 && 5 == i) {
597     fac = std::sqrt(T*0.025);                  << 564     fac = sqrt(T/40.0) ;
598     T = 40.0;                                  << 565     T = 40.0 ;  
599                                                   566 
600     // Free electron gas model                    567     // Free electron gas model
601   } else if ( T < 10.0 ) {                        568   } else if ( T < 10.0 ) { 
602     fac = std::sqrt(T*0.1) ;                   << 569     fac = sqrt(T*0.1) ;
603     T = 10.0;                                  << 570     T =10.0 ;
604   }                                               571   }
605                                                   572 
606   // Main parametrisation                         573   // Main parametrisation
607   G4double x1 = (G4double)(a[i][1]);           << 574   G4double slow  = a[i][1] * pow(T, 0.45) ;
608   G4double x2 = (G4double)(a[i][2]);           << 575   G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
609   G4double x3 = (G4double)(a[i][3]);           << 576   ionloss = slow*shigh*fac / (slow + shigh) ;     
610   G4double x4 = (G4double)(a[i][4]);           << 
611   G4double slow  = x1 * G4Exp(G4Log(T) * 0.45) << 
612   G4double shigh = G4Log( 1.0 + x3/T + x4*T )  << 
613   ionloss = slow*shigh*fac / (slow + shigh);   << 
614                                                   577   
615   ionloss = std::max(ionloss, 0.0);            << 578   if ( ionloss < 0.0) ionloss = 0.0 ;
616                                                   579   
617   return ionloss;                                 580   return ionloss;
618 }                                                 581 }
619                                                   582 
620 //....oooOO0OOooo........oooOO0OOooo........oo    583 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
621                                                   584 
622 G4double G4BraggModel::DEDX(const G4Material*  << 585 G4double G4BraggModel::DEDX(const G4Material* material,
                                                   >> 586                                   G4double kineticEnergy) 
623 {                                                 587 {
624   G4double eloss = 0.0;                           588   G4double eloss = 0.0;
625                                                << 589   const G4int numberOfElements = material->GetNumberOfElements();
626   // check DB                                  << 
627   if(material != currentMaterial) {            << 
628     currentMaterial = material;                << 
629     baseMaterial = material->GetBaseMaterial() << 
630       ? material->GetBaseMaterial() : material << 
631     iPSTAR    = -1;                            << 
632     iMolecula = -1;                            << 
633     iICRU90 = fICRU90 ? fICRU90->GetIndex(base << 
634                                                << 
635     if(iICRU90 < 0) {                          << 
636       iPSTAR = fPSTAR->GetIndex(baseMaterial); << 
637       if(iPSTAR < 0) { HasMaterial(baseMateria << 
638     }                                          << 
639     //G4cout << "%%% " <<material->GetName() < << 
640     //       << iMolecula << "  iPSTAR= " << i << 
641     //       << "  iICRU90= " << iICRU90<< G4e << 
642   }                                            << 
643                                                << 
644   // ICRU90 parameterisation                   << 
645   if(iICRU90 >= 0) {                           << 
646     return fICRU90->GetElectronicDEDXforProton << 
647       *material->GetDensity();                 << 
648   }                                            << 
649   // PSTAR parameterisation                    << 
650   if( iPSTAR >= 0 ) {                          << 
651     return fPSTAR->GetElectronicDEDX(iPSTAR, k << 
652       *material->GetDensity();                 << 
653                                                << 
654   }                                            << 
655   const std::size_t numberOfElements = materia << 
656   const G4double* theAtomicNumDensityVector =     590   const G4double* theAtomicNumDensityVector =
657                                  material->Get    591                                  material->GetAtomicNumDensityVector();
658                                                   592   
                                                   >> 593   // compaund material with parametrisation
                                                   >> 594   G4int iNist = pstar.GetIndex(material);
                                                   >> 595 
                                                   >> 596   if( iNist >= 0 ) {
                                                   >> 597     return pstar.GetElectronicDEDX(iNist, kineticEnergy)*material->GetDensity();
659                                                   598 
660   if(iMolecula >= 0) {                         << 599   } else if( HasMaterial(material) ) {
661     eloss = StoppingPower(baseMaterial, kineti << 600 
                                                   >> 601     eloss = StoppingPower(material, kineticEnergy)*
662                           material->GetDensity    602                           material->GetDensity()/amu;
663                                                   603 
664     // Pure material ICRU49 paralmeterisation  << 604   // pure material
665   } else if(1 == numberOfElements) {              605   } else if(1 == numberOfElements) {
666                                                   606 
667     G4double z = material->GetZ();                607     G4double z = material->GetZ();
668     eloss = ElectronicStoppingPower(z, kinetic    608     eloss = ElectronicStoppingPower(z, kineticEnergy)
669                                * (material->Ge    609                                * (material->GetTotNbOfAtomsPerVolume());
670                                                   610 
671                                                   611 
672   // Experimental data exist only for kinetic     612   // Experimental data exist only for kinetic energy 125 keV
673   } else if( MolecIsInZiegler1988(material) )     613   } else if( MolecIsInZiegler1988(material) ) { 
674                                                   614 
675     // Loop over elements - calculation based  << 615   // Cycle over elements - calculation based on Bragg's rule 
676     G4double eloss125 = 0.0 ;                     616     G4double eloss125 = 0.0 ;
677     const G4ElementVector* theElementVector =     617     const G4ElementVector* theElementVector =
678                            material->GetElemen    618                            material->GetElementVector();
679                                                   619   
680     //  Loop for the elements in the material  << 620     //  loop for the elements in the material
681     for (std::size_t i=0; i<numberOfElements;  << 621     for (G4int i=0; i<numberOfElements; i++) {
682       const G4Element* element = (*theElementV    622       const G4Element* element = (*theElementVector)[i] ;
683       G4double z = element->GetZ() ;              623       G4double z = element->GetZ() ;
684       eloss    += ElectronicStoppingPower(z,ki    624       eloss    += ElectronicStoppingPower(z,kineticEnergy)
685                                     * theAtomi    625                                     * theAtomicNumDensityVector[i] ;
686       eloss125 += ElectronicStoppingPower(z,12    626       eloss125 += ElectronicStoppingPower(z,125.0*keV)
687                                     * theAtomi    627                                     * theAtomicNumDensityVector[i] ;
688     }                                             628     }      
689                                                   629 
690     // Chemical factor is taken into account      630     // Chemical factor is taken into account
691     if (eloss125 > 0.0) {                      << 631     eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
692       eloss *= ChemicalFactor(kineticEnergy, e << 
693     }                                          << 
694                                                   632  
695   // Brugg's rule calculation                     633   // Brugg's rule calculation
696   } else {                                        634   } else {
697     const G4ElementVector* theElementVector =     635     const G4ElementVector* theElementVector =
698                            material->GetElemen    636                            material->GetElementVector() ;
699                                                   637   
700     //  loop for the elements in the material     638     //  loop for the elements in the material
701     for (std::size_t i=0; i<numberOfElements;  << 639     for (G4int i=0; i<numberOfElements; i++)
702     {                                             640     {
703       const G4Element* element = (*theElementV    641       const G4Element* element = (*theElementVector)[i] ;
704       eloss   += ElectronicStoppingPower(eleme    642       eloss   += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
705                                    * theAtomic    643                                    * theAtomicNumDensityVector[i];
706     }                                             644     }      
707   }                                               645   }
708   return eloss*theZieglerFactor;                  646   return eloss*theZieglerFactor;
709 }                                                 647 }
710                                                   648 
711 //....oooOO0OOooo........oooOO0OOooo........oo    649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
712                                                   650 
713 G4bool G4BraggModel::MolecIsInZiegler1988(cons    651 G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material) 
714 {                                                 652 {
715   // The list of molecules from                   653   // The list of molecules from
716   // J.F.Ziegler and J.M.Manoyan, The stopping    654   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
717   // Nucl. Inst. & Meth. in Phys. Res. B35 (19    655   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
718                                                   656   
719   G4String myFormula = G4String(" ") ;            657   G4String myFormula = G4String(" ") ;
720   const G4String chFormula = material->GetChem    658   const G4String chFormula = material->GetChemicalFormula() ;
721   if (myFormula == chFormula ) { return false; << 659   if (myFormula == chFormula ) return false ;
722                                                   660   
723   //  There are no evidence for difference of     661   //  There are no evidence for difference of stopping power depended on
724   //  phase of the compound except for water.     662   //  phase of the compound except for water. The stopping power of the 
725   //  water in gas phase can be predicted usin    663   //  water in gas phase can be predicted using Bragg's rule.
726   //                                              664   //  
727   //  No chemical factor for water-gas            665   //  No chemical factor for water-gas 
728                                                   666    
729   myFormula = G4String("H_2O") ;                  667   myFormula = G4String("H_2O") ;
730   const G4State theState = material->GetState(    668   const G4State theState = material->GetState() ;
731   if( theState == kStateGas && myFormula == ch    669   if( theState == kStateGas && myFormula == chFormula) return false ;
732                                                   670     
                                                   >> 671   const size_t numberOfMolecula = 53 ;
733                                                   672 
734   // The coffecient from Table.4 of Ziegler &     673   // The coffecient from Table.4 of Ziegler & Manoyan
735   static const G4float HeEff = 2.8735f;        << 674   const G4double HeEff = 2.8735 ;
736                                                   675   
737   static const std::size_t numberOfMolecula =  << 676   static G4String nameOfMol[numberOfMolecula] = {
738   static const G4String nameOfMol[numberOfMole << 
739     "H_2O",      "C_2H_4O",    "C_3H_6O",  "C_    677     "H_2O",      "C_2H_4O",    "C_3H_6O",  "C_2H_2",             "C_H_3OH",
740     "C_2H_5OH",  "C_3H_7OH",   "C_3H_4",   "NH    678     "C_2H_5OH",  "C_3H_7OH",   "C_3H_4",   "NH_3",               "C_14H_10",
741     "C_6H_6",    "C_4H_10",    "C_4H_6",   "C_    679     "C_6H_6",    "C_4H_10",    "C_4H_6",   "C_4H_8O",            "CCl_4",
742     "CF_4",      "C_6H_8",     "C_6H_12",  "C_    680     "CF_4",      "C_6H_8",     "C_6H_12",  "C_6H_10O",           "C_6H_10",
743     "C_8H_16",   "C_5H_10",    "C_5H_8",   "C_    681     "C_8H_16",   "C_5H_10",    "C_5H_8",   "C_3H_6-Cyclopropane","C_2H_4F_2",
744     "C_2H_2F_2", "C_4H_8O_2",  "C_2H_6",   "C_    682     "C_2H_2F_2", "C_4H_8O_2",  "C_2H_6",   "C_2F_6",             "C_2H_6O",
745     "C_3H_6O",   "C_4H_10O",   "C_2H_4",   "C_    683     "C_3H_6O",   "C_4H_10O",   "C_2H_4",   "C_2H_4O",            "C_2H_4S",
746     "SH_2",      "CH_4",       "CCLF_3",   "CC    684     "SH_2",      "CH_4",       "CCLF_3",   "CCl_2F_2",           "CHCl_2F",
747     "(CH_3)_2S", "N_2O",       "C_5H_10O", "C_    685     "(CH_3)_2S", "N_2O",       "C_5H_10O", "C_8H_6",             "(CH_2)_N",
748     "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8",   "C_    686     "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8",   "C_3H_6-Propylene",   "C_3H_6O",
749     "C_3H_6S",   "C_4H_4S",    "C_7H_8"           687     "C_3H_6S",   "C_4H_4S",    "C_7H_8"
750   };                                           << 688   } ;
751                                                   689 
752   static const G4float expStopping[numberOfMol << 690   static G4double expStopping[numberOfMolecula] = {
753      66.1f,  190.4f, 258.7f,  42.2f, 141.5f,   << 691      66.1,  190.4, 258.7,  42.2, 141.5,
754     210.9f,  279.6f, 198.8f,  31.0f, 267.5f,   << 692     210.9,  279.6, 198.8,  31.0, 267.5,
755     122.8f,  311.4f, 260.3f, 328.9f, 391.3f,   << 693     122.8,  311.4, 260.3, 328.9, 391.3,
756     206.6f,  374.0f, 422.0f, 432.0f, 398.0f,   << 694     206.6,  374.0, 422.0, 432.0, 398.0,
757     554.0f,  353.0f, 326.0f,  74.6f, 220.5f,   << 695     554.0,  353.0, 326.0,  74.6, 220.5,
758     197.4f,  362.0f, 170.0f, 330.5f, 211.3f,   << 696     197.4,  362.0, 170.0, 330.5, 211.3,
759     262.3f,  349.6f,  51.3f, 187.0f, 236.9f,   << 697     262.3,  349.6,  51.3, 187.0, 236.9,
760     121.9f,   35.8f, 247.0f, 292.6f, 268.0f,   << 698     121.9,   35.8, 247.0, 292.6, 268.0,
761     262.3f,   49.0f, 398.9f, 444.0f,  22.91f,  << 699     262.3,   49.0, 398.9, 444.0,  22.91,
762      68.0f,  155.0f,  84.0f,  74.2f, 254.7f,   << 700      68.0,  155.0,  84.0,  74.2, 254.7,
763     306.8f,  324.4f, 420.0f                    << 701     306.8,  324.4, 420.0
764   } ;                                             702   } ;
765                                                   703 
766   static const G4float expCharge[numberOfMolec << 704   static G4double expCharge[numberOfMolecula] = {
767     HeEff, HeEff, HeEff,  1.0f, HeEff,         << 705     HeEff, HeEff, HeEff,   1.0, HeEff,
768     HeEff, HeEff, HeEff,  1.0f,  1.0f,         << 706     HeEff, HeEff, HeEff,   1.0,   1.0,
769      1.0f, HeEff, HeEff, HeEff, HeEff,         << 707       1.0, HeEff, HeEff, HeEff, HeEff,
770     HeEff, HeEff, HeEff, HeEff, HeEff,            708     HeEff, HeEff, HeEff, HeEff, HeEff,
771     HeEff, HeEff, HeEff,  1.0f, HeEff,         << 709     HeEff, HeEff, HeEff,   1.0, HeEff,
772     HeEff, HeEff, HeEff, HeEff, HeEff,            710     HeEff, HeEff, HeEff, HeEff, HeEff,
773     HeEff, HeEff,  1.0f, HeEff, HeEff,         << 711     HeEff, HeEff,   1.0, HeEff, HeEff,
774     HeEff,  1.0f, HeEff, HeEff, HeEff,         << 712     HeEff,   1.0, HeEff, HeEff, HeEff,
775     HeEff,  1.0f, HeEff, HeEff,  1.0f,         << 713     HeEff,   1.0, HeEff, HeEff,   1.0,
776      1.0f,  1.0f,  1.0f,  1.0f, HeEff,         << 714       1.0,   1.0,   1.0,   1.0, HeEff,
777     HeEff, HeEff, HeEff                           715     HeEff, HeEff, HeEff
778   } ;                                             716   } ;
779                                                   717 
780   static const G4int numberOfAtomsPerMolecula[ << 718   static G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
781     3,  7, 10,  4,  6,                         << 719     3.0,  7.0, 10.0,  4.0,  6.0,
782     9, 12,  7,  4, 24,                         << 720     9.0, 12.0,  7.0,  4.0, 24.0,
783     12,14, 10, 13,  5,                         << 721     12.0, 14.0, 10.0, 13.0,  5.0,
784     5, 14, 18, 17, 17,                         << 722     5.0, 14.0, 18.0, 17.0, 17.0,
785     24,15, 13,  9,  8,                         << 723     24.0, 15.0, 13.0,  9.0,  8.0,
786     6, 14,  8,  8,  9,                         << 724     6.0, 14.0,  8.0,  8.0,  9.0,
787     10,15,  6,  7,  7,                         << 725     10.0, 15.0,  6.0,  7.0,  7.0,
788     3,  5,  5,  5,  5,                         << 726     3.0,  5.0,  5.0,  5.0,  5.0,
789     9,  3, 16, 14,  3,                         << 727     9.0,  3.0, 16.0, 14.0,  3.0,
790     9, 16, 11,  9, 10,                         << 728     9.0, 16.0, 11.0,  9.0, 10.0,
791     10, 9,  15};                               << 729     10.0,  9.0, 15.0
                                                   >> 730   } ;
792                                                   731 
793   // Search for the compaund in the table         732   // Search for the compaund in the table
794   for (std::size_t i=0; i<numberOfMolecula; ++ << 733   for (size_t i=0; i<numberOfMolecula; i++)
795     if(chFormula == nameOfMol[i]) {            << 734     {
796       expStopPower125 = ((G4double)expStopping << 735       if(chFormula == nameOfMol[i]) {
797         * (material->GetTotNbOfAtomsPerVolume( << 736         G4double exp125 = expStopping[i] *
798         ((G4double)(expCharge[i] * numberOfAto << 737                     (material->GetTotNbOfAtomsPerVolume()) /
799       return true;                             << 738                     (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
                                                   >> 739         SetExpStopPower125(exp125);
                                                   >> 740         return true;
                                                   >> 741       }
800     }                                             742     }
801   }                                            << 743   
802   return false;                                   744   return false;
803 }                                                 745 }
804                                                   746 
805 //....oooOO0OOooo........oooOO0OOooo........oo    747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806                                                   748 
807 G4double G4BraggModel::ChemicalFactor(G4double    749 G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy, 
808                                       G4double    750                                       G4double eloss125) const
809 {                                                 751 {
810   // Approximation of Chemical Factor accordin    752   // Approximation of Chemical Factor according to
811   // J.F.Ziegler and J.M.Manoyan, The stopping    753   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
812   // Nucl. Inst. & Meth. in Phys. Res. B35 (19    754   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
813                                                << 
814   static const G4double gamma25  = 1.0 + 25.0* << 
815   static const G4double gamma125 = 1.0 + 125.0 << 
816   static const G4double beta25   = std::sqrt(1 << 
817   static const G4double beta125  = std::sqrt(1 << 
818   static const G4double f12525   = 1.0 + G4Exp << 
819                                                << 
820   G4double gamma = 1.0 + kineticEnergy/proton_ << 
821   G4double beta  = std::sqrt(1.0 - 1.0/(gamma* << 
822                                                   755   
823   G4double factor = 1.0 + (expStopPower125/elo << 756   G4double gamma    = 1.0 + kineticEnergy/proton_mass_c2 ;    
824     (1.0 + G4Exp( 1.48 * ( beta/beta25    - 7. << 757   G4double gamma25  = 1.0 + 25.0*keV /proton_mass_c2 ;
                                                   >> 758   G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
                                                   >> 759   G4double beta     = sqrt(1.0 - 1.0/(gamma*gamma)) ;
                                                   >> 760   G4double beta25   = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
                                                   >> 761   G4double beta125  = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
                                                   >> 762   
                                                   >> 763   G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
                                                   >> 764                    (1.0 + exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
                                                   >> 765                    (1.0 + exp( 1.48 * ( beta/beta25    - 7.0 ) ) ) ;
825                                                   766 
826   return factor ;                                 767   return factor ;
827 }                                                 768 }
828                                                   769 
829 //....oooOO0OOooo........oooOO0OOooo........oo    770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
830                                                   771 
831                                                   772