Geant4 Cross Reference

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


  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: G4BraggIonModel.cc,v 1.30 2010/11/04 17:30:31 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-04 $
                                                   >>  28 //
 26 // -------------------------------------------     29 // -------------------------------------------------------------------
 27 //                                                 30 //
 28 // GEANT4 Class file                               31 // GEANT4 Class file
 29 //                                                 32 //
 30 //                                                 33 //
 31 // File name:   G4BraggIonModel                    34 // File name:   G4BraggIonModel
 32 //                                                 35 //
 33 // Author:        Vladimir Ivanchenko              36 // Author:        Vladimir Ivanchenko
 34 //                                                 37 //
 35 // Creation date: 13.10.2004                       38 // Creation date: 13.10.2004
 36 //                                                 39 //
 37 // Modifications:                                  40 // Modifications:
 38 // 11-05-05 Major optimisation of internal int     41 // 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko)
 39 // 29-11-05 Do not use G4Alpha class (V.Ivantc     42 // 29-11-05 Do not use G4Alpha class (V.Ivantchenko)
 40 // 15-02-06 ComputeCrossSectionPerElectron, Co     43 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 41 // 25-04-06 Add stopping data from ASTAR (V.Iv     44 // 25-04-06 Add stopping data from ASTAR (V.Ivanchenko)
 42 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV     45 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko)
 43 // 12-08-08 Added methods GetParticleCharge, G     46 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 
 44 //          CorrectionsAlongStep needed for io     47 //          CorrectionsAlongStep needed for ions(V.Ivanchenko)
 45 //                                                 48 //
 46                                                    49 
 47 // Class Description:                              50 // Class Description:
 48 //                                                 51 //
 49 // Implementation of energy loss and delta-ele     52 // Implementation of energy loss and delta-electron production by
 50 // slow charged heavy particles                    53 // slow charged heavy particles
 51                                                    54 
 52 // -------------------------------------------     55 // -------------------------------------------------------------------
 53 //                                                 56 //
 54                                                    57 
                                                   >>  58 
 55 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 56 //....oooOO0OOooo........oooOO0OOooo........oo     60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 57                                                    61 
 58 #include "G4BraggIonModel.hh"                      62 #include "G4BraggIonModel.hh"
 59 #include "G4PhysicalConstants.hh"              << 
 60 #include "G4SystemOfUnits.hh"                  << 
 61 #include "Randomize.hh"                            63 #include "Randomize.hh"
 62 #include "G4Electron.hh"                           64 #include "G4Electron.hh"
 63 #include "G4ParticleChangeForLoss.hh"              65 #include "G4ParticleChangeForLoss.hh"
                                                   >>  66 #include "G4LossTableManager.hh"
 64 #include "G4EmCorrections.hh"                      67 #include "G4EmCorrections.hh"
 65 #include "G4DeltaAngle.hh"                     << 
 66 #include "G4ICRU90StoppingData.hh"             << 
 67 #include "G4ASTARStopping.hh"                  << 
 68 #include "G4PSTARStopping.hh"                  << 
 69 #include "G4NistManager.hh"                    << 
 70 #include "G4Log.hh"                            << 
 71 #include "G4Exp.hh"                            << 
 72 #include "G4AutoLock.hh"                       << 
 73                                                    68 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 75                                                    70 
 76 G4ASTARStopping* G4BraggIonModel::fASTAR = nul <<  71 using namespace std;
 77                                                << 
 78 namespace                                      << 
 79 {                                              << 
 80   G4Mutex alphaMutex = G4MUTEX_INITIALIZER;    << 
 81 }                                              << 
 82                                                    72 
 83 G4BraggIonModel::G4BraggIonModel(const G4Parti     73 G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p,
 84                                  const G4Strin     74                                  const G4String& nam)
 85   : G4BraggModel(p, nam)                       <<  75   : G4VEmModel(nam),
 86 {                                              <<  76     corr(0),
 87   HeMass = 3.727417*CLHEP::GeV;                <<  77     particle(0),
 88   massFactor = 1000.*CLHEP::amu_c2/HeMass;     <<  78     fParticleChange(0),
                                                   >>  79     iMolecula(0),
                                                   >>  80     isIon(false),
                                                   >>  81     isInitialised(false)
                                                   >>  82 {
                                                   >>  83   SetHighEnergyLimit(2.0*MeV);
                                                   >>  84 
                                                   >>  85   HeMass           = 3.727417*GeV;
                                                   >>  86   rateMassHe2p     = HeMass/proton_mass_c2;
                                                   >>  87   lowestKinEnergy  = 1.0*keV/rateMassHe2p;
                                                   >>  88   massFactor       = 1000.*amu_c2/HeMass;
                                                   >>  89   theZieglerFactor = eV*cm2*1.0e-15;
                                                   >>  90   theElectron      = G4Electron::Electron();
                                                   >>  91   corrFactor       = 1.0;
                                                   >>  92   if(p) { SetParticle(p); }
                                                   >>  93   else  { SetParticle(theElectron); }
 89 }                                                  94 }
 90                                                    95 
 91 //....oooOO0OOooo........oooOO0OOooo........oo     96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 92                                                    97 
 93 G4BraggIonModel::~G4BraggIonModel()                98 G4BraggIonModel::~G4BraggIonModel()
                                                   >>  99 {}
                                                   >> 100 
                                                   >> 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 102 
                                                   >> 103 G4double G4BraggIonModel::MinEnergyCut(const G4ParticleDefinition*,
                                                   >> 104                                        const G4MaterialCutsCouple* couple)
 94 {                                                 105 {
 95   if(isFirstAlpha) {                           << 106   return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
 96     delete fASTAR;                             << 
 97     fASTAR = nullptr;                          << 
 98   }                                            << 
 99 }                                                 107 }
