| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // ------------------------------------------- <<
27 // 23 //
28 // GEANT4 Class file << 24 // $Id: G4hIonisation.cc,v 1.23 2001/11/09 13:59:47 maire Exp $
>> 25 // GEANT4 tag $Name: geant4-04-00 $
29 // 26 //
>> 27 //---------------- G4hIonisation physics process -------------------------------
>> 28 // by Laszlo Urban, 30 May 1997
>> 29 //------------------------------------------------------------------------------
30 // 30 //
31 // File name: G4hIonisation << 31 // corrected by L.Urban on 24/09/97
>> 32 // several bugs corrected by L.Urban on 13/01/98
>> 33 // 07-04-98 remove 'tracking cut' of the ionizing particle, mma
>> 34 // 22-10-98 cleanup L.Urban
>> 35 // 02-02-99 bugs fixed , L.Urban
>> 36 // 29-07-99 correction in BuildLossTable for low energy, L.Urban
>> 37 // 10-02-00 modifications , new e.m. structure, L.Urban
>> 38 // 10-08-00 V.Ivanchenko change BuildLambdaTable, in order to
>> 39 // simulate energy losses of ions; correction to
>> 40 // cross section for particles with spin 1 is inserted as well
>> 41 // 28-05-01 V.Ivanchenko minor changes to provide ANSI -wall compilation
>> 42 // 10-08-01 new methods Store/Retrieve PhysicsTable (mma)
>> 43 // 14-08-01 new function ComputeRestrictedMeandEdx() + 'cleanup' (mma)
>> 44 // 29-08-01 PostStepDoIt: correction for spin 1/2 (instead of 1) (mma)
>> 45 // 17-09-01 migration of Materials to pure STL (mma)
>> 46 // 25-09-01 completion of RetrievePhysicsTable() (mma)
>> 47 // 29-10-01 all static functions no more inlined
>> 48 // 08-11-01 Charge renamed zparticle; added to the dedx
32 // 49 //
33 // Author: Laszlo Urban << 50 //------------------------------------------------------------------------------
34 // << 51
35 // Creation date: 30.05.1997 << 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36 // << 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
37 // Modified by Laszlo Urban, Michel Maire and <<
38 // <<
39 // ------------------------------------------- <<
40 // <<
41 //....oooOO0OOooo........oooOO0OOooo........oo <<
42 //....oooOO0OOooo........oooOO0OOooo........oo <<
43 54
44 #include "G4hIonisation.hh" 55 #include "G4hIonisation.hh"
45 #include "G4PhysicalConstants.hh" << 56 #include "G4UnitsTable.hh"
46 #include "G4SystemOfUnits.hh" <<
47 #include "G4Electron.hh" <<
48 #include "G4Proton.hh" <<
49 #include "G4AntiProton.hh" <<
50 #include "G4BraggModel.hh" <<
51 #include "G4BetheBlochModel.hh" <<
52 #include "G4EmStandUtil.hh" <<
53 #include "G4PionPlus.hh" <<
54 #include "G4PionMinus.hh" <<
55 #include "G4KaonPlus.hh" <<
56 #include "G4KaonMinus.hh" <<
57 #include "G4ICRU73QOModel.hh" <<
58 #include "G4EmParameters.hh" <<
59 <<
60 //....oooOO0OOooo........oooOO0OOooo........oo <<
61 <<
62 G4hIonisation::G4hIonisation(const G4String& n <<
63 : G4VEnergyLossProcess(name) <<
64 { <<
65 SetProcessSubType(fIonisation); <<
66 SetSecondaryParticle(G4Electron::Electron()) <<
67 eth = 2*CLHEP::MeV; <<
68 } <<
69 57
70 //....oooOO0OOooo........oooOO0OOooo........oo << 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 59
72 G4bool G4hIonisation::IsApplicable(const G4Par << 60 G4double G4hIonisation::LowerBoundLambda = 1.*keV;
>> 61 G4double G4hIonisation::UpperBoundLambda = 100.*TeV;
>> 62 G4int G4hIonisation::NbinLambda = 100;
>> 63
>> 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 65
>> 66 G4hIonisation::G4hIonisation(const G4String& processName)
>> 67 : G4VhEnergyLoss(processName),
>> 68 theMeanFreePathTable(0),
>> 69 Tmincut(1*keV)
>> 70 { }
>> 71
>> 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 73
>> 74 G4hIonisation::~G4hIonisation()
73 { 75 {
74 return true; << 76 if (theMeanFreePathTable) {
>> 77 theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
75 } 78 }
76 79
77 //....oooOO0OOooo........oooOO0OOooo........oo << 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 81
79 G4double G4hIonisation::MinPrimaryEnergy(const << 82 void G4hIonisation::SetLowerBoundLambda(G4double val)
80 const G4Material*, << 83 {LowerBoundLambda = val;}
81 G4double cut) << 84
82 { << 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 G4double x = 0.5*cut/electron_mass_c2; << 86
84 G4double gam = x*ratio + std::sqrt((1. + x)* << 87 void G4hIonisation::SetUpperBoundLambda(G4double val)
85 return mass*(gam - 1.0); << 88 {UpperBoundLambda = val;}
86 } << 89
87 << 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 //....oooOO0OOooo........oooOO0OOooo........oo << 91
89 << 92 void G4hIonisation::SetNbinLambda(G4int n)
90 void G4hIonisation::InitialiseEnergyLossProces << 93 {NbinLambda = n;}
91 const G4ParticleDefinition* part, << 94
92 const G4ParticleDefinition* bpart) << 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 { << 96
94 if(!isInitialised) { << 97 G4double G4hIonisation::GetLowerBoundLambda()
95 << 98 {return LowerBoundLambda;}
96 const G4ParticleDefinition* theBaseParticl << 99
97 G4String pname = part->GetParticleName(); << 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 G4double q = part->GetPDGCharge(); << 101
99 << 102 G4double G4hIonisation::GetUpperBoundLambda()
100 //G4cout << " G4hIonisation::InitialiseEne << 103 {return UpperBoundLambda;}
101 // << " " << bpart << G4endl; << 104
102 << 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 // define base particle << 106
104 if(part == bpart) { << 107 G4int G4hIonisation::GetNbinLambda()
105 theBaseParticle = nullptr; << 108 {return NbinLambda;}
106 } else if(nullptr != bpart) { <<
107 theBaseParticle = bpart; <<
108 <<
109 } else if(pname == "proton" || pname == "a <<
110 pname == "pi+" || pname == "pi-" || <<
111 pname == "kaon+" || pname == "kaon-" | <<
112 pname == "GenericIon" || pname == "alp <<
113 // no base particles <<
114 theBaseParticle = nullptr; <<
115 <<
116 } else { <<
117 // select base particle <<
118 if(part->GetPDGSpin() == 0.0) { <<
119 if(q > 0.0) { theBaseParticle = G4KaonPlus:: <<
120 else { theBaseParticle = G4KaonMinus::KaonMi <<
121 } else { <<
122 if(q > 0.0) { theBaseParticle = G4Proton::Pr <<
123 else { theBaseParticle = G4AntiProton::AntiP <<
124 } <<
125 } <<
126 SetBaseParticle(theBaseParticle); <<
127 109
128 // model limit defined for protons << 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 mass = part->GetPDGMass(); << 111
130 ratio = electron_mass_c2/mass; << 112 void G4hIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
131 eth = 2.0*MeV*mass/proton_mass_c2; << 113 // just call BuildLossTable+BuildLambdaTable
132 << 114 {
133 G4EmParameters* param = G4EmParameters::In << 115 // get bining from EnergyLoss
134 G4double emin = param->MinKinEnergy(); << 116 LowestKineticEnergy = GetLowerBoundEloss();
135 G4double emax = param->MaxKinEnergy(); << 117 HighestKineticEnergy = GetUpperBoundEloss();
136 << 118 TotBin = GetNbinEloss();
137 // define model of energy loss fluctuation << 119
138 if (nullptr == FluctModel()) { << 120 G4double* ElectronCutInRange = G4Electron::Electron()->GetLengthCuts();
139 G4bool ion = (pname == "GenericIon" || p << 121
140 SetFluctModel(G4EmStandUtil::ModelOfFluc << 122 if (aParticleType.GetPDGCharge() > 0.)
>> 123 {
>> 124 if( !EqualCutVectors(ptableElectronCutInRange,ElectronCutInRange)
>> 125 || (theDEDXpTable == NULL))
>> 126 {
>> 127 BuildLossTable(aParticleType);
>> 128 RecorderOfpProcess[CounterOfpProcess] = theLossTable;
>> 129 CounterOfpProcess++;
>> 130 }
>> 131 }
>> 132 else
>> 133 {
>> 134 if( !EqualCutVectors(pbartableElectronCutInRange,ElectronCutInRange)
>> 135 || (theDEDXpbarTable == NULL))
>> 136 {
>> 137 BuildLossTable(aParticleType) ;
>> 138 RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable;
>> 139 CounterOfpbarProcess++;
141 } 140 }
>> 141 }
>> 142
>> 143 BuildLambdaTable(aParticleType);
>> 144
>> 145 BuildDEDXTable(aParticleType);
>> 146
>> 147 if (&aParticleType == G4Proton::Proton()) PrintInfoDefinition();
>> 148 }
>> 149
>> 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 151
>> 152 void G4hIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
>> 153 {
>> 154 // Build tables of dE/dx due to the ionization process
>> 155 // the tables are built for *MATERIALS*
142 156
143 if (nullptr == EmModel(0)) { << 157 // create table
144 if(q > 0.0) { SetEmModel(new G4BraggMode << 158 //
145 else { SetEmModel(new G4ICRU73QOM << 159 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 160 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
>> 161
>> 162 if (theLossTable) {theLossTable->clearAndDestroy(); delete theLossTable;}
>> 163 theLossTable = new G4PhysicsTable(numOfMaterials);
>> 164
>> 165 // get delta cut in energy
>> 166 G4double* DeltaCutInKinEnergy = (G4Electron::Electron())->GetEnergyCuts();
>> 167
>> 168 // loop for materials
>> 169 //
>> 170 for (G4int J=0; J<numOfMaterials; J++)
>> 171 {
>> 172 // create physics vector and fill it
>> 173 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 174 LowestKineticEnergy, HighestKineticEnergy, TotBin);
>> 175
>> 176 const G4Material* material= (*theMaterialTable)[J];
>> 177
>> 178 // get electron cut in kinetic energy for the material
>> 179 G4double DeltaThreshold = G4std::max(DeltaCutInKinEnergy[J],Tmincut);
>> 180
>> 181 // now comes the loop for the kinetic energy values
>> 182 //
>> 183 for (G4int i = 0 ; i < TotBin ; i++)
>> 184 {
>> 185 G4double dEdx = ComputeRestrictedMeandEdx(aParticleType,
>> 186 aVector->GetLowEdgeEnergy(i),
>> 187 material,
>> 188 DeltaThreshold);
>> 189 aVector->PutValue(i,dEdx);
146 } 190 }
147 // to compute ranges correctly we have to << 191 theLossTable->insert(aVector);
148 // model even if activation limit is high << 192 }
149 EmModel(0)->SetLowEnergyLimit(emin); << 193 }
150 << 194
151 // high energy limit may be eth or DBL_MAX << 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 G4double emax1 = (EmModel(0)->HighEnergyLi << 196
153 EmModel(0)->SetHighEnergyLimit(emax1); << 197 void G4hIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
154 AddEmModel(1, EmModel(0), FluctModel()); << 198 {
155 << 199 // Build mean free path tables for the delta ray production process
156 // second model is used if the first does << 200 // tables are built for MATERIALS
157 if(emax1 < emax) { << 201
158 if (nullptr == EmModel(1)) { SetEmModel( << 202 //create table
159 EmModel(1)->SetLowEnergyLimit(emax1); << 203 //
160 << 204
161 // for extremely heavy particles upper l << 205 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
162 // should be increased << 206 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
163 emax = std::max(emax, eth*10); << 207
164 EmModel(1)->SetHighEnergyLimit(emax); << 208 if (theMeanFreePathTable)
165 AddEmModel(2, EmModel(1), FluctModel()); << 209 {theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
>> 210
>> 211 theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
>> 212
>> 213 // get electron cut in kinetic energy
>> 214 G4double* DeltaCutInKinEnergy = (G4Electron::Electron())->GetEnergyCuts() ;
>> 215
>> 216 // loop for materials
>> 217
>> 218 for (G4int J=0 ; J < numOfMaterials; J++)
>> 219 {
>> 220 //create physics vector then fill it ....
>> 221 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 222 LowerBoundLambda,UpperBoundLambda,NbinLambda);
>> 223
>> 224 // compute the (macroscopic) cross section first
>> 225
>> 226 const G4Material* material= (*theMaterialTable)[J];
>> 227 const G4ElementVector* theElementVector = material->GetElementVector();
>> 228 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
>> 229 G4int NumberOfElements = material->GetNumberOfElements();
>> 230
>> 231 // get the electron kinetic energy cut for the actual material,
>> 232 // it will be used in ComputeCrossSectionPerAtom
>> 233 // ( --> it will be the same for all the elements in this material)
>> 234 G4double DeltaThreshold =G4std::max(DeltaCutInKinEnergy[J],Tmincut);
>> 235
>> 236 for ( G4int i = 0 ; i < NbinLambda ; i++ )
>> 237 {
>> 238 G4double LowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
>> 239 G4double sigma = 0.;
>> 240 for (G4int iel=0; iel<NumberOfElements; iel++ )
>> 241 {
>> 242 sigma += NbOfAtomsPerVolume[iel]*
>> 243 ComputeCrossSectionPerAtom(aParticleType,
>> 244 LowEdgeEnergy,
>> 245 (*theElementVector)[iel]->GetZ(),
>> 246 DeltaThreshold);
>> 247 }
>> 248
>> 249 // mean free path = 1./macroscopic cross section
>> 250 G4double Value = sigma > DBL_MIN ? 1./sigma : DBL_MAX;
>> 251 aVector->PutValue(i, Value) ;
>> 252 }
>> 253 theMeanFreePathTable->insert(aVector);
166 } 254 }
167 isInitialised = true; << 255 }
>> 256
>> 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 258
>> 259 G4double G4hIonisation::ComputeRestrictedMeandEdx (
>> 260 const G4ParticleDefinition& aParticleType,
>> 261 G4double KineticEnergy,
>> 262 const G4Material* material,
>> 263 G4double DeltaThreshold)
>> 264 {
>> 265 // calculate the dE/dx due to the ionization process (Geant4 internal units)
>> 266 // Bethe-Bloch formula
>> 267 //
>> 268 G4double particleMass = aParticleType.GetPDGMass();
>> 269 G4double particleZ = aParticleType.GetPDGCharge()/eplus;
>> 270 G4double zsquare = particleZ*particleZ;
>> 271
>> 272 G4double ElectronDensity = material->GetElectronDensity();
>> 273 G4double Eexc = material->GetIonisation()->GetMeanExcitationEnergy();
>> 274 G4double Eexc2 = Eexc*Eexc;
>> 275
>> 276 G4double tau = KineticEnergy/particleMass;
>> 277 G4double gamma = tau + 1., bg2 = tau*(tau+2.), beta2 = bg2/(gamma*gamma);
>> 278 G4double RateMass = electron_mass_c2/particleMass;
>> 279 G4double Tmax=2.*electron_mass_c2*bg2/(1.+2.*gamma*RateMass+RateMass*RateMass);
>> 280
>> 281 G4double taul = material->GetIonisation()->GetTaul();
>> 282
>> 283 G4double dEdx = 0.;
>> 284 //
>> 285 // high energy part , Bethe-Bloch formula
>> 286 //
>> 287 if (tau > taul)
>> 288 {
>> 289 G4double rcut = G4std::min(DeltaThreshold/Tmax, 1.);
>> 290 dEdx = log(2.*electron_mass_c2*bg2*Tmax/Eexc2)
>> 291 +log(rcut)-(1.+rcut)*beta2;
>> 292
>> 293 //density correction
>> 294 G4double Cden = material->GetIonisation()->GetCdensity();
>> 295 G4double Mden = material->GetIonisation()->GetMdensity();
>> 296 G4double Aden = material->GetIonisation()->GetAdensity();
>> 297 G4double X0den = material->GetIonisation()->GetX0density();
>> 298 G4double X1den = material->GetIonisation()->GetX1density();
>> 299
>> 300 const G4double twoln10 = 2.*log(10.);
>> 301 G4double x = log(bg2)/twoln10;
>> 302 G4double delta;
>> 303 if (x < X0den) delta = 0.;
>> 304 else {delta = twoln10*x - Cden;
>> 305 if (x < X1den) delta += Aden*pow((X1den-x),Mden);
>> 306 }
>> 307
>> 308 // shell correction
>> 309 G4double* ShellCorrectionVector = material->GetIonisation()->
>> 310 GetShellCorrectionVector();
>> 311 const G4double bg2lim = 0.0169, taulim = 8.4146e-3;
>> 312 G4double sh = 0., xs = 1.;
>> 313 if (bg2 > bg2lim) for (G4int k=0; k<3; k++)
>> 314 {xs *= bg2; sh += ShellCorrectionVector[k]/xs;}
>> 315 else { for (G4int k=0; k<3; k++)
>> 316 {xs *= bg2lim; sh += ShellCorrectionVector[k]/xs;}
>> 317 sh *= log(tau/taul)/log(taulim/taul);
>> 318 }
>> 319
>> 320 // now you can compute the total ionization loss
>> 321 dEdx -= (delta + sh);
>> 322 dEdx *= twopi_mc2_rcl2*ElectronDensity*zsquare/beta2;
>> 323 if (dEdx < 0.) dEdx = 0.;
>> 324 }
>> 325 //
>> 326 // low energy part , parametrized energy loss formulae
>> 327 //
>> 328 if (tau <= taul)
>> 329 {
>> 330 // get elements in the actual material,
>> 331 const G4ElementVector* theElementVector = material->GetElementVector();
>> 332 const G4double* NbOfAtomsPerVolume=material->GetVecNbOfAtomsPerVolume();
>> 333 G4int NumberOfElements = material->GetNumberOfElements();
>> 334
>> 335 // loop for the elements in the material
>> 336 dEdx = 0.;
>> 337 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 338 {
>> 339 const G4Element* element = (*theElementVector)[iel];
>> 340 if (tau < element->GetIonisation()->GetTau0())
>> 341 dEdx += NbOfAtomsPerVolume[iel]
>> 342 *(element->GetIonisation()->GetAlow()*sqrt(tau)
>> 343 + element->GetIonisation()->GetBlow()*tau);
>> 344 else
>> 345 dEdx += NbOfAtomsPerVolume[iel]
>> 346 * element->GetIonisation()->GetClow()/sqrt(tau);
>> 347 }
>> 348 G4double deltaloss = 0.;
>> 349 if (DeltaThreshold < Tmax)
>> 350 {
>> 351 deltaloss = log(Tmax/DeltaThreshold)-
>> 352 beta2*(1.-DeltaThreshold/Tmax) ;
>> 353 if (aParticleType.GetPDGSpin() == 0.5)
>> 354 deltaloss += 0.25*(Tmax-DeltaThreshold)*(Tmax-DeltaThreshold)/
>> 355 (KineticEnergy*KineticEnergy+proton_mass_c2*proton_mass_c2);
>> 356 deltaloss *= twopi_mc2_rcl2*ElectronDensity*zsquare/beta2;
>> 357 }
>> 358 dEdx -= deltaloss;
>> 359 if (dEdx < 0.) dEdx = 0.;
>> 360 }
>> 361 return dEdx;
>> 362 }
>> 363
>> 364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 365
>> 366 G4double G4hIonisation::ComputeCrossSectionPerAtom(
>> 367 const G4ParticleDefinition& aParticleType,
>> 368 G4double KineticEnergy,
>> 369 G4double AtomicNumber,
>> 370 G4double DeltaThreshold)
>> 371 {
>> 372 // calculates the totalcross section per atom in GEANT4 internal units
>> 373 // ( it is called for elements , AtomicNumber = Z )
>> 374 //
>> 375 // nb: cross section formula is OK for spin=0 and 1/2 only !
>> 376
>> 377 G4double particleMass = aParticleType.GetPDGMass();
>> 378 G4double particleZ = aParticleType.GetPDGCharge()/eplus;
>> 379 G4double zparticle2 = particleZ*particleZ;
>> 380
>> 381 G4double TotalEnergy = KineticEnergy + particleMass;
>> 382
>> 383 G4double betasquare = KineticEnergy*(TotalEnergy+particleMass)
>> 384 /(TotalEnergy*TotalEnergy);
>> 385 G4double tempvar = particleMass+electron_mass_c2;
>> 386 G4double MaxKineticEnergyTransfer = 2.*electron_mass_c2*KineticEnergy
>> 387 *(TotalEnergy+particleMass)
>> 388 /(tempvar*tempvar+2.*electron_mass_c2*KineticEnergy);
>> 389
>> 390 G4double TotalCrossSection = 0.;
>> 391 if (MaxKineticEnergyTransfer > DeltaThreshold)
>> 392 {
>> 393 tempvar = DeltaThreshold/MaxKineticEnergyTransfer;
>> 394 TotalCrossSection = (1.-tempvar*(1.-betasquare*log(tempvar)))
>> 395 /DeltaThreshold;
>> 396
>> 397 G4double spin = aParticleType.GetPDGSpin();
>> 398 if (spin == 0.5) TotalCrossSection += 0.5
>> 399 *(MaxKineticEnergyTransfer-DeltaThreshold)
>> 400 /(TotalEnergy*TotalEnergy);
>> 401 if (spin == 1.) TotalCrossSection +=
>> 402 -log(tempvar)/(3.0*DeltaThreshold) +
>> 403 (MaxKineticEnergyTransfer - DeltaThreshold) *
>> 404 ((5.0+ 1.0/tempvar)*0.25 / (TotalEnergy*TotalEnergy) -
>> 405 betasquare /
>> 406 (MaxKineticEnergyTransfer * DeltaThreshold)) / 3.0;
>> 407
>> 408 TotalCrossSection *= twopi_mc2_rcl2*AtomicNumber*zparticle2/betasquare;
>> 409 }
>> 410 return TotalCrossSection;
>> 411 }
>> 412
>> 413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 414
>> 415 G4VParticleChange* G4hIonisation::PostStepDoIt(const G4Track& trackData,
>> 416 const G4Step& stepData)
>> 417 {
>> 418 aParticleChange.Initialize(trackData);
>> 419
>> 420 G4Material* aMaterial = trackData.GetMaterial();
>> 421 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
>> 422
>> 423 G4double particleMass = aParticle->GetDefinition()->GetPDGMass();
>> 424 G4double KineticEnergy = aParticle->GetKineticEnergy();
>> 425 G4double TotalEnergy = KineticEnergy + particleMass;
>> 426 G4double Psquare = KineticEnergy*(TotalEnergy+particleMass);
>> 427 G4double Esquare = TotalEnergy*TotalEnergy;
>> 428 G4double betasquare=Psquare/Esquare;
>> 429 G4double summass = particleMass + electron_mass_c2;
>> 430 G4double MaxKineticEnergyTransfer = 2.*electron_mass_c2*Psquare
>> 431 /(summass*summass+2.*electron_mass_c2*KineticEnergy);
>> 432 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
>> 433
>> 434 // get electron cut in kinetic energy
>> 435 G4double* DeltaCutInKinEnergy = (G4Electron::Electron())->GetEnergyCuts();
>> 436 G4double DeltaThreshold =
>> 437 G4std::max(DeltaCutInKinEnergy[aMaterial->GetIndex()],Tmincut);
>> 438
>> 439 // sampling kinetic energy of the delta ray
>> 440 //
>> 441 if (MaxKineticEnergyTransfer <= DeltaThreshold)
>> 442 // pathological case (it should not happen, there is no change at all)
>> 443 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 444
>> 445 // normal case
>> 446 G4double xc = DeltaThreshold/MaxKineticEnergyTransfer;
>> 447 G4double rate = MaxKineticEnergyTransfer/TotalEnergy;
>> 448 G4double te2 = 0.;
>> 449 if (aParticle->GetDefinition()->GetPDGSpin() == 0.5) te2=0.5*rate*rate;
>> 450
>> 451 // sampling follows ...
>> 452 G4double x,grej;
>> 453 G4double grejc=1.-betasquare*xc+te2*xc*xc;
>> 454 do { x=xc/(1.-(1.-xc)*G4UniformRand());
>> 455 grej=(1.-x*(betasquare-x*te2))/grejc;
>> 456 } while(G4UniformRand() > grej);
>> 457
>> 458 G4double DeltaKineticEnergy = x * MaxKineticEnergyTransfer;
>> 459
>> 460 if (DeltaKineticEnergy <= 0.)
>> 461 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 462
>> 463 G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
>> 464 2. * electron_mass_c2 ));
>> 465 G4double TotalMomentum = sqrt(Psquare);
>> 466 G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
>> 467 /(DeltaTotalMomentum * TotalMomentum);
>> 468
>> 469 if (costheta < -1.) costheta = -1.;
>> 470 if (costheta > +1.) costheta = +1.;
>> 471
>> 472 // direction of the delta electron
>> 473 //
>> 474 G4double phi = twopi*G4UniformRand();
>> 475 G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
>> 476 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
>> 477
>> 478 G4ThreeVector DeltaDirection(dirx,diry,dirz);
>> 479 DeltaDirection.rotateUz(ParticleDirection);
>> 480
>> 481 // create G4DynamicParticle object for delta ray
>> 482 //
>> 483 G4DynamicParticle *theDeltaRay = new G4DynamicParticle;
>> 484 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
>> 485 theDeltaRay->SetMomentumDirection(
>> 486 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z());
>> 487 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 488
>> 489 // fill aParticleChange
>> 490 //
>> 491 G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy;
>> 492 G4double Edep = 0;
>> 493
>> 494 if (finalKineticEnergy > MinKineticEnergy)
>> 495 {
>> 496 G4double finalPx = TotalMomentum*ParticleDirection.x()
>> 497 - DeltaTotalMomentum*DeltaDirection.x();
>> 498 G4double finalPy = TotalMomentum*ParticleDirection.y()
>> 499 - DeltaTotalMomentum*DeltaDirection.y();
>> 500 G4double finalPz = TotalMomentum*ParticleDirection.z()
>> 501 - DeltaTotalMomentum*DeltaDirection.z();
>> 502 G4double finalMomentum =
>> 503 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
>> 504 finalPx /= finalMomentum;
>> 505 finalPy /= finalMomentum;
>> 506 finalPz /= finalMomentum;
>> 507
>> 508 aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
>> 509 }
>> 510 else
>> 511 {
>> 512 Edep = finalKineticEnergy;
>> 513 finalKineticEnergy = 0.;
>> 514 if (aParticle->GetDefinition()->GetParticleName() == "proton")
>> 515 aParticleChange.SetStatusChange(fStopAndKill);
>> 516 else aParticleChange.SetStatusChange(fStopButAlive);
>> 517 }
>> 518
>> 519 aParticleChange.SetEnergyChange( finalKineticEnergy );
>> 520 aParticleChange.SetNumberOfSecondaries(1);
>> 521 aParticleChange.AddSecondary(theDeltaRay);
>> 522 aParticleChange.SetLocalEnergyDeposit (Edep);
>> 523
>> 524 //ResetNumberOfInteractionLengthLeft();
>> 525 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 526 }
>> 527
>> 528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 529
>> 530 G4bool G4hIonisation::StorePhysicsTable(G4ParticleDefinition* particle,
>> 531 const G4String& directory,
>> 532 G4bool ascii)
>> 533 {
>> 534 G4String particleName = particle->GetParticleName();
>> 535 G4String filename;
>> 536
>> 537 // store stopping power table
>> 538 if ((particleName == "proton")||(particleName == "anti_proton")){
>> 539 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 540 if ( !theLossTable->StorePhysicsTable(filename, ascii) ){
>> 541 G4cout << " FAIL theLossTable->StorePhysicsTable in " << filename
>> 542 << G4endl;
>> 543 return false;
>> 544 }}
>> 545
>> 546 // store mean free path table
>> 547 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 548 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 549 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 550 << G4endl;
>> 551 return false;
168 } 552 }
>> 553
>> 554 G4cout << GetProcessName() << " for " << particleName
>> 555 << ": Success to store the PhysicsTables in "
>> 556 << directory << G4endl;
>> 557 return true;
169 } 558 }
170 559
171 //....oooOO0OOooo........oooOO0OOooo........oo << 560 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 561
>> 562 G4bool G4hIonisation::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 563 const G4String& directory,
>> 564 G4bool ascii)
>> 565 {
>> 566 // delete theLossTable and theMeanFreePathTable
>> 567 if (theLossTable != 0) {
>> 568 theLossTable->clearAndDestroy();
>> 569 delete theLossTable;
>> 570 }
>> 571 if (theMeanFreePathTable != 0) {
>> 572 theMeanFreePathTable->clearAndDestroy();
>> 573 delete theMeanFreePathTable;
>> 574 }
>> 575
>> 576 // get bining from EnergyLoss
>> 577 LowestKineticEnergy = GetLowerBoundEloss();
>> 578 HighestKineticEnergy = GetUpperBoundEloss();
>> 579 TotBin = GetNbinEloss();
>> 580
>> 581 G4String particleName = particle->GetParticleName();
>> 582 G4String filename;
>> 583
>> 584 // retreive stopping power table
>> 585 if ((particleName == "proton")||(particleName == "anti_proton")){
>> 586 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 587 theLossTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 588 if ( !theLossTable->RetrievePhysicsTable(filename, ascii) ){
>> 589 G4cout << " FAIL theLossTable0->RetrievePhysicsTable in " << filename
>> 590 << G4endl;
>> 591 return false;
>> 592 }
>> 593 if (particleName == "proton")
>> 594 RecorderOfpProcess[CounterOfpProcess++] = theLossTable;
>> 595 if (particleName == "anti_proton")
>> 596 RecorderOfpbarProcess[CounterOfpbarProcess++] = theLossTable;
>> 597 }
>> 598
>> 599 // retreive mean free path table
>> 600 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 601 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 602 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 603 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 604 << G4endl;
>> 605 return false;
>> 606 }
>> 607
>> 608
>> 609 G4cout << GetProcessName() << " for " << particleName
>> 610 << ": Success to retrieve the PhysicsTables from "
>> 611 << directory << G4endl;
>> 612
>> 613 BuildDEDXTable(*particle);
>> 614
>> 615 return true;
>> 616 }
>> 617
>> 618 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 619
173 void G4hIonisation::ProcessDescription(std::os << 620 void G4hIonisation::PrintInfoDefinition()
174 { 621 {
175 out << " Hadron ionisation"; << 622 G4String comments = " Knock-on electron cross sections . "
176 G4VEnergyLossProcess::ProcessDescription(out << 623 "\n Good description above the mean excitation energy.\n"
>> 624 " delta ray energy sampled from differential Xsection.";
>> 625
>> 626 G4cout << G4endl << GetProcessName() << ": " << comments
>> 627 << "\n PhysicsTables from " << G4BestUnit(LowestKineticEnergy,
>> 628 "Energy")
>> 629 << " to " << G4BestUnit(HighestKineticEnergy,"Energy")
>> 630 << " in " << TotBin << " bins. \n";
177 } 631 }
178 632
179 //....oooOO0OOooo........oooOO0OOooo........oo << 633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 634