| 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: G4eIonisation.cc,v 1.11.2.2 2001/06/28 20:19:50 gunter Exp $
>> 25 // GEANT4 tag $Name: $
29 // 26 //
>> 27 //
>> 28 // -------------------------------------------------------------
>> 29 // GEANT 4 class implementation file
30 // 30 //
31 // File name: G4eIonisation << 31 // History: based on object model of
>> 32 // 2nd December 1995, G.Cosmo
>> 33 // ---------- G4eIonisation physics process -----------
>> 34 // by Laszlo Urban, 20 March 1997
>> 35 // **************************************************************
>> 36 // It is the first implementation of the NEW IONISATION PROCESS.
>> 37 // It calculates the ionisation of e+/e-.
>> 38 // **************************************************************
32 // 39 //
33 // Author: Laszlo Urban << 40 // 07-04-98: remove 'tracking cut' of the ionizing particle, MMa
34 // << 41 // 04-09-98: new methods SetBining() PrintInfo()
35 // Creation date: 20.03.1997 << 42 // 07-09-98: Cleanup
36 // << 43 // 02/02/99: correction inDoIt , L.Urban
37 // Modified by Michel Maire and Vladimir Ivanc << 44 // 10/02/00 modifications , new e.m. structure, L.Urban
38 // << 45 // 28/05/01 V.Ivanchenko minor changes to provide ANSI -wall compilation
39 // ------------------------------------------- << 46 // --------------------------------------------------------------
40 // << 47
41 //....oooOO0OOooo........oooOO0OOooo........oo <<
42 //....oooOO0OOooo........oooOO0OOooo........oo <<
43 <<
44 #include "G4eIonisation.hh" 48 #include "G4eIonisation.hh"
45 #include "G4Electron.hh" << 49 #include "G4EnergyLossTables.hh"
46 #include "G4MollerBhabhaModel.hh" << 50 #include "G4ios.hh"
47 #include "G4EmParameters.hh" << 51 #include "G4UnitsTable.hh"
48 #include "G4EmStandUtil.hh" << 52
>> 53 G4double G4eIonisation::LowerBoundLambda = 1.*keV ;
>> 54 G4double G4eIonisation::UpperBoundLambda = 100.*TeV ;
>> 55 G4int G4eIonisation::NbinLambda = 100 ;
>> 56
>> 57 // constructor and destructor
>> 58
>> 59 G4eIonisation::G4eIonisation(const G4String& processName)
>> 60 : G4VeEnergyLoss(processName),
>> 61 theMeanFreePathTable(NULL)
>> 62 { }
49 63
50 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 65
52 G4eIonisation::G4eIonisation(const G4String& n << 66 G4eIonisation::~G4eIonisation()
53 : G4VEnergyLossProcess(name), <<
54 theElectron(G4Electron::Electron()), <<
55 isElectron(true), <<
56 isInitialised(false) <<
57 { 67 {
58 SetProcessSubType(fIonisation); << 68 if (theMeanFreePathTable) {
59 SetSecondaryParticle(theElectron); << 69 theMeanFreePathTable->clearAndDestroy();
>> 70 delete theMeanFreePathTable;
>> 71 }
>> 72
60 } 73 }
61 74
62 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 76
64 G4eIonisation::~G4eIonisation() = default; << 77 void G4eIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
>> 78 // just call BuildLossTable+BuildLambdaTable
>> 79 {
>> 80 // get bining from EnergyLoss
>> 81 LowestKineticEnergy = GetLowerBoundEloss() ;
>> 82 HighestKineticEnergy = GetUpperBoundEloss() ;
>> 83 TotBin = GetNbinEloss() ;
>> 84
>> 85 BuildLossTable(aParticleType) ;
>> 86
>> 87 if(&aParticleType==G4Electron::Electron())
>> 88 {
>> 89 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable ;
>> 90 CounterOfElectronProcess++;
>> 91 }
>> 92 else
>> 93 {
>> 94 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable ;
>> 95 CounterOfPositronProcess++;
>> 96 }
>> 97
>> 98 BuildLambdaTable(aParticleType) ;
>> 99
>> 100 BuildDEDXTable(aParticleType) ;
>> 101
>> 102 if(&aParticleType==G4Electron::Electron())
>> 103 PrintInfoDefinition();
>> 104 }
65 105
66 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 107
68 G4double G4eIonisation::MinPrimaryEnergy(const << 108 void G4eIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
69 const G4Material*, <<
70 G4double cut) <<
71 { 109 {
72 G4double x = cut; << 110 // Build tables for the ionization energy loss
73 if(isElectron) x += cut; << 111 // the tables are built for *MATERIALS*
74 return x; << 112
>> 113 const G4double twoln10 = 2.*log(10.);
>> 114 const G4double Factor = twopi_mc2_rcl2;
>> 115
>> 116 static const G4double Tl = 0.2*keV ;
>> 117
>> 118 G4double LowEdgeEnergy, ionloss;
>> 119
>> 120 // material properties
>> 121 G4double ElectronDensity,Eexc,Eexcm2,Cden,Mden,Aden,X0den,X1den ;
>> 122 // some local variables
>> 123 G4double tau,Tmax,gamma,gamma2,bg2,beta2,d,d2,d3,d4,delta,x,y ;
>> 124
>> 125 ParticleMass = aParticleType.GetPDGMass();
>> 126 G4double* ParticleCutInKineticEnergy = aParticleType.GetEnergyCuts() ;
>> 127
>> 128 // create table
>> 129
>> 130 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 131 G4int numOfMaterials = theMaterialTable->length();
>> 132
>> 133 if (theLossTable) { theLossTable->clearAndDestroy();
>> 134 delete theLossTable;
>> 135 }
>> 136
>> 137 theLossTable = new G4PhysicsTable(numOfMaterials);
>> 138
>> 139 // loop for materials
>> 140
>> 141 for (G4int J=0; J<numOfMaterials; J++)
>> 142 {
>> 143 // create physics vector and fill it
>> 144
>> 145 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 146 LowestKineticEnergy, HighestKineticEnergy, TotBin);
>> 147
>> 148 // get material parameters needed for the energy loss calculation
>> 149 const G4Material* material= (*theMaterialTable)[J];
>> 150
>> 151 ElectronDensity = material->GetElectronDensity();
>> 152 Eexc = material->GetIonisation()->GetMeanExcitationEnergy();
>> 153 Eexc /= ParticleMass; Eexcm2 = Eexc*Eexc;
>> 154 Cden = material->GetIonisation()->GetCdensity();
>> 155 Mden = material->GetIonisation()->GetMdensity();
>> 156 Aden = material->GetIonisation()->GetAdensity();
>> 157 X0den = material->GetIonisation()->GetX0density();
>> 158 X1den = material->GetIonisation()->GetX1density();
>> 159
>> 160 // for the lowenergy extrapolation
>> 161 G4double Zeff = material->GetTotNbOfElectPerVolume()/
>> 162 material->GetTotNbOfAtomsPerVolume() ;
>> 163 G4double Th = 0.25*sqrt(Zeff)*keV ;
>> 164
>> 165 G4double Tsav ;
>> 166
>> 167 // now comes the loop for the kinetic energy values
>> 168
>> 169 for (G4int i = 0 ; i < TotBin ; i++)
>> 170 {
>> 171 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 172 // low energy ?
>> 173 if(LowEdgeEnergy < Th)
>> 174 {
>> 175 Tsav = LowEdgeEnergy ;
>> 176 LowEdgeEnergy = Th ;
>> 177 }
>> 178 else
>> 179 Tsav = 0. ;
>> 180
>> 181 tau = LowEdgeEnergy/ParticleMass ;
>> 182
>> 183 // Seltzer-Berger formula
>> 184 gamma = tau + 1.; gamma2 = gamma*gamma;
>> 185 bg2 = tau*(tau+2.);
>> 186 beta2 = bg2/gamma2;
>> 187
>> 188 // electron
>> 189 if (&aParticleType==G4Electron::Electron())
>> 190 {
>> 191 Tmax = LowEdgeEnergy/2.;
>> 192 d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
>> 193 ionloss = log(2.*(tau+2.)/Eexcm2)-1.-beta2
>> 194 + log((tau-d)*d)+tau/(tau-d)
>> 195 + (0.5*d*d+(2.*tau+1.)*log(1.-d/tau))/gamma2;
>> 196 }
>> 197 else //positron
>> 198 {
>> 199 Tmax = LowEdgeEnergy ;
>> 200 d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass;
>> 201 d2=d*d/2.; d3=d*d*d/3.; d4=d*d*d*d/4.;
>> 202 y=1./(1.+gamma);
>> 203 ionloss = log(2.*(tau+2.)/Eexcm2)+log(tau*d)
>> 204 - beta2*(tau+2.*d-y*(3.*d2+y*(d-d3+y*(d2-tau*d3+d4))))/tau;
>> 205 }
>> 206
>> 207 //density correction
>> 208 x = log(bg2)/twoln10;
>> 209 if (x < X0den) delta = 0.;
>> 210 else { delta = twoln10*x - Cden;
>> 211 if (x < X1den) delta += Aden*pow((X1den-x),Mden);
>> 212 }
>> 213
>> 214 //now you can compute the total ionization loss
>> 215 ionloss -= delta ;
>> 216 ionloss *= Factor*ElectronDensity/beta2 ;
>> 217 if (ionloss <= 0.) ionloss = 0.;
>> 218
>> 219 // low energy ?
>> 220 if(Tsav > 0.)
>> 221 {
>> 222 if(Tsav >= Tl)
>> 223 ionloss *= sqrt(LowEdgeEnergy/Tsav) ;
>> 224 else
>> 225 ionloss *= sqrt(LowEdgeEnergy*Tsav)/Tl ;
>> 226 }
>> 227
>> 228 aVector->PutValue(i,ionloss) ;
>> 229 }
>> 230 theLossTable->insert(aVector);
>> 231 }
75 } 232 }
76 233
77 //....oooOO0OOooo........oooOO0OOooo........oo 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 235
79 G4bool G4eIonisation::IsApplicable(const G4Par << 236 void G4eIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
80 { 237 {
81 return (&p == theElectron || &p == G4Positro << 238 // Build mean free path tables for the delta ray production process
>> 239 // tables are built for MATERIALS
>> 240
>> 241 G4double LowEdgeEnergy, Value, SIGMA;
>> 242
>> 243 //create table
>> 244 const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
>> 245 G4int numOfMaterials = theMaterialTable->length();
>> 246
>> 247 if (theMeanFreePathTable) { theMeanFreePathTable->clearAndDestroy();
>> 248 delete theMeanFreePathTable;
>> 249 }
>> 250
>> 251 theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
>> 252
>> 253 // get electron cuts in kinetic energy
>> 254 // The electron cuts needed in the case of the positron , too!
>> 255 // This is the reason why SetCut has to be called for electron first !!!!!!!
>> 256
>> 257 if((G4Electron::Electron()->GetCutsInEnergy() == 0) &&
>> 258 ( &aParticleType == G4Positron::Positron()))
>> 259 {
>> 260 G4cout << " The ELECTRON energy cuts needed to compute energy loss/mean free path "
>> 261 " for POSITRON , too. " << G4endl;
>> 262 G4Exception(" Call SetCut for e- first !!!!!!") ;
>> 263 }
>> 264
>> 265 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
>> 266
>> 267 // loop for materials
>> 268
>> 269 for (G4int J=0 ; J < numOfMaterials; J++)
>> 270 {
>> 271 //create physics vector then fill it ....
>> 272
>> 273 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 274 LowerBoundLambda, UpperBoundLambda, NbinLambda);
>> 275
>> 276 // compute the (macroscopic) cross section first
>> 277
>> 278 const G4Material* material= (*theMaterialTable)[J];
>> 279 const
>> 280 G4ElementVector* theElementVector = material->GetElementVector();
>> 281 const
>> 282 G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
>> 283 const
>> 284 G4int NumberOfElements = material->GetNumberOfElements() ;
>> 285
>> 286 // get the electron kinetic energy cut for the actual material,
>> 287 // it will be used in ComputeMicroscopicCrossSection
>> 288 // (--> it will be the same for all the elements in this material )
>> 289 G4double DeltaThreshold = DeltaCutInKineticEnergy[J] ;
>> 290
>> 291 for (G4int i = 0 ; i < NbinLambda ; i++)
>> 292 {
>> 293 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 294 SIGMA = 0.;
>> 295 for (G4int iel=0; iel<NumberOfElements; iel++ )
>> 296 {
>> 297 SIGMA += theAtomicNumDensityVector[iel]*
>> 298 ComputeMicroscopicCrossSection( aParticleType,
>> 299 LowEdgeEnergy,
>> 300 (*theElementVector)(iel)->GetZ(),
>> 301 DeltaThreshold);
>> 302 }
>> 303
>> 304 // mean free path = 1./macroscopic cross section
>> 305 Value = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
>> 306 aVector->PutValue(i, Value) ;
>> 307 }
>> 308 theMeanFreePathTable->insert(aVector);
>> 309 }
82 } 310 }
83 311
84 //....oooOO0OOooo........oooOO0OOooo........oo 312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 313
86 void G4eIonisation::InitialiseEnergyLossProces << 314 G4double G4eIonisation::ComputeMicroscopicCrossSection(
87 const G4ParticleDefinition* part, << 315 const G4ParticleDefinition& aParticleType,
88 const G4ParticleDefinition*) << 316 G4double KineticEnergy,
>> 317 G4double AtomicNumber ,
>> 318 G4double DeltaThreshold)
89 { 319 {
90 if(!isInitialised) { << 320 // calculates the microscopic cross section
91 if(part != theElectron) { isElectron = fal << 321 //(it is called for elements , AtomicNumber = Z )
92 if (nullptr == EmModel(0)) { SetEmModel(ne << 322
93 G4EmParameters* param = G4EmParameters::In << 323 G4double MaxKineticEnergyTransfer, TotalCrossSection(0.);
94 EmModel(0)->SetLowEnergyLimit(param->MinKi << 324
95 EmModel(0)->SetHighEnergyLimit(param->MaxK << 325 ParticleMass = aParticleType.GetPDGMass();
96 if (nullptr == FluctModel()) { << 326 G4double TotalEnergy = KineticEnergy + ParticleMass;
97 SetFluctModel(G4EmStandUtil::ModelOfFluc << 327
>> 328 G4double betasquare = KineticEnergy*(TotalEnergy+ParticleMass)
>> 329 /(TotalEnergy*TotalEnergy);
>> 330 G4double gamma = TotalEnergy/ParticleMass, gamma2 = gamma*gamma;
>> 331 G4double x=DeltaThreshold/KineticEnergy, x2 = x*x;
>> 332
>> 333 if (&aParticleType==G4Electron::Electron())
>> 334 MaxKineticEnergyTransfer = 0.5*KineticEnergy;
>> 335 else MaxKineticEnergyTransfer = KineticEnergy;
>> 336
>> 337 // now you can calculate the total cross section
>> 338
>> 339 if (MaxKineticEnergyTransfer > DeltaThreshold)
>> 340 {
>> 341 if (&aParticleType==G4Electron::Electron()) //Moller (e-e-) scattering
>> 342 {
>> 343 TotalCrossSection = (gamma-1.)*(gamma-1.)*(0.5-x)/gamma2 + 1./x
>> 344 - 1./(1.-x)-(2.*gamma-1.)*log((1.-x)/x)/gamma2;
>> 345 TotalCrossSection /= betasquare;
>> 346 }
>> 347 else //Bhabha (e+e-) scattering
>> 348 {
>> 349 G4double y=1./(1.+gamma), y2 =y*y, y12=1.-2.*y;
>> 350 G4double b1=2.-y2, b2=y12*(3.+y2), b4=y12*y12*y12, b3=b4+y12*y12;
>> 351 TotalCrossSection = (1./x-1.)/betasquare+b1*log(x)+b2*(1.-x)
>> 352 - b3*(1.-x2)/2.+b4*(1.-x2*x)/3.;
>> 353 }
>> 354 TotalCrossSection *= (twopi_mc2_rcl2*AtomicNumber/KineticEnergy);
>> 355 }
>> 356 return TotalCrossSection ;
>> 357 }
>> 358
>> 359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 360
>> 361 G4VParticleChange* G4eIonisation::PostStepDoIt( const G4Track& trackData,
>> 362 const G4Step& stepData)
>> 363 {
>> 364 aParticleChange.Initialize(trackData) ;
>> 365
>> 366 G4Material* aMaterial = trackData.GetMaterial() ;
>> 367 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle() ;
>> 368
>> 369 ParticleMass = aParticle->GetDefinition()->GetPDGMass();
>> 370 G4double KineticEnergy = aParticle->GetKineticEnergy();
>> 371 G4double TotalEnergy = KineticEnergy + ParticleMass;
>> 372 G4double Psquare = KineticEnergy*(TotalEnergy+ParticleMass);
>> 373 G4double TotalMomentum = sqrt(Psquare);
>> 374 //G4double Esquare=TotalEnergy*TotalEnergy;
>> 375 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
>> 376
>> 377 // get kinetic energy cut for the electron
>> 378 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ;
>> 379 G4double DeltaThreshold = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
>> 380
>> 381 // some kinematics
>> 382 G4double MaxKineticEnergyTransfer;
>> 383 if (Charge < 0.) MaxKineticEnergyTransfer = 0.5*KineticEnergy;
>> 384 else MaxKineticEnergyTransfer = KineticEnergy;
>> 385
>> 386 // sampling kinetic energy of the delta ray
>> 387
>> 388 if (MaxKineticEnergyTransfer <= DeltaThreshold) // pathological case (should not happen,
>> 389 // there is no change at all)
>> 390 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 391
>> 392
>> 393 // normal case
>> 394 G4double cc,y,y2,c2,b0,b1,b2,b3,b4,x,x1,grej,grejc;
>> 395
>> 396 G4double tau = KineticEnergy/ParticleMass;
>> 397 G4double gamma = tau+1., gamma2=gamma*gamma;
>> 398 G4double xc = DeltaThreshold/KineticEnergy, xc1=1.-xc;
>> 399
>> 400 if (Charge < 0.) // Moller (e-e-) scattering
>> 401 {
>> 402 b1=4./(9.*gamma2-10.*gamma+5.);
>> 403 b2=tau*tau*b1; b3=(2.*gamma2+2.*gamma-1.)*b1;
>> 404 cc=1.-2.*xc;
>> 405 do {
>> 406 x = xc/(1.-cc*G4UniformRand()); x1 = 1.-x;
>> 407 grej = b2*x*x-b3*x/x1+b1*gamma2/(x1*x1);
>> 408 } while (G4UniformRand()>grej) ;
98 } 409 }
99 AddEmModel(1, EmModel(), FluctModel()); << 410 else // Bhabha (e+e-) scattering
100 isInitialised = true; << 411 {
101 } << 412 y=1./(gamma+1.); y2=y*y; cc=1.-2.*y;
>> 413 b1=2.-y2; b2=cc*(3.+y2);
>> 414 c2=cc*cc; b4=c2*cc; b3=c2+b4;
>> 415 b0=gamma2/(gamma2-1.);
>> 416 grejc=(((b4*xc-b3)*xc+b2)*xc-b1)*xc+b0;
>> 417 do {
>> 418 x = xc/(1.-xc1*G4UniformRand());
>> 419 grej = ((((b4*x-b3)*x+b2)*x-b1)*x+b0)/grejc;
>> 420 } while (G4UniformRand()>grej);
>> 421 }
>> 422
>> 423 G4double DeltaKineticEnergy = x * KineticEnergy;
>> 424
>> 425 // protection :do not produce a secondary with 0. kinetic energy !
>> 426 if (DeltaKineticEnergy <= 0.)
>> 427 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 428
>> 429 G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
>> 430 2. * electron_mass_c2 ));
>> 431
>> 432 G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
>> 433 /(DeltaTotalMomentum * TotalMomentum);
>> 434
>> 435 if (costheta < -1.) costheta = -1.;
>> 436 if (costheta > +1.) costheta = +1.;
>> 437
>> 438 // direction of the delta electron
>> 439
>> 440 G4double phi = twopi * G4UniformRand();
>> 441 G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
>> 442 G4double dirx = sintheta * cos(phi), diry = sintheta * sin(phi), dirz = costheta;
>> 443
>> 444 G4ThreeVector DeltaDirection(dirx,diry,dirz);
>> 445 DeltaDirection.rotateUz(ParticleDirection);
>> 446
>> 447 // create G4DynamicParticle object for delta ray
>> 448
>> 449 G4DynamicParticle* theDeltaRay = new G4DynamicParticle;
>> 450 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
>> 451 theDeltaRay->SetMomentumDirection(
>> 452 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z());
>> 453 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 454
>> 455 // fill aParticleChange
>> 456 // changed energy and momentum of the actual particle
>> 457 G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy;
>> 458
>> 459 G4double Edep = 0. ;
>> 460
>> 461 if (finalKineticEnergy > MinKineticEnergy)
>> 462 {
>> 463 G4double finalPx = TotalMomentum*ParticleDirection.x()
>> 464 - DeltaTotalMomentum*DeltaDirection.x();
>> 465 G4double finalPy = TotalMomentum*ParticleDirection.y()
>> 466 - DeltaTotalMomentum*DeltaDirection.y();
>> 467 G4double finalPz = TotalMomentum*ParticleDirection.z()
>> 468 - DeltaTotalMomentum*DeltaDirection.z();
>> 469 G4double finalMomentum =
>> 470 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
>> 471 finalPx /= finalMomentum ;
>> 472 finalPy /= finalMomentum ;
>> 473 finalPz /= finalMomentum ;
>> 474
>> 475 aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
>> 476 }
>> 477 else
>> 478 {
>> 479 finalKineticEnergy = 0.;
>> 480 Edep = finalKineticEnergy ;
>> 481 if (Charge < 0.) aParticleChange.SetStatusChange(fStopAndKill);
>> 482 else aParticleChange.SetStatusChange(fStopButAlive);
>> 483 }
>> 484
>> 485 aParticleChange.SetEnergyChange( finalKineticEnergy );
>> 486 aParticleChange.SetNumberOfSecondaries(1);
>> 487 aParticleChange.AddSecondary( theDeltaRay );
>> 488 aParticleChange.SetLocalEnergyDeposit (Edep);
>> 489
>> 490 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
102 } 491 }
103 492
104 //....oooOO0OOooo........oooOO0OOooo........oo 493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 494
106 void G4eIonisation::ProcessDescription(std::os << 495 void G4eIonisation::PrintInfoDefinition()
107 { 496 {
108 out << " Ionisation"; << 497 G4String comments = "delta cross sections from Moller+Bhabha. ";
109 G4VEnergyLossProcess::ProcessDescription(out << 498 comments += "Good description from 1 KeV to 100 GeV.\n";
110 } << 499 comments += " delta ray energy sampled from differential Xsection.";
>> 500
>> 501 G4cout << G4endl << GetProcessName() << ": " << comments
>> 502 << "\n PhysicsTables from " << G4BestUnit(LowerBoundLambda,"Energy")
>> 503 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 504 << " in " << NbinLambda << " bins. \n";
>> 505 }
111 506
112 //....oooOO0OOooo........oooOO0OOooo........oo << 507 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 508