100                                                   108 
101 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102                                                   110 
103 void G4BraggIonModel::Initialise(const G4Parti    111 void G4BraggIonModel::Initialise(const G4ParticleDefinition* p,
104                                  const G4DataV << 112                                  const G4DataVector&)
105 {                                                 113 {
106   G4BraggModel::Initialise(p, ref);            << 114   if(p != particle) { SetParticle(p); }
107   const G4String& pname = particle->GetParticl << 115 
108   if(pname == "alpha") { isAlpha = true; }     << 116   corrFactor = chargeSquare;
109   if(isAlpha && fASTAR == nullptr) {           << 117 
110     G4AutoLock l(&alphaMutex);                 << 118   // always false before the run
111     if(fASTAR == nullptr) {                    << 119   SetDeexcitationFlag(false);
112       isFirstAlpha = true;                     << 120 
113       fASTAR = new G4ASTARStopping();          << 121   if(!isInitialised) {
114     }                                          << 122     isInitialised = true;
115     l.unlock();                                << 123 
116   }                                            << 124     G4String pname = particle->GetParticleName();
117   if(isFirstAlpha) {                           << 125     if(particle->GetParticleType() == "nucleus" &&
118     fASTAR->Initialise();                      << 126        pname != "deuteron" && pname != "triton") { isIon = true; }
                                                   >> 127 
                                                   >> 128     corr = G4LossTableManager::Instance()->EmCorrections();
                                                   >> 129 
                                                   >> 130     fParticleChange = GetParticleChangeForLoss();
119   }                                               131   }
120 }                                                 132 }
121                                                   133 
122                                                << 
123 //....oooOO0OOooo........oooOO0OOooo........oo    134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124                                                   135 
125 G4double G4BraggIonModel::GetChargeSquareRatio    136 G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
126                                                << 137                  const G4Material* mat,
127                                                << 138                  G4double kineticEnergy)
                                                   >> 139 {
                                                   >> 140   //G4cout << "G4BraggIonModel::GetChargeSquareRatio e= " <<  kineticEnergy << G4endl;
                                                   >> 141   // this method is called only for ions
                                                   >> 142   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
                                                   >> 143   corrFactor  = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy); 
                                                   >> 144   return corrFactor;
                                                   >> 145 }
                                                   >> 146 
                                                   >> 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 148 
                                                   >> 149 G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p,
                                                   >> 150               const G4Material* mat,
                                                   >> 151               G4double kineticEnergy)
                                                   >> 152 {
                                                   >> 153   //G4cout << "G4BraggIonModel::GetParticleCharge e= " <<  kineticEnergy << G4endl;
                                                   >> 154   // this method is called only for ions
                                                   >> 155   return corr->GetParticleCharge(p,mat,kineticEnergy);
                                                   >> 156 }
                                                   >> 157 
                                                   >> 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 159 
                                                   >> 160 G4double G4BraggIonModel::ComputeCrossSectionPerElectron(
                                                   >> 161                                            const G4ParticleDefinition* p,
                                                   >> 162                                                  G4double kineticEnergy,
                                                   >> 163                                                  G4double cutEnergy,
                                                   >> 164                                                  G4double maxKinEnergy)
