| 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.25 2002/04/09 17:34:44 vnivanch Exp $
>> 25 // GEANT4 tag $Name: geant4-04-01 $
29 // 26 //
>> 27 //--------------- G4eIonisation physics process --------------------------------
>> 28 // by Laszlo Urban, 20 March 1997
>> 29 //------------------------------------------------------------------------------
30 // 30 //
31 // File name: G4eIonisation << 31 // 07-04-98 remove 'tracking cut' of the ionizing particle, mma
32 // << 32 // 04-09-98 new methods SetBining() PrintInfo()
33 // Author: Laszlo Urban << 33 // 07-09-98 Cleanup
34 // << 34 // 02-02-99 correction inDoIt , L.Urban
35 // Creation date: 20.03.1997 << 35 // 10-02-00 modifications , new e.m. structure, L.Urban
36 // << 36 // 28-05-01 V.Ivanchenko minor changes to provide ANSI -wall compilation
37 // Modified by Michel Maire and Vladimir Ivanc << 37 // 09-08-01 new methods Store/Retrieve PhysicsTable (mma)
38 // << 38 // 13-08-01 new function ComputeRestrictedMeandEdx() (mma)
39 // ------------------------------------------- << 39 // 17-09-01 migration of Materials to pure STL (mma)
40 // << 40 // 21-09-01 completion of RetrievePhysicsTable() (mma)
41 //....oooOO0OOooo........oooOO0OOooo........oo << 41 // 29-10-01 all static functions no more inlined (mma)
42 //....oooOO0OOooo........oooOO0OOooo........oo << 42 // 07-11-01 particleMass and Charge become local variables
43 << 43 // 26-03-02 change access to cuts in BuildLossTables (V.Ivanchenko)
>> 44 //------------------------------------------------------------------------------
>> 45
>> 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 48
44 #include "G4eIonisation.hh" 49 #include "G4eIonisation.hh"
45 #include "G4Electron.hh" << 50 #include "G4UnitsTable.hh"
46 #include "G4MollerBhabhaModel.hh" <<
47 #include "G4EmParameters.hh" <<
48 #include "G4EmStandUtil.hh" <<
49 51
50 //....oooOO0OOooo........oooOO0OOooo........oo << 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 53
52 G4eIonisation::G4eIonisation(const G4String& n << 54 G4double G4eIonisation::LowerBoundLambda = 1.*keV;
53 : G4VEnergyLossProcess(name), << 55 G4double G4eIonisation::UpperBoundLambda = 100.*TeV;
54 theElectron(G4Electron::Electron()), << 56 G4int G4eIonisation::NbinLambda = 100;
55 isElectron(true), << 57
56 isInitialised(false) << 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 59
>> 60 G4eIonisation::G4eIonisation(const G4String& processName)
>> 61 : G4VeEnergyLoss(processName),
>> 62 theMeanFreePathTable(NULL)
>> 63 {}
>> 64
>> 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 66
>> 67 G4eIonisation::~G4eIonisation()
57 { 68 {
58 SetProcessSubType(fIonisation); << 69 if (theMeanFreePathTable)
59 SetSecondaryParticle(theElectron); << 70 {theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
60 } 71 }
61 72
62 //....oooOO0OOooo........oooOO0OOooo........oo << 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 74
64 G4eIonisation::~G4eIonisation() = default; << 75 void G4eIonisation::SetLowerBoundLambda(G4double val)
>> 76 {LowerBoundLambda = val;}
>> 77
>> 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 79
>> 80 void G4eIonisation::SetUpperBoundLambda(G4double val)
>> 81 {UpperBoundLambda = val;}
>> 82
>> 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 84
>> 85 void G4eIonisation::SetNbinLambda(G4int n)
>> 86 {NbinLambda = n;}
>> 87
>> 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 89
>> 90 G4double G4eIonisation::GetLowerBoundLambda()
>> 91 {return LowerBoundLambda;}
>> 92
>> 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 94
>> 95 G4double G4eIonisation::GetUpperBoundLambda()
>> 96 {return UpperBoundLambda;}
>> 97
>> 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 99
>> 100 G4int G4eIonisation::GetNbinLambda()
>> 101 {return NbinLambda;}
>> 102
>> 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 104
>> 105 void G4eIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
>> 106 // just call BuildLossTable+BuildLambdaTable
>> 107 {
>> 108 // get bining from EnergyLoss
>> 109 LowestKineticEnergy = GetLowerBoundEloss();
>> 110 HighestKineticEnergy = GetUpperBoundEloss();
>> 111 TotBin = GetNbinEloss();
>> 112
>> 113 BuildLossTable(aParticleType);
>> 114
>> 115 if (&aParticleType==G4Electron::Electron())
>> 116 {
>> 117 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 118 CounterOfElectronProcess++;
>> 119 }
>> 120 else
>> 121 {
>> 122 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 123 CounterOfPositronProcess++;
>> 124 }
>> 125
>> 126 BuildLambdaTable(aParticleType);
>> 127
>> 128 BuildDEDXTable(aParticleType);
>> 129
>> 130 if (&aParticleType==G4Electron::Electron()) PrintInfoDefinition();
>> 131 }
65 132
66 //....oooOO0OOooo........oooOO0OOooo........oo << 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 134
68 G4double G4eIonisation::MinPrimaryEnergy(const << 135 void G4eIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
69 const G4Material*, <<
70 G4double cut) <<
71 { 136 {
72 G4double x = cut; << 137 // Build tables of dE/dx due to the ionization process
73 if(isElectron) x += cut; << 138 // the tables are built for *MATERIALS*
74 return x; << 139
>> 140 // create table
>> 141 //
>> 142 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 143 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
>> 144
>> 145 if (theLossTable) {theLossTable->clearAndDestroy(); delete theLossTable;}
>> 146 theLossTable = new G4PhysicsTable(numOfMaterials);
>> 147
>> 148
>> 149 // get electron cuts in kinetic energy
>> 150 // The electron cuts needed in the case of the positron , too!
>> 151 // This is the reason why SetCut has to be called for electron first !!
>> 152
>> 153 if((G4Electron::Electron()->GetEnergyCuts() == 0) &&
>> 154 (&aParticleType == G4Positron::Positron()))
>> 155 {
>> 156 G4cout << " The ELECTRON energy cuts needed to compute energy loss"
>> 157 " and mean free path; and for POSITRON, too. " << G4endl;
>> 158 G4Exception(" Call SetCut for e- first !!");
>> 159 }
>> 160
>> 161 // get DeltaCut in energy
>> 162 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetEnergyCuts();
>> 163
>> 164 // loop for materials
>> 165 //
>> 166 for (G4int J=0; J<numOfMaterials; J++)
>> 167 {
>> 168 // create physics vector and fill it
>> 169 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 170 LowestKineticEnergy, HighestKineticEnergy, TotBin);
>> 171
>> 172 const G4Material* material = (*theMaterialTable)[J];
>> 173 G4double DeltaThreshold = DeltaCutInKineticEnergy[J];
>> 174
>> 175 // now comes the loop for the kinetic energy values
>> 176 //
>> 177 for (G4int i = 0; i < TotBin; i++)
>> 178 {
>> 179 G4double dEdx = ComputeRestrictedMeandEdx(aParticleType,
>> 180 aVector->GetLowEdgeEnergy(i),
>> 181 material,
>> 182 DeltaThreshold);
>> 183 if(1 < verboseLevel) {
>> 184 G4cout << "Material= " << material->GetName()
>> 185 << " E(MeV)= " << aVector->GetLowEdgeEnergy(i)/MeV
>> 186 << " dEdx(MeV/mm)= " << dEdx*mm/MeV
>> 187 << G4endl;
>> 188 }
>> 189 aVector->PutValue(i,dEdx);
>> 190 }
>> 191 theLossTable->insert(aVector);
>> 192 }
75 } 193 }
76 194
77 //....oooOO0OOooo........oooOO0OOooo........oo << 195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 196
79 G4bool G4eIonisation::IsApplicable(const G4Par << 197 void G4eIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
80 { 198 {
81 return (&p == theElectron || &p == G4Positro << 199 // Build mean free path tables for the delta ray production process
>> 200 // tables are built for MATERIALS
>> 201
>> 202 //create table
>> 203 //
>> 204 const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
>> 205 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
>> 206
>> 207 if (theMeanFreePathTable)
>> 208 { theMeanFreePathTable->clearAndDestroy(); delete theMeanFreePathTable;}
>> 209
>> 210 theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
>> 211
>> 212 // get electron cuts in kinetic energy
>> 213 // The electron cuts needed in the case of the positron , too!
>> 214 // This is the reason why SetCut has to be called for electron first !!
>> 215
>> 216 if((G4Electron::Electron()->GetEnergyCuts() == 0) &&
>> 217 (&aParticleType == G4Positron::Positron()))
>> 218 {
>> 219 G4cout << " The ELECTRON energy cuts needed to compute energy loss"
>> 220 " and mean free path; and for POSITRON, too. " << G4endl;
>> 221 G4Exception(" Call SetCut for e- first !!");
>> 222 }
>> 223
>> 224 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetEnergyCuts();
>> 225
>> 226 // loop for materials
>> 227
>> 228 for (G4int J=0 ; J < numOfMaterials; J++)
>> 229 {
>> 230 //create physics vector then fill it ....
>> 231
>> 232 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 233 LowerBoundLambda, UpperBoundLambda, NbinLambda);
>> 234
>> 235 // compute the (macroscopic) cross section first
>> 236
>> 237 const G4Material* material= (*theMaterialTable)[J];
>> 238 const G4ElementVector* theElementVector = material->GetElementVector();
>> 239 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
>> 240 G4int NumberOfElements = material->GetNumberOfElements();
>> 241
>> 242 // get the electron kinetic energy cut for the actual material,
>> 243 // it will be used in ComputeCrossSectionPerAtom
>> 244 // (--> it will be the same for all the elements in this material )
>> 245 G4double DeltaThreshold = DeltaCutInKineticEnergy[J];
>> 246
>> 247 for (G4int i = 0 ; i < NbinLambda ; i++)
>> 248 {
>> 249 G4double LowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
>> 250 G4double SIGMA = 0.;
>> 251 for (G4int iel=0; iel<NumberOfElements; iel++ )
>> 252 {
>> 253 SIGMA += NbOfAtomsPerVolume[iel]*
>> 254 ComputeCrossSectionPerAtom(aParticleType,
>> 255 LowEdgeEnergy,
>> 256 (*theElementVector)[iel]->GetZ(),
>> 257 DeltaThreshold);
>> 258 }
>> 259
>> 260 // mean free path = 1./macroscopic cross section
>> 261 G4double Value = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
>> 262 aVector->PutValue(i, Value);
>> 263 }
>> 264 theMeanFreePathTable->insert(aVector);
>> 265 }
82 } 266 }
83 267
84 //....oooOO0OOooo........oooOO0OOooo........oo << 268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 269
86 void G4eIonisation::InitialiseEnergyLossProces << 270 G4double G4eIonisation::ComputeRestrictedMeandEdx (
87 const G4ParticleDefinition* part, << 271 const G4ParticleDefinition& aParticleType,
88 const G4ParticleDefinition*) << 272 G4double KineticEnergy,
>> 273 const G4Material* material,
>> 274 G4double DeltaThreshold)
89 { 275 {
90 if(!isInitialised) { << 276 // calculate the dE/dx due to the ionization process (Geant4 internal units)
91 if(part != theElectron) { isElectron = fal << 277 // Seltzer-Berger formula
92 if (nullptr == EmModel(0)) { SetEmModel(ne << 278 //
93 G4EmParameters* param = G4EmParameters::In << 279 G4double particleMass = aParticleType.GetPDGMass();
94 EmModel(0)->SetLowEnergyLimit(param->MinKi << 280
95 EmModel(0)->SetHighEnergyLimit(param->MaxK << 281 G4double ElectronDensity = material->GetElectronDensity();
96 if (nullptr == FluctModel()) { << 282 G4double Eexc = material->GetIonisation()->GetMeanExcitationEnergy();
97 SetFluctModel(G4EmStandUtil::ModelOfFluc << 283 Eexc /= particleMass; G4double Eexcm2 = Eexc*Eexc;
>> 284
>> 285 // for the lowenergy extrapolation
>> 286 G4double Zeff = material->GetTotNbOfElectPerVolume()/
>> 287 material->GetTotNbOfAtomsPerVolume();
>> 288 G4double Th = 0.25*sqrt(Zeff)*keV;
>> 289 G4double Tsav = 0.;
>> 290 if (KineticEnergy < Th) {Tsav = KineticEnergy; KineticEnergy = Th;}
>> 291
>> 292 G4double tau = KineticEnergy/particleMass;
>> 293 G4double gamma = tau + 1., gamma2 = gamma*gamma, bg2 = tau*(tau+2.);
>> 294 G4double beta2 = bg2/gamma2;
>> 295
>> 296 G4double Tmax,d,dEdx;
>> 297
>> 298 // electron
>> 299 if (&aParticleType==G4Electron::Electron())
>> 300 {
>> 301 Tmax = KineticEnergy/2.;
>> 302 d = G4std::min(DeltaThreshold, Tmax)/particleMass;
>> 303 dEdx = log(2.*(tau+2.)/Eexcm2)-1.-beta2
>> 304 + log((tau-d)*d)+tau/(tau-d)
>> 305 + (0.5*d*d+(2.*tau+1.)*log(1.-d/tau))/gamma2;
>> 306 }
>> 307
>> 308 else //positron
>> 309 {
>> 310 Tmax = KineticEnergy;
>> 311 d = G4std::min(DeltaThreshold, Tmax)/particleMass;
>> 312 G4double d2=d*d/2., d3=d*d*d/3., d4=d*d*d*d/4.;
>> 313 G4double y=1./(1.+gamma);
>> 314 dEdx = log(2.*(tau+2.)/Eexcm2)+log(tau*d)
>> 315 - beta2*(tau+2.*d-y*(3.*d2+y*(d-d3+y*(d2-tau*d3+d4))))/tau;
>> 316 }
>> 317
>> 318 //density correction
>> 319 G4double Cden = material->GetIonisation()->GetCdensity();
>> 320 G4double Mden = material->GetIonisation()->GetMdensity();
>> 321 G4double Aden = material->GetIonisation()->GetAdensity();
>> 322 G4double X0den = material->GetIonisation()->GetX0density();
>> 323 G4double X1den = material->GetIonisation()->GetX1density();
>> 324
>> 325 const G4double twoln10 = 2.*log(10.);
>> 326 G4double x = log(bg2)/twoln10;
>> 327 G4double delta;
>> 328 if (x < X0den) delta = 0.;
>> 329 else {delta = twoln10*x - Cden;
>> 330 if (x < X1den) delta += Aden*pow((X1den-x),Mden);
>> 331 }
>> 332
>> 333 //now you can compute the total ionization loss
>> 334 dEdx -= delta;
>> 335 dEdx *= twopi_mc2_rcl2*ElectronDensity/beta2;
>> 336 if (dEdx <= 0.) dEdx = 0.;
>> 337
>> 338 // low energy ?
>> 339 const G4double Tl = 0.2*keV;
>> 340 if (Tsav > 0.)
>> 341 {
>> 342 if (Tsav >= Tl) dEdx *= sqrt(KineticEnergy/Tsav);
>> 343 else dEdx *= sqrt(KineticEnergy*Tsav)/Tl;
>> 344 }
>> 345 return dEdx;
>> 346 }
>> 347
>> 348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 349
>> 350 G4double G4eIonisation::ComputeCrossSectionPerAtom(
>> 351 const G4ParticleDefinition& aParticleType,
>> 352 G4double KineticEnergy,
>> 353 G4double AtomicNumber ,
>> 354 G4double DeltaThreshold)
>> 355 {
>> 356 // calculates the cross section per atom (Geant4 internal units)
>> 357 //(it is called for elements , AtomicNumber = Z )
>> 358
>> 359 G4double particleMass = aParticleType.GetPDGMass();
>> 360 G4double TotalEnergy = KineticEnergy + particleMass;
>> 361
>> 362 G4double betasquare = KineticEnergy*(TotalEnergy+particleMass)
>> 363 /(TotalEnergy*TotalEnergy);
>> 364 G4double gamma = TotalEnergy/particleMass, gamma2 = gamma*gamma;
>> 365 G4double x=DeltaThreshold/KineticEnergy, x2 = x*x;
>> 366
>> 367 G4double MaxKineticEnergyTransfer;
>> 368 if (&aParticleType==G4Electron::Electron())
>> 369 MaxKineticEnergyTransfer = 0.5*KineticEnergy;
>> 370 else MaxKineticEnergyTransfer = KineticEnergy;
>> 371
>> 372 // now you can calculate the total cross section
>> 373 //
>> 374 G4double TotalCrossSection = 0.;
>> 375 if (MaxKineticEnergyTransfer > DeltaThreshold)
>> 376 {
>> 377 if (&aParticleType==G4Electron::Electron()) //Moller (e-e-) scattering
>> 378 {
>> 379 TotalCrossSection = (gamma-1.)*(gamma-1.)*(0.5-x)/gamma2 + 1./x
>> 380 - 1./(1.-x)-(2.*gamma-1.)*log((1.-x)/x)/gamma2;
>> 381 TotalCrossSection /= betasquare;
>> 382 }
>> 383 else //Bhabha (e+e-) scattering
>> 384 {
>> 385 G4double y=1./(1.+gamma), y2 =y*y, y12=1.-2.*y;
>> 386 G4double b1=2.-y2, b2=y12*(3.+y2), b4=y12*y12*y12, b3=b4+y12*y12;
>> 387 TotalCrossSection = (1./x-1.)/betasquare+b1*log(x)+b2*(1.-x)
>> 388 - b3*(1.-x2)/2.+b4*(1.-x2*x)/3.;
>> 389 }
>> 390 TotalCrossSection *= (twopi_mc2_rcl2*AtomicNumber/KineticEnergy);
>> 391 }
>> 392 return TotalCrossSection ;
>> 393 }
>> 394
>> 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 396
>> 397 G4VParticleChange* G4eIonisation::PostStepDoIt( const G4Track& trackData,
>> 398 const G4Step& stepData)
>> 399 {
>> 400 aParticleChange.Initialize(trackData);
>> 401
>> 402 G4Material* aMaterial = trackData.GetMaterial();
>> 403 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
>> 404
>> 405 G4double particleMass = aParticle->GetDefinition()->GetPDGMass();
>> 406 G4double Charge = aParticle->GetDefinition()->GetPDGCharge();
>> 407 G4double KineticEnergy = aParticle->GetKineticEnergy();
>> 408 G4double TotalEnergy = KineticEnergy + particleMass;
>> 409 G4double Psquare = KineticEnergy*(TotalEnergy+particleMass);
>> 410 G4double TotalMomentum = sqrt(Psquare);
>> 411 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
>> 412
>> 413 // get kinetic energy cut for the electron
>> 414 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetEnergyCuts() ;
>> 415 G4double DeltaThreshold = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
>> 416
>> 417 // some kinematics
>> 418 G4double MaxKineticEnergyTransfer;
>> 419 if (Charge < 0.) MaxKineticEnergyTransfer = 0.5*KineticEnergy;
>> 420 else MaxKineticEnergyTransfer = KineticEnergy;
>> 421
>> 422 // sampling kinetic energy of the delta ray
>> 423
>> 424 if (MaxKineticEnergyTransfer <= DeltaThreshold)
>> 425 // pathological case (should not happen, there is no change at all)
>> 426 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 427
>> 428
>> 429 // normal case
>> 430 G4double cc,y,y2,c2,b0,b1,b2,b3,b4,x,x1,grej,grejc;
>> 431
>> 432 G4double tau = KineticEnergy/particleMass;
>> 433 G4double gamma = tau+1., gamma2=gamma*gamma;
>> 434 G4double xc = DeltaThreshold/KineticEnergy, xc1=1.-xc;
>> 435
>> 436 if (Charge < 0.) // Moller (e-e-) scattering
>> 437 {
>> 438 b1=4./(9.*gamma2-10.*gamma+5.);
>> 439 b2=tau*tau*b1; b3=(2.*gamma2+2.*gamma-1.)*b1;
>> 440 cc=1.-2.*xc;
>> 441 do {
>> 442 x = xc/(1.-cc*G4UniformRand()); x1 = 1.-x;
>> 443 grej = b2*x*x-b3*x/x1+b1*gamma2/(x1*x1);
>> 444 } while (G4UniformRand()>grej);
98 } 445 }
99 AddEmModel(1, EmModel(), FluctModel()); << 446 else // Bhabha (e+e-) scattering
100 isInitialised = true; << 447 {
>> 448 y=1./(gamma+1.); y2=y*y; cc=1.-2.*y;
>> 449 b1=2.-y2; b2=cc*(3.+y2);
>> 450 c2=cc*cc; b4=c2*cc; b3=c2+b4;
>> 451 b0=gamma2/(gamma2-1.);
>> 452 grejc=(((b4*xc-b3)*xc+b2)*xc-b1)*xc+b0;
>> 453 do {
>> 454 x = xc/(1.-xc1*G4UniformRand());
>> 455 grej = ((((b4*x-b3)*x+b2)*x-b1)*x+b0)/grejc;
>> 456 } while (G4UniformRand()>grej);
>> 457 }
>> 458
>> 459 G4double DeltaKineticEnergy = x * KineticEnergy;
>> 460
>> 461 // protection :do not produce a secondary with 0. kinetic energy !
>> 462 if (DeltaKineticEnergy <= 0.)
>> 463 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 464
>> 465 G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy*(DeltaKineticEnergy +
>> 466 2.*electron_mass_c2 ));
>> 467
>> 468 G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
>> 469 /(DeltaTotalMomentum * TotalMomentum);
>> 470
>> 471 if (costheta < -1.) costheta = -1.;
>> 472 if (costheta > +1.) costheta = +1.;
>> 473
>> 474 // direction of the delta electron
>> 475
>> 476 G4double phi = twopi * G4UniformRand();
>> 477 G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
>> 478 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
>> 479
>> 480 G4ThreeVector DeltaDirection(dirx,diry,dirz);
>> 481 DeltaDirection.rotateUz(ParticleDirection);
>> 482
>> 483 // create G4DynamicParticle object for delta ray
>> 484
>> 485 G4DynamicParticle* theDeltaRay = new G4DynamicParticle;
>> 486 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
>> 487 theDeltaRay->SetMomentumDirection(
>> 488 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z());
>> 489 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 490
>> 491 // fill aParticleChange
>> 492 // changed energy and momentum of the actual particle
>> 493 G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy;
>> 494
>> 495 G4double Edep = 0.;
>> 496
>> 497 if (finalKineticEnergy > MinKineticEnergy)
>> 498 {
>> 499 G4double finalPx = TotalMomentum*ParticleDirection.x()
>> 500 - DeltaTotalMomentum*DeltaDirection.x();
>> 501 G4double finalPy = TotalMomentum*ParticleDirection.y()
>> 502 - DeltaTotalMomentum*DeltaDirection.y();
>> 503 G4double finalPz = TotalMomentum*ParticleDirection.z()
>> 504 - DeltaTotalMomentum*DeltaDirection.z();
>> 505 G4double finalMomentum =
>> 506 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
>> 507 finalPx /= finalMomentum;
>> 508 finalPy /= finalMomentum;
>> 509 finalPz /= finalMomentum;
>> 510
>> 511 aParticleChange.SetMomentumChange(finalPx, finalPy, finalPz);
>> 512 }
>> 513 else
>> 514 {
>> 515 Edep = finalKineticEnergy;
>> 516 finalKineticEnergy = 0.;
>> 517 if (Charge < 0.) aParticleChange.SetStatusChange(fStopAndKill);
>> 518 else aParticleChange.SetStatusChange(fStopButAlive);
>> 519 }
>> 520
>> 521 aParticleChange.SetEnergyChange(finalKineticEnergy);
>> 522 aParticleChange.SetNumberOfSecondaries(1);
>> 523 aParticleChange.AddSecondary(theDeltaRay);
>> 524 aParticleChange.SetLocalEnergyDeposit(Edep);
>> 525
>> 526 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 527 }
>> 528
>> 529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 530
>> 531 G4bool G4eIonisation::StorePhysicsTable(G4ParticleDefinition* particle,
>> 532 const G4String& directory,
>> 533 G4bool ascii)
>> 534 {
>> 535 G4String filename;
>> 536
>> 537 // store stopping power table
>> 538 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 539 if ( !theLossTable->StorePhysicsTable(filename, ascii) ){
>> 540 G4cout << " FAIL theLossTable->StorePhysicsTable in " << filename
>> 541 << G4endl;
>> 542 return false;
101 } 543 }
>> 544 // store mean free path table
>> 545 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 546 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 547 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 548 << G4endl;
>> 549 return false;
>> 550 }
>> 551
>> 552 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 553 << ": Success to store the PhysicsTables in "
>> 554 << directory << G4endl;
>> 555 return true;
102 } 556 }
103 557
104 //....oooOO0OOooo........oooOO0OOooo........oo << 558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 559
106 void G4eIonisation::ProcessDescription(std::os << 560 G4bool G4eIonisation::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 561 const G4String& directory,
>> 562 G4bool ascii)
107 { 563 {
108 out << " Ionisation"; << 564 // delete theLossTable and theMeanFreePathTable
109 G4VEnergyLossProcess::ProcessDescription(out << 565 if (theLossTable != 0) {
>> 566 theLossTable->clearAndDestroy();
>> 567 delete theLossTable;
>> 568 }
>> 569 if (theMeanFreePathTable != 0) {
>> 570 theMeanFreePathTable->clearAndDestroy();
>> 571 delete theMeanFreePathTable;
>> 572 }
>> 573
>> 574 // get bining from EnergyLoss
>> 575 LowestKineticEnergy = GetLowerBoundEloss();
>> 576 HighestKineticEnergy = GetUpperBoundEloss();
>> 577 TotBin = GetNbinEloss();
>> 578
>> 579 G4String filename;
>> 580
>> 581 // retreive stopping power table
>> 582 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 583 theLossTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 584 if ( !theLossTable->RetrievePhysicsTable(filename, ascii) ){
>> 585 G4cout << " FAIL theLossTable->RetrievePhysicsTable in " << filename
>> 586 << G4endl;
>> 587 return false;
>> 588 }
>> 589
>> 590 // retreive mean free path table
>> 591 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 592 theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
>> 593 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 594 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 595 << G4endl;
>> 596 return false;
>> 597 }
>> 598
>> 599 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 600 << ": Success to retrieve the PhysicsTables from "
>> 601 << directory << G4endl;
>> 602
>> 603 if (particle==G4Electron::Electron())
>> 604 {
>> 605 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
>> 606 CounterOfElectronProcess++;
>> 607 }
>> 608 else
>> 609 {
>> 610 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
>> 611 CounterOfPositronProcess++;
>> 612 }
>> 613
>> 614 BuildDEDXTable(*particle);
>> 615
>> 616 return true;
110 } 617 }
>> 618
>> 619 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 620
>> 621 void G4eIonisation::PrintInfoDefinition()
>> 622 {
>> 623 G4String comments = "delta cross sections from Moller+Bhabha. "
>> 624 "Good description from 1 KeV to 100 GeV.\n"
>> 625 " delta ray energy sampled from differential Xsection.";
>> 626
>> 627 G4cout << G4endl << GetProcessName() << ": " << comments
>> 628 << "\n PhysicsTables from "
>> 629 << G4BestUnit(LowerBoundLambda,"Energy")
>> 630 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 631 << " in " << NbinLambda << " bins. \n";
>> 632 }
111 633
112 //....oooOO0OOooo........oooOO0OOooo........oo 634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 635