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 11.1)


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