128 {                                                 165 {
129   // this method is called only for ions, so n << 166   G4double cross     = 0.0;
130   if(isAlpha) { return 1.0; }                  << 167   G4double tmax      = MaxSecondaryEnergy(p, kineticEnergy);
131   return G4BraggModel::GetChargeSquareRatio(p, << 168   G4double maxEnergy = std::min(tmax,maxKinEnergy);
                                                   >> 169   if(cutEnergy < tmax) {
                                                   >> 170 
                                                   >> 171     G4double energy  = kineticEnergy + mass;
                                                   >> 172     G4double energy2 = energy*energy;
                                                   >> 173     G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
                                                   >> 174     cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
                                                   >> 175 
                                                   >> 176     cross *= twopi_mc2_rcl2*chargeSquare/beta2;
                                                   >> 177   }
                                                   >> 178  //   G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy 
                                                   >> 179  //          << " tmax= " << tmax << " cross= " << cross << G4endl;
                                                   >> 180  
                                                   >> 181   return cross;
132 }                                                 182 }
133                                                   183 
134 //....oooOO0OOooo........oooOO0OOooo........oo    184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135                                                   185 
136 G4double G4BraggIonModel::ComputeCrossSectionP    186 G4double G4BraggIonModel::ComputeCrossSectionPerAtom(
137                                            con    187                                            const G4ParticleDefinition* p,
138                                                << 188                                                  G4double kineticEnergy,
139                                                << 189              G4double Z, G4double,
140                                                   190                                                  G4double cutEnergy,
141                                                   191                                                  G4double maxEnergy)
142 {                                                 192 {
143   G4double sigma =                             << 193   G4double cross = Z*ComputeCrossSectionPerElectron
144     Z*ComputeCrossSectionPerElectron(p,kinEner << 194                                          (p,kineticEnergy,cutEnergy,maxEnergy);
145   if(isAlpha) {                                << 195   return cross;
146     sigma *= (HeEffChargeSquare(Z, kinEnergy/C << 
147   }                                            << 
148   return sigma;                                << 
149 }                                                 196 }
150                                                   197 
151 //....oooOO0OOooo........oooOO0OOooo........oo    198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152                                                   199 
153 G4double G4BraggIonModel::CrossSectionPerVolum    200 G4double G4BraggIonModel::CrossSectionPerVolume(
154                                            con << 201              const G4Material* material,
155                                            con    202                                            const G4ParticleDefinition* p,
156                                                << 203                                                  G4double kineticEnergy,
157                                                   204                                                  G4double cutEnergy,
158                                                   205                                                  G4double maxEnergy)
159 {                                                 206 {
160   G4double sigma = material->GetElectronDensit << 207   G4double eDensity = material->GetElectronDensity();
161     ComputeCrossSectionPerElectron(p,kinEnergy << 208   G4double cross = eDensity*ComputeCrossSectionPerElectron
162   if(isAlpha) {                                << 209                                          (p,kineticEnergy,cutEnergy,maxEnergy);
163     const G4double zeff = material->GetTotNbOf << 210   return cross;
164       material->GetTotNbOfAtomsPerVolume();    << 
165     sigma *= (HeEffChargeSquare(zeff, kinEnerg << 
166   }                                            << 
167   return sigma;                                << 
168 }                                                 211 }
169                                                   212 
170 //....oooOO0OOooo........oooOO0OOooo........oo    213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171                                                   214 
172 G4double G4BraggIonModel::ComputeDEDXPerVolume    215 G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material,
173                                                << 216                  const G4ParticleDefinition* p,
174                                                << 217                  G4double kineticEnergy,
175                                                << 218                  G4double cutEnergy)
176 {                                              << 219 {
177   const G4double tmax = MaxSecondaryEnergy(p,  << 220   G4double tmax  = MaxSecondaryEnergy(p, kineticEnergy);
178   const G4double tlim = lowestKinEnergy*massRa << 221   G4double tmin  = min(cutEnergy, tmax);
179   const G4double tmin = std::max(std::min(cut, << 222   G4double tkin  = kineticEnergy/massRate;
180   G4double dedx = 0.0;                         << 223   G4double dedx  = 0.0;
                                                   >> 224   if(tkin > lowestKinEnergy) dedx = DEDX(material, tkin);
                                                   >> 225   else      dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
181                                                   226 
182   if(kineticEnergy < tlim) {                   << 227   if (cutEnergy < tmax) {
183     dedx = HeDEDX(material, tlim)*std::sqrt(ki << 
184   } else {                                     << 
185     dedx = HeDEDX(material, kineticEnergy);    << 
186                                                   228 
187     if (tmin < tmax) {                         << 229     G4double tau   = kineticEnergy/mass;
188       const G4double tau = kineticEnergy/mass; << 230     G4double gam   = tau + 1.0;
189       const G4double x   = tmin/tmax;          << 231     G4double bg2   = tau * (tau+2.0);
190                                                << 232     G4double beta2 = bg2/(gam*gam);
191       G4double del =                           << 233     G4double x     = tmin/tmax;
192         (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * << 234 
193   CLHEP::twopi_mc2_rcl2*material->GetElectronD << 235     dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
194       if(isAlpha) {                            << 236           * (material->GetElectronDensity())/beta2;
195   const G4double zeff = material->GetTotNbOfEl << 
196     material->GetTotNbOfAtomsPerVolume();      << 
197   heChargeSquare = HeEffChargeSquare(zeff, kin << 
198   del *= heChargeSquare;                       << 
199       }                                        << 
200       dedx += del;                             << 
201     }                                          << 
202   }                                               237   }
203   dedx = std::max(dedx, 0.0);                  << 238 
204   /*                                           << 239   // now compute the total ionization loss
205     G4cout << "BraggIon: " << material->GetNam << 240 
206            << " E(MeV)=" << kineticEnergy/MeV  << 241   if (dedx < 0.0) dedx = 0.0 ;
207            << " Tmin(MeV)=" << tmin << " dedx( << 242 
208            << dedx*gram/(MeV*cm2*material->Get << 243   dedx *= chargeSquare;
209            << " q2=" << chargeSquare << G4endl << 244 
210   */                                           << 245   //G4cout << " tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = " 
                                                   >> 246   //       << dedx*gram/(MeV*cm2*material->GetDensity()) 
                                                   >> 247   //       << " q2 = " << chargeSquare <<  G4endl;
                                                   >> 248 
211   return dedx;                                    249   return dedx;
212 }                                                 250 }
213                                                   251 
214 //....oooOO0OOooo........oooOO0OOooo........oo    252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215                                                   253 
216 void G4BraggIonModel::CorrectionsAlongStep(con    254 void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
217                                            con << 255              const G4DynamicParticle* dp,
218                                            con << 256              G4double& eloss,
219                                            G4d << 257              G4double&,
220 {                                              << 258              G4double /*length*/)
221   // no correction for alpha                   << 259 {
222   if(isAlpha) { return; }                      << 260   // this method is called only for ions
223                                                << 
224   // no correction at a small step at the last << 
225   const G4double preKinEnergy = dp->GetKinetic << 
226   if(eloss >= preKinEnergy || eloss < preKinEn << 
227                                                << 
228   // corrections only for ions                 << 
229   const G4ParticleDefinition* p = dp->GetDefin    261   const G4ParticleDefinition* p = dp->GetDefinition();
230   if(p != particle) { SetParticle(p); }        << 
231                                                << 
232   // effective energy and charge at a step     << 
233   const G4Material* mat = couple->GetMaterial(    262   const G4Material* mat = couple->GetMaterial();
234   const G4double e = std::max(preKinEnergy - e << 263   G4double preKinEnergy = dp->GetKineticEnergy();
235   const G4double q20 = corr->EffectiveChargeSq << 264   G4double e = preKinEnergy - eloss*0.5;
236   const G4double q2 = corr->EffectiveChargeSqu << 265   if(e < 0.0) e = preKinEnergy*0.5;
237   const G4double qfactor = q2/q20;             << 266 
238   /*                                           << 267   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
239     G4cout << "G4BraggIonModel::CorrectionsAlo << 268   GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
240     << preKinEnergy << " Eeff(MeV)=" << e      << 269   G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor; 
241     << " eloss=" << eloss << " elossnew=" << e << 270   eloss *= qfactor; 
242     << " qfactor=" << qfactor << " Qpre=" << q << 271 
243     << p->GetParticleName() <<G4endl;          << 272   //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " <<  e 
244   */                                           << 273   //   << " qfactor= " << qfactor << " " << p->GetParticleName() <<G4endl;
245   eloss *= qfactor;                            << 
246 }                                                 274 }
247                                                   275 
248 //....oooOO0OOooo........oooOO0OOooo........oo    276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249                                                   277 
250 G4int G4BraggIonModel::HasMaterialForHe(const  << 278 void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
                                                   >> 279           const G4MaterialCutsCouple*,
                                                   >> 280           const G4DynamicParticle* dp,
                                                   >> 281           G4double xmin,
                                                   >> 282           G4double maxEnergy)
251 {                                                 283 {
252   const G4String& chFormula = mat->GetChemical << 284   G4double tmax = MaxSecondaryKinEnergy(dp);
253   if(chFormula.empty()) { return -1; }         << 285   G4double xmax = std::min(tmax, maxEnergy);
                                                   >> 286   if(xmin >= xmax) return;
                                                   >> 287 
                                                   >> 288   G4double kineticEnergy = dp->GetKineticEnergy();
                                                   >> 289   G4double energy  = kineticEnergy + mass;
                                                   >> 290   G4double energy2 = energy*energy;
                                                   >> 291   G4double beta2   = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
                                                   >> 292   G4double grej    = 1.0;
                                                   >> 293   G4double deltaKinEnergy, f;
                                                   >> 294 
                                                   >> 295   G4ThreeVector direction = dp->GetMomentumDirection();
                                                   >> 296 
                                                   >> 297   // sampling follows ...
                                                   >> 298   do {
                                                   >> 299     G4double q = G4UniformRand();
                                                   >> 300     deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
                                                   >> 301 
                                                   >> 302     f = 1.0 - beta2*deltaKinEnergy/tmax;
                                                   >> 303 
                                                   >> 304     if(f > grej) {
                                                   >> 305         G4cout << "G4BraggIonModel::SampleSecondary Warning! "
                                                   >> 306                << "Majorant " << grej << " < "
                                                   >> 307                << f << " for e= " << deltaKinEnergy
                                                   >> 308                << G4endl;
                                                   >> 309     }
                                                   >> 310 
                                                   >> 311   } while( grej*G4UniformRand() >= f );
                                                   >> 312 
                                                   >> 313   G4double deltaMomentum =
                                                   >> 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));
                                                   >> 320 
                                                   >> 321   G4double phi = twopi * G4UniformRand() ;
                                                   >> 322 
                                                   >> 323   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
                                                   >> 324   deltaDirection.rotateUz(direction);
                                                   >> 325 
                                                   >> 326   // create G4DynamicParticle object for delta ray
                                                   >> 327   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
                                                   >> 328                deltaKinEnergy);
                                                   >> 329 
                                                   >> 330   vdp->push_back(delta);
                                                   >> 331 
                                                   >> 332   // Change kinematics of primary particle
                                                   >> 333   kineticEnergy       -= deltaKinEnergy;
                                                   >> 334   G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
                                                   >> 335   finalP               = finalP.unit();
                                                   >> 336 
                                                   >> 337   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
                                                   >> 338   fParticleChange->SetProposedMomentumDirection(finalP);
                                                   >> 339 }
                                                   >> 340 
                                                   >> 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 342 
                                                   >> 343 G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
                                                   >> 344                G4double kinEnergy)
                                                   >> 345 {
                                                   >> 346   if(pd != particle) SetParticle(pd);
                                                   >> 347   G4double tau  = kinEnergy/mass;
                                                   >> 348   G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
                                                   >> 349                   (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
                                                   >> 350   return tmax;
                                                   >> 351 }
                                                   >> 352 
                                                   >> 353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 354 
                                                   >> 355 G4bool G4BraggIonModel::HasMaterial(const G4Material* material)
                                                   >> 356 {
                                                   >> 357   const size_t numberOfMolecula = 11 ;
                                                   >> 358   SetMoleculaNumber(numberOfMolecula) ;
                                                   >> 359   G4String chFormula = material->GetChemicalFormula() ;
254                                                   360 
255   // ICRU Report N49, 1993. Ziegler model for     361   // ICRU Report N49, 1993. Ziegler model for He.
256                                                << 362   static G4String molName[numberOfMolecula] = {
257   static const G4int numberOfMolecula = 11;    << 
258   static const G4String molName[numberOfMolecu << 
259     "CaF_2",  "Cellulose_Nitrate",  "LiF", "Po    363     "CaF_2",  "Cellulose_Nitrate",  "LiF", "Policarbonate",  
260     "(C_2H_4)_N-Polyethylene",  "(C_2H_4)_N-Po    364     "(C_2H_4)_N-Polyethylene",  "(C_2H_4)_N-Polymethly_Methacralate",
261     "Polysterene", "SiO_2", "NaI", "H_2O",        365     "Polysterene", "SiO_2", "NaI", "H_2O",
262     "Graphite" };                              << 366     "Graphite" } ;
263                                                   367 
264   // Search for the material in the table         368   // Search for the material in the table
265   for (G4int i=0; i<numberOfMolecula; ++i) {   << 369   for (size_t i=0; i<numberOfMolecula; i++) {
266     if (chFormula == molName[i]) {             << 370       if (chFormula == molName[i]) {
267       return i;                                << 371         SetMoleculaNumber(i) ;
268     }                                          << 372   return true ;
                                                   >> 373       }
269   }                                               374   }
270   return -1;                                   << 375   return false ;
271 }                                                 376 }
272                                                   377 
273 //....oooOO0OOooo........oooOO0OOooo........oo    378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274                                                   379 
275 G4double G4BraggIonModel::HeStoppingPower(cons << 380 G4double G4BraggIonModel::StoppingPower(const G4Material* material,
                                                   >> 381           G4double kineticEnergy) 
276 {                                                 382 {
277   G4double ionloss = 0.0;                      << 383   G4double ionloss = 0.0 ;
278   if (iMolecula >= 0) {                        << 384 
                                                   >> 385   if (iMolecula < 11) {
279                                                   386   
280     // The data and the fit from:                 387     // The data and the fit from: 
281     // ICRU Report N49, 1993. Ziegler's model     388     // ICRU Report N49, 1993. Ziegler's model for alpha
282     // He energy in internal units of parametr    389     // He energy in internal units of parametrisation formula (MeV)
283     // Input scaled energy of a proton or Gene << 
284     G4double T = kineticEnergy/(massRate*CLHEP << 
285                                                   390 
286     static const G4float a[11][5] = {          << 391     G4double T = kineticEnergy*rateMassHe2p/MeV ;
287        {9.43672f, 0.54398f, 84.341f,  1.3705f, << 392 
288        {67.1503f, 0.41409f, 404.512f, 148.97f, << 393     static G4double a[11][5] = {
289        {5.11203f, 0.453f,   36.718f,  50.6f,   << 394        {9.43672, 0.54398, 84.341, 1.3705, 57.422},
290        {61.793f,  0.48445f, 361.537f, 57.889f, << 395        {67.1503, 0.41409, 404.512, 148.97, 20.99},
291        {7.83464f, 0.49804f, 160.452f, 3.192f,  << 396        {5.11203, 0.453,  36.718,  50.6,  28.058}, 
292        {19.729f,  0.52153f, 162.341f, 58.35f,  << 397        {61.793, 0.48445, 361.537, 57.889, 50.674},
293        {26.4648f, 0.50112f, 188.913f, 30.079f, << 398        {7.83464, 0.49804, 160.452, 3.192, 0.71922},
294        {7.8655f,  0.5205f,  63.96f,   51.32f,  << 399        {19.729, 0.52153, 162.341, 58.35, 25.668}, 
295        {8.8965f,  0.5148f,  339.36f,  1.7205f, << 400        {26.4648, 0.50112, 188.913, 30.079, 16.509},
296        {2.959f,   0.53255f, 34.247f,  60.655f, << 401        {7.8655, 0.5205, 63.96, 51.32, 67.775},
297        {3.80133f, 0.41590f, 12.9966f, 117.83f, << 402        {8.8965, 0.5148, 339.36, 1.7205, 0.70423},
298                                                << 403        {2.959, 0.53255, 34.247, 60.655, 15.153}, 
299     static const G4double atomicWeight[11] = { << 404        {3.80133, 0.41590, 12.9966, 117.83, 242.28} };   
300        101.96128f, 44.0098f, 16.0426f, 28.0536 << 405 
301        104.1512f,  44.665f,  60.0843f, 18.0152 << 406     static G4double atomicWeight[11] = {
302                                                << 407        101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
303     const G4int i = iMolecula;                 << 408        104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};       
304                                                << 409 
305     G4double slow = (G4double)(a[i][0]);       << 410     G4int i = iMolecula;
306                                                << 
307     G4double x1 = (G4double)(a[i][1]);         << 
308     G4double x2 = (G4double)(a[i][2]);         << 
309     G4double x3 = (G4double)(a[i][3]);         << 
310     G4double x4 = (G4double)(a[i][4]);         << 
311                                                   411 
312     // Free electron gas model                    412     // Free electron gas model
313     if ( T < 0.001 ) {                            413     if ( T < 0.001 ) {
314       G4double shigh = G4Log( 1.0 + x3*1000.0  << 414       G4double slow  = a[i][0] ;
                                                   >> 415       G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
                                                   >> 416    * a[i][2]*1000.0 ;
315       ionloss  = slow*shigh / (slow + shigh) ;    417       ionloss  = slow*shigh / (slow + shigh) ;
316       ionloss *= std::sqrt(T*1000.0) ;         << 418       ionloss *= sqrt(T*1000.0) ;
317                                                   419 
318       // Main parametrisation                     420       // Main parametrisation
319     } else {                                      421     } else {
320       slow  *= G4Exp(G4Log(T*1000.0)*x1) ;     << 422       G4double slow  = a[i][0] * pow((T*1000.0), a[i][1]) ;
321       G4double shigh = G4Log( 1.0 + x3/T + x4* << 423       G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
322       ionloss = slow*shigh / (slow + shigh) ;     424       ionloss = slow*shigh / (slow + shigh) ;
323        /*                                         425        /*
324          G4cout << "## " << i << ". T= " << T  << 426    G4cout << "## " << i << ". T= " << T << " slow= " << slow
325          << " a0= " << a[i][0] << " a1= " << a << 427    << " a0= " << a[i][0] << " a1= " << a[i][1] 
326          << " shigh= " << shigh                << 428    << " shigh= " << shigh 
327          << " dedx= " << ionloss << " q^2= " < << 429    << " dedx= " << ionloss << " q^2= " <<  HeEffChargeSquare(z, T*MeV)
328          << G4endl;                            << 430    << G4endl;
329        */                                         431        */
330     }                                             432     }
331     ionloss = std::max(ionloss, 0.0) * atomicW << 433     if ( ionloss < 0.0) ionloss = 0.0 ;
                                                   >> 434 
                                                   >> 435     // He effective charge
                                                   >> 436     G4double aa = atomicWeight[iMolecula];
                                                   >> 437     ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
                                                   >> 438 
                                                   >> 439   // pure material (normally not the case for this function)
                                                   >> 440   } else if(1 == (material->GetNumberOfElements())) {
                                                   >> 441     G4double z = material->GetZ() ;
                                                   >> 442     ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;  
332   }                                               443   }
                                                   >> 444   
333   return ionloss;                                 445   return ionloss;
334 }                                                 446 }
335                                                   447 
336 //....oooOO0OOooo........oooOO0OOooo........oo    448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337                                                   449 
338 G4double G4BraggIonModel::HeElectronicStopping << 450 G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
339                           const G4double kinet << 451                                                   G4double kineticEnergy) const
340 {                                                 452 {
341   G4double ionloss ;                              453   G4double ionloss ;
342   G4int i = std::min(z-1, 91);  // index of at << 454   G4int i = G4int(z)-1 ;  // index of atom
343   //G4cout << "ElectronicStoppingPower z=" <<  << 455   if(i < 0)  i = 0 ;
344   // << " E=" << kineticEnergy << G4endl;      << 456   if(i > 91) i = 91 ;
                                                   >> 457 
345   // The data and the fit from:                   458   // The data and the fit from:
346   // ICRU Report 49, 1993. Ziegler's type of p    459   // ICRU Report 49, 1993. Ziegler's type of parametrisations.
347   // Proton kinetic energy for parametrisation    460   // Proton kinetic energy for parametrisation (keV/amu)
348   // He energy in internal units of parametris << 461 
349   //G4double T = kineticEnergy*rateMassHe2p/CL << 462    // He energy in internal units of parametrisation formula (MeV)
350   G4double T = kineticEnergy/CLHEP::MeV;       << 463   G4double T = kineticEnergy*rateMassHe2p/MeV ;
351                                                << 464 
352   static const G4float a[92][5] = {            << 465   static G4double a[92][5] = {
353     {  0.35485f, 0.6456f, 6.01525f,  20.8933f, << 466     {0.35485, 0.6456, 6.01525,  20.8933, 4.3515
354    },{ 0.58f,    0.59f,   6.3f,      130.0f,   << 467    },{ 0.58,    0.59,   6.3,     130.0,   44.07
355    },{ 1.42f,    0.49f,   12.25f,    32.0f,    << 468    },{ 1.42,    0.49,   12.25,    32.0,    9.161
356    },{ 2.206f,   0.51f,   15.32f,    0.25f,    << 469    },{ 2.206,   0.51,   15.32,    0.25,    8.995 //Be Ziegler77
357        // },{ 2.1895f,  0.47183,7.2362f,   134 << 470        // },{ 2.1895,  0.47183,7.2362,   134.30,  197.96 //Be from ICRU
358    },{ 3.691f,   0.4128f, 18.48f,    50.72f,   << 471    },{ 3.691,   0.4128, 18.48,    50.72,   9.0
359    },{ 3.83523f, 0.42993f,12.6125f,  227.41f,  << 472    },{ 3.83523, 0.42993,12.6125,  227.41,  188.97
360        // },{ 1.9259f,  0.5550f, 27.15125f, 26 << 473    },{ 1.9259,  0.5550, 27.15125, 26.0665, 6.2768
361    },{ 1.9259f,  0.5550f, 27.1513f,  26.0665f, << 474    },{ 2.81015, 0.4759, 50.0253,  10.556,  1.0382
362    },{ 2.81015f, 0.4759f, 50.0253f,  10.556f,  << 475    },{ 1.533,   0.531,  40.44,    18.41,   2.718
363    },{ 1.533f,   0.531f,  40.44f,    18.41f,   << 476    },{ 2.303,   0.4861, 37.01,    37.96,   5.092
364    },{ 2.303f,   0.4861f, 37.01f,    37.96f,   << 
365        // Z= 11-20                                477        // Z= 11-20
366    },{ 9.894f,   0.3081f, 23.65f,    0.384f,   << 478    },{ 9.894,   0.3081, 23.65,    0.384,   92.93
367    },{ 4.3f,     0.47f,   34.3f,     3.3f,     << 479    },{ 4.3,     0.47,   34.3,     3.3,     12.74
368    },{ 2.5f,     0.625f,  45.7f,     0.1f,     << 480    },{ 2.5,     0.625,  45.7,     0.1,     4.359
369    },{ 2.1f,     0.65f,   49.34f,    1.788f,   << 481    },{ 2.1,     0.65,   49.34,    1.788,   4.133
370    },{ 1.729f,   0.6562f, 53.41f,    2.405f,   << 482    },{ 1.729,   0.6562, 53.41,    2.405,   3.845
371    },{ 1.402f,   0.6791f, 58.98f,    3.528f,   << 483    },{ 1.402,   0.6791, 58.98,    3.528,   3.211
372    },{ 1.117f,   0.7044f, 69.69f,    3.705f,   << 484    },{ 1.117,   0.7044, 69.69,    3.705,    2.156
373    },{ 2.291f,   0.6284f, 73.88f,    4.478f,   << 485    },{ 2.291,   0.6284, 73.88,    4.478,    2.066
374    },{ 8.554f,   0.3817f, 83.61f,    11.84f,   << 486    },{ 8.554,   0.3817, 83.61,    11.84,    1.875
375    },{ 6.297f,   0.4622f, 65.39f,    10.14f,   << 487    },{ 6.297,   0.4622, 65.39,    10.14,    5.036
376        // Z= 21-30                                488        // Z= 21-30     
377    },{ 5.307f,   0.4918f, 61.74f,    12.4f,    << 489    },{ 5.307,   0.4918, 61.74,    12.4,    6.665
378    },{ 4.71f,    0.5087f, 65.28f,    8.806f,   << 490    },{ 4.71,    0.5087, 65.28,    8.806,    5.948
379    },{ 6.151f,   0.4524f, 83.0f,     18.31f,   << 491    },{ 6.151,   0.4524, 83.0,    18.31,    2.71
380    },{ 6.57f,    0.4322f, 84.76f,    15.53f,   << 492    },{ 6.57,    0.4322, 84.76,    15.53,    2.779
381    },{ 5.738f,   0.4492f, 84.6f,     14.18f,   << 493    },{ 5.738,   0.4492, 84.6,    14.18,    3.101
382    },{ 5.013f,   0.4707f, 85.8f,     16.55f,   << 494    },{ 5.013,   0.4707, 85.8,    16.55,    3.211
383    },{ 4.32f,    0.4947f, 76.14f,    10.85f,   << 495    },{ 4.32,    0.4947, 76.14,    10.85,    5.441
384    },{ 4.652f,   0.4571f, 80.73f,    22.0f,    << 496    },{ 4.652,   0.4571, 80.73,    22.0,    4.952
385    },{ 3.114f,   0.5236f, 76.67f,    7.62f,    << 497    },{ 3.114,   0.5236, 76.67,    7.62,    6.385
386    },{ 3.114f,   0.5236f, 76.67f,    7.62f,    << 498    },{ 3.114,   0.5236, 76.67,    7.62,    7.502
387        // Z= 31-40                                499        // Z= 31-40
388    },{ 3.114f,   0.5236f, 76.67f,    7.62f,    << 500    },{ 3.114,   0.5236, 76.67,    7.62,    8.514
389    },{ 5.746f,   0.4662f, 79.24f,    1.185f,   << 501    },{ 5.746,   0.4662, 79.24,    1.185,    7.993
390    },{ 2.792f,   0.6346f, 106.1f,    0.2986f,  << 502    },{ 2.792,   0.6346, 106.1,    0.2986,   2.331
391    },{ 4.667f,   0.5095f, 124.3f,    2.102f,   << 503    },{ 4.667,   0.5095, 124.3,    2.102,    1.667
392    },{ 2.44f,    0.6346f, 105.0f,    0.83f,    << 504    },{ 2.44,    0.6346, 105.0,    0.83,    2.851
393    },{ 1.413f,   0.7377f, 147.9f,    1.466f,   << 505    },{ 1.413,   0.7377, 147.9,    1.466,    1.016
394    },{ 11.72f,   0.3826f, 102.8f,    9.231f,   << 506    },{ 11.72,   0.3826, 102.8,    9.231,    4.371
395    },{ 7.126f,   0.4804f, 119.3f,    5.784f,   << 507    },{ 7.126,   0.4804, 119.3,    5.784,    2.454
396    },{ 11.61f,   0.3955f, 146.7f,    7.031f,   << 508    },{ 11.61,   0.3955, 146.7,    7.031,    1.423
397    },{ 10.99f,   0.41f,   163.9f,    7.1f,     << 509    },{ 10.99,   0.41,   163.9,   7.1,      1.052
398        // Z= 41-50                                510        // Z= 41-50
399    },{ 9.241f,   0.4275f, 163.1f,    7.954f,   << 511    },{ 9.241,   0.4275, 163.1,    7.954,    1.102
400    },{ 9.276f,   0.418f,  157.1f,    8.038f,   << 512    },{ 9.276,   0.418,  157.1,   8.038,    1.29
401    },{ 3.999f,   0.6152f, 97.6f,     1.297f,   << 513    },{ 3.999,   0.6152, 97.6,    1.297,    5.792
402    },{ 4.306f,   0.5658f, 97.99f,    5.514f,   << 514    },{ 4.306,   0.5658, 97.99,    5.514,    5.754
403    },{ 3.615f,   0.6197f, 86.26f,    0.333f,   << 515    },{ 3.615,   0.6197, 86.26,    0.333,    8.689
404    },{ 5.8f,     0.49f,   147.2f,    6.903f,   << 516    },{ 5.8,     0.49,   147.2,   6.903,    1.289
405    },{ 5.6f,     0.49f,   130.0f,    10.0f,    << 517    },{ 5.6,     0.49,   130.0,   10.0,     2.844
406    },{ 3.55f,    0.6068f, 124.7f,    1.112f,   << 518    },{ 3.55,    0.6068, 124.7,    1.112,    3.119
407    },{ 3.6f,     0.62f,   105.8f,    0.1692f,  << 519    },{ 3.6,     0.62,   105.8,   0.1692,   6.026
408    },{ 5.4f,     0.53f,   103.1f,    3.931f,   << 520    },{ 5.4,     0.53,   103.1,   3.931,    7.767
409        // Z= 51-60                                521        // Z= 51-60
410    },{ 3.97f,    0.6459f, 131.8f,    0.2233f,  << 522    },{ 3.97,    0.6459, 131.8,    0.2233,   2.723
411    },{ 3.65f,    0.64f,   126.8f,    0.6834f,  << 523    },{ 3.65,    0.64,   126.8,   0.6834,   3.411
412    },{ 3.118f,   0.6519f, 164.9f,    1.208f,   << 524    },{ 3.118,   0.6519, 164.9,    1.208,    1.51
413    },{ 3.949f,   0.6209f, 200.5f,    1.878f,   << 525    },{ 3.949,   0.6209, 200.5,    1.878,    0.9126
414    },{ 14.4f,    0.3923f, 152.5f,    8.354f,   << 526    },{ 14.4,    0.3923, 152.5,    8.354,    2.597
415    },{ 10.99f,   0.4599f, 138.4f,    4.811f,   << 527    },{ 10.99,   0.4599, 138.4,    4.811,    3.726
416    },{ 16.6f,    0.3773f, 224.1f,    6.28f,    << 528    },{ 16.6,    0.3773, 224.1,    6.28,    0.9121
417    },{ 10.54f,   0.4533f, 159.3f,    4.832f,   << 529    },{ 10.54,   0.4533, 159.3,   4.832,    2.529
418    },{ 10.33f,   0.4502f, 162.0f,    5.132f,   << 530    },{ 10.33,   0.4502, 162.0,   5.132,    2.444
419    },{ 10.15f,   0.4471f, 165.6f,    5.378f,   << 531    },{ 10.15,   0.4471, 165.6,   5.378,    2.328
420        // Z= 61-70                                532        // Z= 61-70
421    },{ 9.976f,   0.4439f, 168.0f,    5.721f,   << 533    },{ 9.976,   0.4439, 168.0,   5.721,    2.258
422    },{ 9.804f,   0.4408f, 176.2f,    5.675f,   << 534    },{ 9.804,   0.4408, 176.2,   5.675,    1.997
423    },{ 14.22f,   0.363f,  228.4f,    7.024f,   << 535    },{ 14.22,   0.363,  228.4,   7.024,    1.016
424    },{ 9.952f,   0.4318f, 233.5f,    5.065f,   << 536    },{ 9.952,   0.4318, 233.5,   5.065,    0.9244
425    },{ 9.272f,   0.4345f, 210.0f,    4.911f,   << 537    },{ 9.272,   0.4345, 210.0,   4.911,    1.258
426    },{ 10.13f,   0.4146f, 225.7f,    5.525f,   << 538    },{ 10.13,   0.4146, 225.7,   5.525,    1.055
427    },{ 8.949f,   0.4304f, 213.3f,    5.071f,   << 539    },{ 8.949,   0.4304, 213.3,   5.071,    1.221
428    },{ 11.94f,   0.3783f, 247.2f,    6.655f,   << 540    },{ 11.94,   0.3783, 247.2,   6.655,    0.849
429    },{ 8.472f,   0.4405f, 195.5f,    4.051f,   << 541    },{ 8.472,   0.4405, 195.5,   4.051,    1.604
430    },{ 8.301f,   0.4399f, 203.7f,    3.667f,   << 542    },{ 8.301,   0.4399, 203.7,   3.667,    1.459
431        // Z= 71-80                                543        // Z= 71-80
432    },{ 6.567f,   0.4858f, 193.0f,    2.65f,    << 544    },{ 6.567,   0.4858, 193.0,   2.65,     1.66
433    },{ 5.951f,   0.5016f, 196.1f,    2.662f,   << 545    },{ 5.951,   0.5016, 196.1,   2.662,    1.589
434    },{ 7.495f,   0.4523f, 251.4f,    3.433f,   << 546    },{ 7.495,   0.4523, 251.4,   3.433,    0.8619
435    },{ 6.335f,   0.4825f, 255.1f,    2.834f,   << 547    },{ 6.335,   0.4825, 255.1,   2.834,    0.8228
436    },{ 4.314f,   0.5558f, 214.8f,    2.354f,   << 548    },{ 4.314,   0.5558, 214.8,   2.354,    1.263
437    },{ 4.02f,    0.5681f, 219.9f,    2.402f,   << 549    },{ 4.02,    0.5681, 219.9,   2.402,    1.191
438    },{ 3.836f,   0.5765f, 210.2f,    2.742f,   << 550    },{ 3.836,   0.5765, 210.2,   2.742,    1.305
439    },{ 4.68f,    0.5247f, 244.7f,    2.749f,   << 551    },{ 4.68,    0.5247, 244.7,   2.749,    0.8962
440    },{ 2.892f,   0.6204f, 208.6f,    2.415f,   << 552    },{ 2.892,   0.6204, 208.6,   2.415,    1.416 //Au Z77
441        // },{ 3.223f,   0.5883f, 232.7f,   2.9 << 553        // },{ 3.223,   0.5883, 232.7,   2.954,    1.05  //Au ICRU
442    },{ 2.892f,   0.6204f, 208.6f,    2.415f,   << 554    },{ 2.892,   0.6204, 208.6,   2.415,    1.416
443        // Z= 81-90                                555        // Z= 81-90
444    },{ 4.728f,   0.5522f, 217.0f,    3.091f,   << 556    },{ 4.728,   0.5522, 217.0,   3.091,    1.386
445    },{ 6.18f,    0.52f,   170.0f,    4.0f,     << 557    },{ 6.18,    0.52,   170.0,   4.0,      3.224
446    },{ 9.0f,     0.47f,   198.0f,    3.8f,     << 558    },{ 9.0,     0.47,   198.0,   3.8,      2.032
447    },{ 2.324f,   0.6997f, 216.0f,    1.599f,   << 559    },{ 2.324,   0.6997, 216.0,   1.599,    1.399
448    },{ 1.961f,   0.7286f, 223.0f,    1.621f,   << 560    },{ 1.961,   0.7286, 223.0,   1.621,    1.296
449    },{ 1.75f,    0.7427f, 350.1f,    0.9789f,  << 561    },{ 1.75,    0.7427, 350.1,   0.9789,   0.5507
450    },{ 10.31f,   0.4613f, 261.2f,    4.738f,   << 562    },{ 10.31,   0.4613, 261.2,   4.738,    0.9899
451    },{ 7.962f,   0.519f,  235.7f,    4.347f,   << 563    },{ 7.962,   0.519,  235.7,   4.347,    1.313
452    },{ 6.227f,   0.5645f, 231.9f,    3.961f,   << 564    },{ 6.227,   0.5645, 231.9,   3.961,    1.379
453    },{ 5.246f,   0.5947f, 228.6f,    4.027f,   << 565    },{ 5.246,   0.5947, 228.6,   4.027,    1.432
454        // Z= 91-92                                566        // Z= 91-92
455    },{ 5.408f,   0.5811f, 235.7f,    3.961f,   << 567    },{ 5.408,   0.5811, 235.7,   3.961,    1.358
456    },{ 5.218f,   0.5828f, 245.0f,    3.838f,   << 568    },{ 5.218,   0.5828, 245.0,   3.838,    1.25}
457   };                                              569   };
458                                                   570 
459   G4double slow = (G4double)(a[i][0]);         << 
460                                                << 
461   G4double x1 = (G4double)(a[i][1]);           << 
462   G4double x2 = (G4double)(a[i][2]);           << 
463   G4double x3 = (G4double)(a[i][3]);           << 
464   G4double x4 = (G4double)(a[i][4]);           << 
465                                                << 
466   // Free electron gas model                      571   // Free electron gas model
467   if ( T < 0.001 ) {                              572   if ( T < 0.001 ) {
468     G4double shigh = G4Log( 1.0 + x3*1000.0 +  << 573     G4double slow  = a[i][0] ;
469     ionloss  = slow*shigh*std::sqrt(T*1000.0)  << 574     G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
                                                   >> 575                    * a[i][2]*1000.0 ;
                                                   >> 576     ionloss  = slow*shigh / (slow + shigh) ;
                                                   >> 577     ionloss *= sqrt(T*1000.0) ;
470                                                   578 
471   // Main parametrisation                         579   // Main parametrisation
472   } else {                                        580   } else {
473     slow  *= G4Exp(G4Log(T*1000.0)*x1);        << 581     G4double slow  = a[i][0] * pow((T*1000.0), a[i][1]) ;
474     G4double shigh = G4Log( 1.0 + x3/T + x4*T  << 582     G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
475     ionloss = slow*shigh / (slow + shigh) ;       583     ionloss = slow*shigh / (slow + shigh) ;
476     /*                                            584     /*
477     G4cout << "## " << i << ". T= " << T << "     585     G4cout << "## " << i << ". T= " << T << " slow= " << slow
478            << " a0= " << a[i][0] << " a1= " <<    586            << " a0= " << a[i][0] << " a1= " << a[i][1] 
479            << " shigh= " << shigh                 587            << " shigh= " << shigh 
480            << " dedx= " << ionloss << " q^2= " << 588            << " dedx= " << ionloss << " q^2= " <<  HeEffChargeSquare(z, T*MeV) 
481            << G4endl;                          << 589      << G4endl;
482     */                                            590     */
483   }                                               591   }
484   ionloss = std::max(ionloss, 0.0);            << 592   if ( ionloss < 0.0) ionloss = 0.0 ;
                                                   >> 593 
                                                   >> 594   // He effective charge
                                                   >> 595   ionloss /= HeEffChargeSquare(z, T);
                                                   >> 596 
485   return ionloss;                                 597   return ionloss;
486 }                                                 598 }
487                                                   599 
488 //....oooOO0OOooo........oooOO0OOooo........oo    600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489                                                   601 
490 G4double G4BraggIonModel::HeDEDX(const G4Mater << 602 G4double G4BraggIonModel::DEDX(const G4Material* material,
491                                const G4double  << 603                                      G4double kineticEnergy)
492 {                                                 604 {
493   // aEnergy is energy of alpha                << 
494   G4double eloss = 0.0;                           605   G4double eloss = 0.0;
495   // check DB                                  << 606   const G4int numberOfElements = material->GetNumberOfElements();
496   if(material != currentMaterial) {            << 
497     currentMaterial = material;                << 
498     baseMaterial = material->GetBaseMaterial() << 
499       ? material->GetBaseMaterial() : material << 
500     iPSTAR    = -1;                            << 
501     iASTAR    = -1;                            << 
502     iMolecula = -1;                            << 
503     iICRU90 = (nullptr != fICRU90) ? fICRU90-> << 
504                                                << 
505     if(iICRU90 < 0) {                          << 
506       if(isAlpha) {                            << 
507   iASTAR = fASTAR->GetIndex(baseMaterial);     << 
508   if(iASTAR < 0) { iMolecula = HasMaterialForH << 
509       } else {                                 << 
510   iPSTAR = fPSTAR->GetIndex(baseMaterial);     << 
511       }                                        << 
512     }                                          << 
513     /*                                         << 
514     G4cout << "%%% " <<material->GetName() <<  << 
515            << iMolecula << "  iASTAR= " << iAS << 
516            << "  iICRU90= " << iICRU90<< G4end << 
517     */                                         << 
518   }                                            << 
519   // ICRU90                                    << 
520   if(iICRU90 >= 0) {                           << 
521     eloss = (isAlpha)                          << 
522       ? fICRU90->GetElectronicDEDXforAlpha(iIC << 
523       : fICRU90->GetElectronicDEDXforProton(iI << 
524     if(eloss > 0.0) { return eloss*material->G << 
525   }                                            << 
526   // PSTAR parameterisation                    << 
527   if( iPSTAR >= 0 ) {                          << 
528     return fPSTAR->GetElectronicDEDX(iPSTAR, a << 
529       *material->GetDensity();                 << 
530   }                                            << 
531   // ASTAR                                     << 
532   if( iASTAR >= 0 ) {                          << 
533     eloss = fASTAR->GetElectronicDEDX(iASTAR,  << 
534     /*                                         << 
535     G4cout << "ASTAR:  E=" << aEnergy          << 
536      << " dedx=" << eloss*material->GetDensity << 
537      << "  " << particle->GetParticleName() << << 
538     */                                         << 
539     if(eloss > 0.0) { return eloss*material->G << 
540   }                                            << 
541                                                << 
542   const std::size_t numberOfElements = materia << 
543   const G4ElementVector* theElmVector = materi << 
544   const G4double* theAtomicNumDensityVector =     607   const G4double* theAtomicNumDensityVector =
545     material->GetAtomicNumDensityVector();     << 608                                  material->GetAtomicNumDensityVector();
                                                   >> 609 
                                                   >> 610   // compaund material with parametrisation
                                                   >> 611   G4int iNist = astar.GetIndex(material);
546                                                   612 
547   // molecular data use proton stopping power  << 613   if( iNist >= 0 ) {
548   // element data from ICRU49 include data for << 614     G4double T = kineticEnergy*rateMassHe2p;
549   if(iMolecula >= 0) {                         << 615     return astar.GetElectronicDEDX(iNist, T)*material->GetDensity()/
550     const G4double zeff = material->GetTotNbOf << 616       HeEffChargeSquare(astar.GetEffectiveZ(iNist), T/MeV);
551       material->GetTotNbOfAtomsPerVolume();    << 
552     heChargeSquare = HeEffChargeSquare(zeff, a << 
553     eloss = HeStoppingPower(aEnergy)*heChargeS << 
554                                                   617 
555     // pure material                           << 618   } else if( HasMaterial(material) ) {
                                                   >> 619 
                                                   >> 620     eloss = StoppingPower(material, kineticEnergy)*
                                                   >> 621       material->GetDensity()/amu;
                                                   >> 622 
                                                   >> 623   // pure material
556   } else if(1 == numberOfElements) {              624   } else if(1 == numberOfElements) {
557                                                   625 
558     const G4Element* element = (*theElmVector) << 626     G4double z = material->GetZ();
559     eloss = HeElectronicStoppingPower(element- << 627     eloss = ElectronicStoppingPower(z, kineticEnergy)
560       * (material->GetTotNbOfAtomsPerVolume()) << 628                                * (material->GetTotNbOfAtomsPerVolume());
561                                                   629 
562   // Brugg's rule calculation                     630   // Brugg's rule calculation
563   } else {                                        631   } else {
                                                   >> 632     const G4ElementVector* theElementVector =
                                                   >> 633                            material->GetElementVector() ;
                                                   >> 634 
564     //  loop for the elements in the material     635     //  loop for the elements in the material
565     for (std::size_t i=0; i<numberOfElements;  << 636     for (G4int i=0; i<numberOfElements; i++)
566       const G4Element* element = (*theElmVecto << 637     {
567       eloss += HeElectronicStoppingPower(eleme << 638       const G4Element* element = (*theElementVector)[i] ;
568   * theAtomicNumDensityVector[i];              << 639       eloss   += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
                                                   >> 640                                    * theAtomicNumDensityVector[i];
569     }                                             641     }
570   }                                               642   }
571   return eloss*theZieglerFactor;                  643   return eloss*theZieglerFactor;
572 }                                                 644 }
573                                                   645 
574 //....oooOO0OOooo........oooOO0OOooo........oo    646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
575                                                   647 
576 G4double                                       << 648 G4double G4BraggIonModel::HeEffChargeSquare(G4double z, 
577 G4BraggIonModel::HeEffChargeSquare(const G4dou << 649                                             G4double kinEnergyHeInMeV) const
578                                    const G4dou << 
579 {                                                 650 {
580   // The aproximation of He effective charge f    651   // The aproximation of He effective charge from:
581   // J.F.Ziegler, J.P. Biersack, U. Littmark      652   // J.F.Ziegler, J.P. Biersack, U. Littmark
582   // The Stopping and Range of Ions in Matter,    653   // The Stopping and Range of Ions in Matter,
583   // Vol.1, Pergamon Press, 1985                  654   // Vol.1, Pergamon Press, 1985
584                                                   655 
585   static const G4double c[6] = {0.2865,  0.126 << 656   static G4double c[6] = {0.2865,  0.1266, -0.001429,
586                                 0.02402,-0.011 << 657                           0.02402,-0.01135, 0.001475};
587                                                   658 
588   G4double e = std::max(0.0, G4Log(kinEnergyHe << 659   G4double e = std::max(0.0,std::log(kinEnergyHeInMeV*massFactor));
589   G4double x = c[0] ;                             660   G4double x = c[0] ;
590   G4double y = 1.0 ;                              661   G4double y = 1.0 ;
591   for (G4int i=1; i<6; ++i) {                  << 662   for (G4int i=1; i<6; i++) {
592     y *= e;                                    << 663     y *= e ;
593     x += y * c[i];                             << 664     x += y * c[i] ;
594   }                                               665   }
595                                                   666 
596   G4double w = 7.6 -  e ;                         667   G4double w = 7.6 -  e ;
597   w = 1.0 + (0.007 + 0.00005*z) * G4Exp( -w*w  << 668   w = 1.0 + (0.007 + 0.00005*z) * exp( -w*w ) ;
598   w = 4.0 * (1.0 - G4Exp(-x)) * w * w ;        << 669   w = 4.0 * (1.0 - exp(-x)) * w * w ;
599                                                   670 
600   return w;                                       671   return w;
601 }                                                 672 }
602                                                   673 
603 //....oooOO0OOooo........oooOO0OOooo........oo    674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
604                                                   675 
605                                                   676