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