| 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 << 7 //
6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4ionIonisation.cc,v 1.3.8.1.2.1 1999/12/08 17:34:27 gunter Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-01-01 $
8 // * LICENSE and available at http://cern.ch/ << 10 //
9 // * include a list of copyright holders. << 11 // -------------------------------------------------------------
10 // * << 12 // GEANT 4 class implementation file
11 // * Neither the authors of this software syst << 13 //
12 // * institutes,nor the agencies providing fin << 14 // For information related to this code contact:
13 // * work make any representation or warran << 15 // CERN, IT Division, ASD group
14 // * regarding this software system or assum << 16 // History: based on object model of
15 // * use. Please see the license in the file << 17 // 2nd December 1995, G.Cosmo
16 // * for the full disclaimer and the limitatio << 18 // ---------- G4ionIonisation physics process -----------
17 // * << 19 // by Laszlo Urban, 08 Dec 1998
18 // * This code implementation is the result << 20 // **************************************************************
19 // * technical work of the GEANT4 collaboratio << 21 // It is the first implementation of the ionisation for IONS
20 // * By using, copying, modifying or distri << 22 // --------------------------------------------------------------
21 // * any work based on the software) you ag << 23
22 // * use in resulting scientific publicati <<
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* <<
25 // <<
26 // <<
27 // ------------------------------------------- <<
28 // <<
29 // GEANT4 Class file <<
30 // <<
31 // <<
32 // File name: G4ionIonisation <<
33 // <<
34 // Author: Vladimir Ivanchenko <<
35 // <<
36 // Creation date: 07.05.2002 <<
37 // <<
38 // Modifications: <<
39 // <<
40 // 23-12-02 Change interface in order to move <<
41 // 26-12-02 Secondary production moved to deri <<
42 // 13-02-03 SubCutoff regime is assigned to a <<
43 // 18-04-03 Use IonFluctuations (V.Ivanchenko) <<
44 // 03-08-03 Add effective charge (V.Ivanchenko <<
45 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossPro <<
46 // 27-05-04 Set integral to be a default regim <<
47 // 08-11-04 Migration to new interface of Stor <<
48 // 08-04-05 Major optimisation of internal int <<
49 // 10-01-06 SetStepLimits -> SetStepFunction ( <<
50 // 10-05-06 Add a possibility to download user <<
51 // 13-05-06 Add data for light ion stopping in <<
52 // 14-01-07 use SetEmModel() and SetFluctModel <<
53 // 16-05-07 Add data for light ion stopping on <<
54 // 07-11-07 Fill non-ionizing energy loss (V.I <<
55 // 12-09-08 Removed InitialiseMassCharge and C <<
56 // <<
57 // <<
58 // ------------------------------------------- <<
59 // <<
60 //....oooOO0OOooo........oooOO0OOooo........oo <<
61 //....oooOO0OOooo........oooOO0OOooo........oo <<
62 24
63 #include "G4ionIonisation.hh" 25 #include "G4ionIonisation.hh"
64 #include "G4PhysicalConstants.hh" << 26 #include "G4UnitsTable.hh"
65 #include "G4SystemOfUnits.hh" <<
66 #include "G4Electron.hh" <<
67 #include "G4GenericIon.hh" <<
68 #include "G4BraggModel.hh" <<
69 #include "G4BraggIonModel.hh" <<
70 #include "G4BetheBlochModel.hh" <<
71 #include "G4LossTableManager.hh" <<
72 #include "G4EmParameters.hh" <<
73 #include "G4EmStandUtil.hh" <<
74 27
75 //....oooOO0OOooo........oooOO0OOooo........oo << 28 // constructor and destructor
>> 29
>> 30 G4ionIonisation::G4ionIonisation(const G4String& processName)
>> 31 : G4VContinuousDiscreteProcess(processName),
>> 32 ParticleMass(proton_mass_c2),Charge(eplus),
>> 33 dEdx(1.*MeV/mm),MinKineticEnergy(1.*keV)
>> 34 { PrintInfoDefinition() ; }
>> 35
>> 36
>> 37 G4ionIonisation::~G4ionIonisation()
>> 38 { }
76 39
77 G4ionIonisation::G4ionIonisation(const G4Strin << 40 G4double G4ionIonisation::GetConstraints(const G4DynamicParticle *aParticle,
78 : G4VEnergyLossProcess(name) << 41 G4Material *aMaterial)
79 { 42 {
80 SetLinearLossLimit(0.02); << 43 // returns the Step limit
81 SetProcessSubType(fIonisation); << 44 // dRoverRange is the max. allowed relative range loss in one step
82 SetSecondaryParticle(G4Electron::Electron()) << 45 // it calculates dEdx and the range as well....
83 eth = 2*CLHEP::MeV; << 46 const G4double minstep=0.01*mm ;
84 } <<
85 47
86 //....oooOO0OOooo........oooOO0OOooo........oo << 48 G4double KineticEnergy,StepLimit;
87 49
88 G4bool G4ionIonisation::IsApplicable(const G4P << 50 Charge = aParticle->GetDefinition()->GetPDGCharge()/eplus ;
89 { <<
90 return true; <<
91 } <<
92 51
93 //....oooOO0OOooo........oooOO0OOooo........oo << 52 KineticEnergy = aParticle->GetKineticEnergy();
94 53
95 G4double G4ionIonisation::MinPrimaryEnergy(con << 54 G4double massratio=proton_mass_c2/
96 const G4Material*, << 55 aParticle->GetDefinition()->GetPDGMass() ;
97 G4double cut) <<
98 { <<
99 return p->GetPDGMass()*(std::sqrt(1. + 0.5*c <<
100 } <<
101 56
102 //....oooOO0OOooo........oooOO0OOooo........oo << 57 G4double Tscaled= KineticEnergy*massratio ;
>> 58 G4double ChargeSquare = Charge*Charge ;
103 59
104 void G4ionIonisation::InitialiseEnergyLossProc << 60 dEdx=ComputedEdx(aParticle,aMaterial) ;
105 const G4ParticleDefinition* part, << 61 StepLimit = 0.2*KineticEnergy/dEdx ;
106 const G4ParticleDefinition* bpart) << 62 if(StepLimit < minstep)
>> 63 StepLimit = minstep ;
>> 64
>> 65 return StepLimit ;
>> 66 }
>> 67
>> 68 G4VParticleChange* G4ionIonisation::AlongStepDoIt(
>> 69 const G4Track& trackData,const G4Step& stepData)
>> 70 // compute the energy loss after a step
107 { 71 {
108 const G4ParticleDefinition* ion = G4GenericI << 72 const G4DynamicParticle* aParticle;
>> 73 G4Material* aMaterial;
>> 74 G4double E,finalT,Step,ChargeSquare,MeanLoss ;
>> 75
>> 76 aParticleChange.Initialize(trackData) ;
>> 77 aMaterial = trackData.GetMaterial() ;
>> 78
>> 79 // get the actual (true) Step length from stepData
>> 80 Step = stepData.GetStepLength() ;
>> 81
>> 82 aParticle = trackData.GetDynamicParticle() ;
>> 83 G4double massratio=proton_mass_c2/
>> 84 aParticle->GetDefinition()->GetPDGMass() ;
>> 85 ChargeSquare = Charge*Charge ;
>> 86
>> 87 G4int index = aMaterial->GetIndex() ;
>> 88 E = aParticle->GetKineticEnergy() ;
>> 89
>> 90 if(E < MinKineticEnergy) MeanLoss = E ;
>> 91 else
>> 92 {
>> 93 MeanLoss = Step*dEdx ;
>> 94 MeanLoss /= (massratio*ChargeSquare) ;
>> 95 }
>> 96 finalT = E - MeanLoss ;
109 97
110 if(!isInitialised) { << 98 if(finalT < MinKineticEnergy) finalT = 0. ;
111 theParticle = part; <<
112 99
113 // define base particle << 100 // kill the particle if the kinetic energy <= 0
114 const G4ParticleDefinition* theBaseParticl << 101 if (finalT <= 0. )
115 const G4int pdg = part->GetPDGEncoding(); << 102 {
>> 103 finalT = 0.;
>> 104 aParticleChange.SetStatusChange(fStopAndKill);
>> 105 }
116 106
117 if(part == bpart) { << 107 aParticleChange.SetEnergyChange( finalT ) ;
118 theBaseParticle = nullptr; << 108 aParticleChange.SetLocalEnergyDeposit(E-finalT) ;
119 } else if(nullptr != bpart) { << 109
120 theBaseParticle = bpart; << 110 return &aParticleChange ;
121 } else if(part == ion || pdg == 1000020040 << 111 }
122 theBaseParticle = nullptr; << 112
123 } else { << 113
124 theBaseParticle = ion; << 114 G4double G4ionIonisation::GetMeanFreePath(
>> 115 const G4Track& trackData,
>> 116 G4double previousStepSize,
>> 117 G4ForceCondition* condition)
>> 118 {
>> 119 const G4DynamicParticle* aParticle ;
>> 120 G4Material* aMaterial ;
>> 121 G4double MeanFreePath;
>> 122
>> 123 *condition = NotForced ;
>> 124
>> 125 aParticle = trackData.GetDynamicParticle() ;
>> 126 aMaterial = trackData.GetMaterial() ;
>> 127
>> 128 G4double KineticEnergy = aParticle->GetKineticEnergy() ;
>> 129 Charge=(aParticle->GetDefinition()->GetPDGCharge())/eplus;
>> 130 G4double ChargeSquare=Charge*Charge ;
>> 131
>> 132 // compute the (macroscopic) cross section first
>> 133
>> 134 const G4ElementVector* theElementVector=
>> 135 aMaterial->GetElementVector() ;
>> 136 const G4double* theAtomicNumDensityVector =
>> 137 aMaterial->GetAtomicNumDensityVector();
>> 138 const G4int NumberOfElements=
>> 139 aMaterial->GetNumberOfElements() ;
>> 140 G4int index = aMaterial->GetIndex() ;
>> 141
>> 142 DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy();
>> 143 DeltaCutInKineticEnergyNow = DeltaCutInKineticEnergy[index] ;
>> 144
>> 145 G4double sigma = 0. ;
>> 146 for (G4int iel=0; iel<NumberOfElements; iel++ )
>> 147 {
>> 148 sigma += theAtomicNumDensityVector[iel]*
>> 149 ComputeMicroscopicCrossSection(aParticle,
>> 150 KineticEnergy,
>> 151 (*theElementVector)(iel)->GetZ() ) ;
125 } 152 }
126 SetBaseParticle(theBaseParticle); <<
127 153
128 // model limit defined for protons << 154 sigma *= twopi_mc2_rcl2 * ChargeSquare ;
129 eth = 2*CLHEP::MeV*part->GetPDGMass()/CLHE << 155
>> 156 // mean free path = 1./macroscopic cross section
>> 157 MeanFreePath = sigma<=0 ? DBL_MAX : 1./sigma ;
130 158
131 G4EmParameters* param = G4EmParameters::In << 159 return MeanFreePath ;
132 G4double emin = param->MinKinEnergy(); << 160 }
133 G4double emax = param->MaxKinEnergy(); << 161
134 << 162 G4double G4ionIonisation::ComputeMicroscopicCrossSection(
135 // define model of energy loss fluctuation << 163 const G4DynamicParticle* aParticle,
136 if (nullptr == FluctModel()) { << 164 G4double KineticEnergy,
137 SetFluctModel(G4EmStandUtil::ModelOfFluc << 165 G4double AtomicNumber)
>> 166 {
>> 167 G4double TotalEnergy,
>> 168 betasquare,
>> 169 MaxKineticEnergyTransfer,TotalCrossSection,tempvar;
>> 170
>> 171 // get particle data ...................................
>> 172 ParticleMass = aParticle->GetDefinition()->GetPDGMass() ;
>> 173 TotalEnergy=KineticEnergy + ParticleMass;
>> 174
>> 175 // some kinematics......................
>> 176 betasquare = KineticEnergy*(TotalEnergy+ParticleMass)
>> 177 /(TotalEnergy*TotalEnergy);
>> 178 tempvar = ParticleMass+electron_mass_c2;
>> 179 MaxKineticEnergyTransfer = 2.*electron_mass_c2*KineticEnergy
>> 180 *(TotalEnergy+ParticleMass)
>> 181 /(tempvar*tempvar+2.*electron_mass_c2*KineticEnergy);
>> 182
>> 183 // now you can calculate the total cross section ------------------
>> 184 if( MaxKineticEnergyTransfer > DeltaCutInKineticEnergyNow )
>> 185 {
>> 186 tempvar=DeltaCutInKineticEnergyNow/MaxKineticEnergyTransfer;
>> 187 TotalCrossSection = (1.-tempvar*(1.-betasquare*log(tempvar)))
>> 188 /DeltaCutInKineticEnergyNow;
>> 189
>> 190 TotalCrossSection *= AtomicNumber/betasquare;
138 } 191 }
>> 192 else
>> 193 TotalCrossSection= 0. ;
>> 194
>> 195 return TotalCrossSection ;
>> 196 }
139 197
140 if (nullptr == EmModel(0)) { << 198 G4double G4ionIonisation::ComputedEdx(const G4DynamicParticle* aParticle,
141 if (pdg == 1000020040) { << 199 G4Material* material)
142 SetEmModel(new G4BraggIonModel()); << 200 {
143 } else { << 201 // cuts for electron ....................
144 SetEmModel(new G4BraggModel()); << 202 DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
>> 203
>> 204 G4double KineticEnergy , ionloss ;
>> 205 G4double RateMass ;
>> 206 G4bool isOutRange ;
>> 207 const G4double twoln10 = 2.*log(10.) ;
>> 208 const G4double Factor = twopi_mc2_rcl2 ;
>> 209 const G4double bg2lim = 0.0169 , taulim = 8.4146e-3 ;
>> 210
>> 211 RateMass = electron_mass_c2/proton_mass_c2 ;
>> 212
>> 213 // get material parameters needed for the energy loss calculation
>> 214
>> 215 G4double ElectronDensity,Eexc,Eexc2,Cden,Mden,Aden,X0den,X1den,taul ;
>> 216 G4double* ShellCorrectionVector;
>> 217
>> 218 ElectronDensity = material->GetElectronDensity();
>> 219 Eexc = material->GetIonisation()->GetMeanExcitationEnergy();
>> 220 Eexc2 = Eexc*Eexc ;
>> 221 Cden = material->GetIonisation()->GetCdensity();
>> 222 Mden = material->GetIonisation()->GetMdensity();
>> 223 Aden = material->GetIonisation()->GetAdensity();
>> 224 X0den = material->GetIonisation()->GetX0density();
>> 225 X1den = material->GetIonisation()->GetX1density();
>> 226 taul = material->GetIonisation()->GetTaul() ;
>> 227 ShellCorrectionVector = material->GetIonisation()->
>> 228 GetShellCorrectionVector();
>> 229
>> 230 // get elements in the actual material,
>> 231 // they are needed for the low energy part ....
>> 232
>> 233 const G4ElementVector* theElementVector=
>> 234 material->GetElementVector() ;
>> 235 const G4double* theAtomicNumDensityVector=
>> 236 material->GetAtomicNumDensityVector() ;
>> 237 const G4int NumberOfElements=
>> 238 material->GetNumberOfElements() ;
>> 239
>> 240 // get electron cut in kin. energy for the material
>> 241 DeltaCutInKineticEnergyNow =
>> 242 DeltaCutInKineticEnergy[material->GetIndex()] ;
>> 243
>> 244 // some local variables -------------------
>> 245 G4double tau,tau0,Tmax,gamma,bg2,beta2,rcut,delta,x,sh ;
>> 246
>> 247 KineticEnergy=aParticle->GetKineticEnergy();
>> 248
>> 249
>> 250 tau = KineticEnergy/proton_mass_c2 ;
>> 251
>> 252 if ( tau < taul )
>> 253 // low energy part , parametrized energy loss formulae
>> 254 {
>> 255 ionloss = 0. ;
>> 256 // loop for the elements in the material
>> 257 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 258 {
>> 259 const G4Element* element = (*theElementVector)(iel);
>> 260
>> 261 if ( tau < element->GetIonisation()->GetTau0())
>> 262 ionloss += theAtomicNumDensityVector[iel]
>> 263 *( element->GetIonisation()->GetAlow()*sqrt(tau)
>> 264 +element->GetIonisation()->GetBlow()*tau) ;
>> 265 else
>> 266 ionloss += theAtomicNumDensityVector[iel]
>> 267 * element->GetIonisation()->GetClow()/sqrt(tau) ;
>> 268 }
145 } 269 }
146 } << 270 else
147 // to compute ranges correctly we have to << 271 // high energy part , Bethe-Bloch formula
148 // model even if activation limit is high << 272 {
149 EmModel(0)->SetLowEnergyLimit(emin); << 273 gamma = tau +1. ;
150 << 274 bg2 = tau*(tau+2.) ;
151 // high energy limit may be eth or DBL_MAX << 275 beta2 = bg2/(gamma*gamma) ;
152 G4double emax1 = (EmModel(0)->HighEnergyLi << 276 Tmax = 2.*electron_mass_c2*bg2
153 EmModel(0)->SetHighEnergyLimit(emax1); << 277 /(1.+2.*gamma*RateMass+RateMass*RateMass) ;
154 AddEmModel(1, EmModel(0), FluctModel()); << 278
155 << 279 if ( DeltaCutInKineticEnergyNow < Tmax)
156 // second model is used if the first does << 280 rcut = DeltaCutInKineticEnergyNow/Tmax ;
157 if(emax1 < emax) { << 281 else
158 if (nullptr == EmModel(1)) { SetEmModel( << 282 rcut = 1.;
159 EmModel(1)->SetLowEnergyLimit(emax1); << 283
160 << 284 ionloss = log(2.*electron_mass_c2*bg2*Tmax/Eexc2)
161 // for extremely heavy particles upper l << 285
162 // should be increased << 286 +log(rcut)-(1.+rcut)*beta2 ;
163 emax = std::max(emax, eth*10); << 287
164 EmModel(1)->SetHighEnergyLimit(emax); << 288 // density correction
165 AddEmModel(2, EmModel(1), FluctModel()); << 289
166 } << 290 x = log(bg2)/twoln10 ;
167 isInitialised = true; << 291 if ( x < X0den )
>> 292 delta = 0. ;
>> 293 else
>> 294 {
>> 295 delta = twoln10*x - Cden ;
>> 296 if ( x < X1den )
>> 297 delta += Aden*pow((X1den-x),Mden) ;
>> 298 }
>> 299
>> 300 // shell correction
>> 301
>> 302 if ( bg2 > bg2lim ) {
>> 303 sh = 0. ;
>> 304 x = 1. ;
>> 305 for (G4int k=0; k<=2; k++) {
>> 306 x *= bg2 ;
>> 307 sh += ShellCorrectionVector[k]/x;
>> 308 }
>> 309 }
>> 310 else {
>> 311 sh = 0. ;
>> 312 x = 1. ;
>> 313 for (G4int k=0; k<=2; k++) {
>> 314 x *= bg2lim ;
>> 315 sh += ShellCorrectionVector[k]/x;
>> 316 }
>> 317 sh *= log(tau/taul)/log(taulim/taul) ;
>> 318 }
>> 319
>> 320 // now you can compute the total ionization loss
>> 321
>> 322 ionloss -= delta + sh ;
>> 323 ionloss *= Factor*ElectronDensity/beta2 ;
>> 324 }
>> 325 if ( ionloss <= 0.)
>> 326 ionloss = 0. ;
>> 327
>> 328 dEdx = ionloss ;
>> 329 return dEdx ;
>> 330 }
>> 331
>> 332
>> 333 G4VParticleChange* G4ionIonisation::PostStepDoIt(
>> 334 const G4Track& trackData,
>> 335 const G4Step& stepData)
>> 336 {
>> 337 const G4DynamicParticle* aParticle ;
>> 338 G4Material* aMaterial;
>> 339 G4double KineticEnergy,TotalEnergy,TotalMomentum,
>> 340 betasquare,MaxKineticEnergyTransfer,
>> 341 DeltaKineticEnergy,DeltaTotalMomentum,costheta,sintheta,phi,
>> 342 dirx,diry,dirz,finalKineticEnergy,finalPx,finalPy,finalPz,
>> 343 x,xc,grej,Psquare,Esquare,summass,rate,grejc,finalMomentum ;
>> 344
>> 345 aParticleChange.Initialize(trackData) ;
>> 346 aMaterial = trackData.GetMaterial() ;
>> 347
>> 348 aParticle = trackData.GetDynamicParticle() ;
>> 349
>> 350 KineticEnergy=aParticle->GetKineticEnergy();
>> 351 TotalEnergy=KineticEnergy + ParticleMass ;
>> 352 Psquare=KineticEnergy*(TotalEnergy+ParticleMass) ;
>> 353 Esquare=TotalEnergy*TotalEnergy ;
>> 354 summass = ParticleMass + electron_mass_c2 ;
>> 355 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection() ;
>> 356
>> 357 DeltaCutInKineticEnergyNow = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
>> 358
>> 359 // some kinematics......................
>> 360
>> 361 betasquare=Psquare/Esquare ;
>> 362 MaxKineticEnergyTransfer = 2.*electron_mass_c2*Psquare
>> 363 /(summass*summass+2.*electron_mass_c2*KineticEnergy);
>> 364
>> 365 // sampling kinetic energy of the delta ray
>> 366
>> 367 if( MaxKineticEnergyTransfer <= DeltaCutInKineticEnergyNow )
>> 368 {
>> 369 // there is no change at all).....
>> 370 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
168 } 371 }
>> 372 else
>> 373 {
>> 374 // normal case ......................................
>> 375 xc=DeltaCutInKineticEnergyNow/MaxKineticEnergyTransfer ;
>> 376 rate=MaxKineticEnergyTransfer/TotalEnergy ;
>> 377
>> 378 // sampling follows ...
>> 379 grejc=1.-betasquare*xc ;
>> 380
>> 381 do {
>> 382 x=xc/(1.-(1.-xc)*G4UniformRand());
>> 383 grej=(1.-x*betasquare)/grejc ;
>> 384 } while( G4UniformRand()>grej );
>> 385 }
>> 386
>> 387 DeltaKineticEnergy = x * MaxKineticEnergyTransfer ;
>> 388
>> 389 if(DeltaKineticEnergy <= 0.)
>> 390 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 391
>> 392 DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
>> 393 2. * electron_mass_c2 )) ;
>> 394 TotalMomentum = sqrt(Psquare) ;
>> 395 costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
>> 396 /(DeltaTotalMomentum * TotalMomentum) ;
>> 397
>> 398 // protection against costheta > 1 or < -1 ---------------
>> 399 if ( costheta < -1. )
>> 400 costheta = -1. ;
>> 401 if ( costheta > +1. )
>> 402 costheta = +1. ;
>> 403
>> 404 // direction of the delta electron ........
>> 405 phi = twopi * G4UniformRand() ;
>> 406 sintheta = sqrt((1.+costheta)*(1.-costheta));
>> 407 dirx = sintheta * cos(phi) ;
>> 408 diry = sintheta * sin(phi) ;
>> 409 dirz = costheta ;
>> 410
>> 411 G4ThreeVector DeltaDirection(dirx,diry,dirz) ;
>> 412 DeltaDirection.rotateUz(ParticleDirection) ;
>> 413
>> 414 // create G4DynamicParticle object for delta ray
>> 415 G4DynamicParticle *theDeltaRay = new G4DynamicParticle;
>> 416 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
>> 417 theDeltaRay->SetMomentumDirection(
>> 418 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z());
>> 419 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 420
>> 421 // fill aParticleChange
>> 422 finalKineticEnergy = KineticEnergy - DeltaKineticEnergy ;
>> 423 if (finalKineticEnergy > 0.)
>> 424 {
>> 425 // changed energy and momentum of the actual particle
>> 426 finalMomentum=sqrt(finalKineticEnergy*
>> 427 (finalKineticEnergy+2.*ParticleMass)) ;
>> 428
>> 429 finalPx = (TotalMomentum*ParticleDirection.x()
>> 430 -DeltaTotalMomentum*DeltaDirection.x())/finalMomentum ;
>> 431 finalPy = (TotalMomentum*ParticleDirection.y()
>> 432 -DeltaTotalMomentum*DeltaDirection.y())/finalMomentum ;
>> 433 finalPz = (TotalMomentum*ParticleDirection.z()
>> 434 -DeltaTotalMomentum*DeltaDirection.z())/finalMomentum ;
>> 435
>> 436 aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
>> 437 }
>> 438 else
>> 439 {
>> 440 finalKineticEnergy = 0. ;
>> 441 aParticleChange.SetStatusChange(fStopAndKill);
>> 442 }
>> 443
>> 444 aParticleChange.SetEnergyChange( finalKineticEnergy );
>> 445 aParticleChange.SetNumberOfSecondaries(1);
>> 446 aParticleChange.AddSecondary( theDeltaRay );
>> 447
>> 448 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
169 } 449 }
170 450
171 //....oooOO0OOooo........oooOO0OOooo........oo << 451 void G4ionIonisation::PrintInfoDefinition()
172 <<
173 void G4ionIonisation::ProcessDescription(std:: <<
174 { 452 {
175 out << " Ion ionisation"; << 453 G4String comments = " Knock-on electron cross sections . ";
176 G4VEnergyLossProcess::ProcessDescription(out << 454 comments += "\n MeanFreePath is computed at tracking time.\n";
>> 455 comments += " delta ray energy sampled from differential Xsection.";
>> 456
>> 457 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl;
177 } 458 }
178 459
179 //....oooOO0OOooo........oooOO0OOooo........oo <<
180 460