| 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: G4MuIonisation.cc,v 1.11 2000/05/23 09:58:47 urban Exp $
7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-03-00 $
8 // * LICENSE and available at http://cern.ch/ << 10 //
9 // * include a list of copyright holders. << 11 //
10 // * << 12 // --------------------------------------------------------------
11 // * Neither the authors of this software syst << 13 // GEANT 4 class implementation file
12 // * institutes,nor the agencies providing fin << 14 //
13 // * work make any representation or warran << 15 // For information related to this code contact:
14 // * regarding this software system or assum << 16 // CERN, CN Division, ASD group
15 // * use. Please see the license in the file << 17 // History: first implementation, based on object model of
16 // * for the full disclaimer and the limitatio << 18 // 2nd December 1995, G.Cosmo
17 // * << 19 // ------------ G4MuIonisation physics process -------------
18 // * This code implementation is the result << 20 // by Laszlo Urban, September 1997
19 // * technical work of the GEANT4 collaboratio << 21 // ------------------------------------------------------------------
20 // * By using, copying, modifying or distri << 22 // It is the implementation of the NEW IONISATION PROCESS.
21 // * any work based on the software) you ag << 23 // It calculates the ionisation of muons.
22 // * use in resulting scientific publicati << 24 // **************************************************************
23 // * acceptance of all terms of the Geant4 Sof << 25 // 08-04-98: remove 'tracking cut' of the ionizing particle, MMa
24 // ******************************************* << 26 // 26/10/98: new stuff from R.Kokoulin + cleanup , L.Urban
25 // << 27 // 10/02/00 modifications , new e.m. structure, L.Urban
26 // ------------------------------------------- << 28 // --------------------------------------------------------------
27 // << 29
28 // GEANT4 Class file <<
29 // <<
30 // <<
31 // File name: G4MuIonisation <<
32 // <<
33 // Author: Laszlo Urban <<
34 // <<
35 // Creation date: 30.09.1997 <<
36 // <<
37 // Modifications: <<
38 // <<
39 // 08-04-98 remove 'tracking cut' of the ioniz <<
40 // 26-10-98 new stuff from R.Kokoulin + cleanu <<
41 // 10-02-00 modifications , new e.m. structure <<
42 // 23-03-01 R.Kokoulin's correction is comment <<
43 // 29-05-01 V.Ivanchenko minor changes to prov <<
44 // 10-08-01 new methods Store/Retrieve Physics <<
45 // 28-08-01 new function ComputeRestrictedMean <<
46 // 17-09-01 migration of Materials to pure STL <<
47 // 26-09-01 completion of RetrievePhysicsTable <<
48 // 29-10-01 all static functions no more inlin <<
49 // 07-11-01 correction(Tmax+xsection computati <<
50 // 08-11-01 particleMass becomes a local varia <<
51 // 10-05-02 V.Ivanchenko update to new design <<
52 // 04-12-02 V.Ivanchenko the low energy limit <<
53 // 23-12-02 Change interface in order to move <<
54 // 26-12-02 Secondary production moved to deri <<
55 // 13-02-03 SubCutoff regime is assigned to a <<
56 // 23-05-03 Define default integral + BohrFluc <<
57 // 03-06-03 Add SetIntegral method to choose f <<
58 // 03-06-03 Fix initialisation problem for STD <<
59 // 04-08-03 Set integral=false to be default ( <<
60 // 08-08-03 STD substitute standard (V.Ivanch <<
61 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossPro <<
62 // 10-02-04 Calculation of radiative correctio <<
63 // 27-05-04 Set integral to be a default regim <<
64 // 17-08-04 Utilise mu+ tables for mu- (V.Ivan <<
65 // 08-11-04 Migration to new interface of Stor <<
66 // 08-04-05 Major optimisation of internal int <<
67 // 12-08-05 SetStepLimits(0.2, 0.1*mm) (mma) <<
68 // 02-09-05 SetStepLimits(0.2, 1*mm) (V.Ivantc <<
69 // 12-08-05 SetStepLimits(0.2, 0.1*mm) + integ <<
70 // 10-01-06 SetStepLimits -> SetStepFunction ( <<
71 // <<
72 // ------------------------------------------- <<
73 // <<
74 //....oooOO0OOooo........oooOO0OOooo........oo <<
75 //....oooOO0OOooo........oooOO0OOooo........oo <<
76 30
77 #include "G4MuIonisation.hh" 31 #include "G4MuIonisation.hh"
78 #include "G4PhysicalConstants.hh" << 32 #include "G4UnitsTable.hh"
79 #include "G4SystemOfUnits.hh" <<
80 #include "G4Electron.hh" <<
81 #include "G4BraggModel.hh" <<
82 #include "G4BetheBlochModel.hh" <<
83 #include "G4MuBetheBlochModel.hh" <<
84 #include "G4EmStandUtil.hh" <<
85 #include "G4ICRU73QOModel.hh" <<
86 #include "G4EmParameters.hh" <<
87 33
88 //....oooOO0OOooo........oooOO0OOooo........oo << 34 #include "G4ios.hh"
89 35
90 G4MuIonisation::G4MuIonisation(const G4String& <<
91 : G4VEnergyLossProcess(name) <<
92 { <<
93 SetProcessSubType(fIonisation); <<
94 SetSecondaryParticle(G4Electron::Electron()) <<
95 } <<
96 36
97 //....oooOO0OOooo........oooOO0OOooo........oo << 37 G4double G4MuIonisation::LowerBoundLambda = 1.*keV ;
>> 38 G4double G4MuIonisation::UpperBoundLambda = 1000000.*TeV ;
>> 39 G4int G4MuIonisation::NbinLambda = 150 ;
98 40
99 G4bool G4MuIonisation::IsApplicable(const G4Pa <<
100 { <<
101 return (p.GetPDGCharge() != 0.0); <<
102 } <<
103 41
104 //....oooOO0OOooo........oooOO0OOooo........oo << 42 // constructor and destructor
>> 43
>> 44 G4MuIonisation::G4MuIonisation(const G4String& processName)
>> 45 : G4VMuEnergyLoss(processName),
>> 46 theMeanFreePathTable(NULL)
>> 47 { }
105 48
106 G4double G4MuIonisation::MinPrimaryEnergy(cons << 49 G4MuIonisation::~G4MuIonisation()
107 const G4Material*, <<
108 G4double cut) <<
109 { 50 {
110 G4double x = 0.5*cut/CLHEP::electron_mass_c2 << 51 if (theMeanFreePathTable) {
111 G4double gam = x*ratio + std::sqrt((1. + x)* << 52 theMeanFreePathTable->clearAndDestroy();
112 return mass*(gam - 1.0); << 53 delete theMeanFreePathTable;
113 } << 54 }
114 55
115 //....oooOO0OOooo........oooOO0OOooo........oo << 56 }
116 57
117 void << 58 void G4MuIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
118 G4MuIonisation::InitialiseEnergyLossProcess(co << 59 // just call BuildLossTable+BuildLambdaTable
119 co <<
120 { 60 {
121 if(!isInitialised) { << 61 // get bining from EnergyLoss
>> 62 LowestKineticEnergy = GetLowerBoundEloss() ;
>> 63 HighestKineticEnergy = GetUpperBoundEloss() ;
>> 64 TotBin = GetNbinEloss() ;
>> 65
>> 66 BuildLossTable(aParticleType) ;
>> 67 G4double Charge = aParticleType.GetPDGCharge();
>> 68
>> 69 if(Charge>0.)
>> 70 {
>> 71 RecorderOfmuplusProcess[CounterOfmuplusProcess] = (*this).theLossTable ;
>> 72 CounterOfmuplusProcess++;
>> 73 }
>> 74 else
>> 75 {
>> 76 RecorderOfmuminusProcess[CounterOfmuminusProcess] = (*this).theLossTable ;
>> 77 CounterOfmuminusProcess++;
>> 78 }
>> 79
>> 80 G4double electronCutInRange = G4Electron::Electron()->GetCuts();
>> 81 if(electronCutInRange != lastelectronCutInRange)
>> 82 BuildLambdaTable(aParticleType) ;
>> 83
>> 84 G4VMuEnergyLoss::BuildDEDXTable(aParticleType) ;
122 85
123 theParticle = part; << 86 if(&aParticleType == theMuonPlus)
124 theBaseParticle = bpart; << 87 PrintInfoDefinition() ;
>> 88 }
125 89
126 mass = theParticle->GetPDGMass(); << 90 void G4MuIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType)
127 ratio = CLHEP::electron_mass_c2/mass; << 91 {
128 G4double q = theParticle->GetPDGCharge(); << 92 DeltaCutInKineticEnergy = theElectron->GetCutsInEnergy() ;
129 93
130 G4EmParameters* param = G4EmParameters::In << 94 G4double LowEdgeEnergy , ionloss ;
131 G4double elow = 0.2*CLHEP::MeV; << 95 G4double deltaloss ;
132 G4double emax = param->MaxKinEnergy(); << 96 G4double RateMass ;
>> 97 G4bool isOutRange ;
>> 98 static const G4MaterialTable* theMaterialTable=
>> 99 G4Material::GetMaterialTable();
>> 100 const G4double twoln10 = 2.*log(10.) ;
>> 101 const G4double Factor = twopi_mc2_rcl2 ;
>> 102 const G4double bg2lim = 0.0169 , taulim = 8.4146e-3 ;
>> 103
>> 104 ParticleMass = aParticleType.GetPDGMass() ;
>> 105 RateMass = electron_mass_c2/ParticleMass ;
>> 106
>> 107 G4int numOfMaterials = theMaterialTable->length();
>> 108
>> 109 if ( theLossTable) {
>> 110 theLossTable->clearAndDestroy();
>> 111 delete theLossTable;
>> 112 }
>> 113 theLossTable = new G4PhysicsTable(numOfMaterials);
133 114
134 // Bragg peak model << 115 for (G4int J=0; J<numOfMaterials; J++)
135 if (nullptr == EmModel(0)) { << 116 {
136 if(q > 0.0) { SetEmModel(new G4BraggMode << 117 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
137 else { SetEmModel(new G4ICRU73QOM << 118 LowestKineticEnergy, HighestKineticEnergy, TotBin);
>> 119
>> 120 G4double ElectronDensity,Eexc,Eexc2,Cden,Mden,Aden,X0den,X1den,taul ;
>> 121 G4double* ShellCorrectionVector;
>> 122
>> 123 const G4Material* material= (*theMaterialTable)[J];
>> 124
>> 125 ElectronDensity = material->GetElectronDensity();
>> 126 Eexc = material->GetIonisation()->GetMeanExcitationEnergy();
>> 127 Eexc2 = Eexc*Eexc ;
>> 128 Cden = material->GetIonisation()->GetCdensity();
>> 129 Mden = material->GetIonisation()->GetMdensity();
>> 130 Aden = material->GetIonisation()->GetAdensity();
>> 131 X0den = material->GetIonisation()->GetX0density();
>> 132 X1den = material->GetIonisation()->GetX1density();
>> 133 taul = material->GetIonisation()->GetTaul() ;
>> 134 ShellCorrectionVector = material->GetIonisation()
>> 135 ->GetShellCorrectionVector();
>> 136
>> 137 const G4ElementVector* theElementVector=
>> 138 material->GetElementVector() ;
>> 139 const G4double* theAtomicNumDensityVector=
>> 140 material->GetAtomicNumDensityVector() ;
>> 141 const G4int NumberOfElements=
>> 142 material->GetNumberOfElements() ;
>> 143 DeltaCutInKineticEnergyNow = DeltaCutInKineticEnergy[J] ;
>> 144
>> 145 G4double tau,tau0,Tmax,gamma,bg2,beta2,rcut,delta,x,sh ;
>> 146 for (G4int i = 0 ; i < TotBin ; i++)
>> 147 {
>> 148 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 149 tau = LowEdgeEnergy/ParticleMass ;
>> 150 gamma = tau +1. ;
>> 151 bg2 = tau*(tau+2.) ;
>> 152 beta2 = bg2/(gamma*gamma) ;
>> 153 Tmax = 2.*electron_mass_c2*bg2
>> 154 /(1.+2.*gamma*RateMass+RateMass*RateMass) ;
>> 155
>> 156 if ( tau < taul )
>> 157 // low energy part , parametrized energy loss formulae
>> 158 {
>> 159 ionloss = 0. ;
>> 160 deltaloss = 0. ;
>> 161
>> 162 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 163 {
>> 164 const G4Element* element = (*theElementVector)(iel);
>> 165
>> 166 if ( tau < element->GetIonisation()->GetTau0())
>> 167 ionloss += theAtomicNumDensityVector[iel]
>> 168 *( element->GetIonisation()->GetAlow()*sqrt(tau)
>> 169 +element->GetIonisation()->GetBlow()*tau) ;
>> 170 else
>> 171 ionloss += theAtomicNumDensityVector[iel]
>> 172 * element->GetIonisation()->GetClow()/sqrt(tau) ;
>> 173 }
>> 174 if ( DeltaCutInKineticEnergyNow < Tmax)
>> 175 {
>> 176 deltaloss = log(Tmax/DeltaCutInKineticEnergyNow)-
>> 177 beta2*(1.-DeltaCutInKineticEnergyNow/Tmax) ;
>> 178 if(aParticleType.GetPDGSpin() == 0.5)
>> 179 deltaloss += 0.25*(Tmax-DeltaCutInKineticEnergyNow)*
>> 180 (Tmax-DeltaCutInKineticEnergyNow)/
>> 181 (LowEdgeEnergy*LowEdgeEnergy+proton_mass_c2*proton_mass_c2) ;
>> 182 deltaloss *= Factor*ElectronDensity/beta2 ;
>> 183 }
>> 184 ionloss -= deltaloss ;
>> 185 }
>> 186 else
>> 187 // high energy part , Bethe-Bloch formula
>> 188 {
>> 189
>> 190 if ( DeltaCutInKineticEnergyNow < Tmax)
>> 191 rcut = DeltaCutInKineticEnergyNow/Tmax ;
>> 192 else
>> 193 rcut = 1.;
>> 194
>> 195 ionloss = log(2.*electron_mass_c2*bg2*Tmax/Eexc2)
>> 196 +log(rcut)-(1.+rcut)*beta2 ;
>> 197
>> 198 // density correction
>> 199 x = log(bg2)/twoln10 ;
>> 200 if ( x < X0den )
>> 201 delta = 0. ;
>> 202 else
>> 203 {
>> 204 delta = twoln10*x - Cden ;
>> 205 if ( x < X1den )
>> 206 delta += Aden*pow((X1den-x),Mden) ;
>> 207 }
>> 208
>> 209 // shell correction
>> 210 if ( bg2 > bg2lim ) {
>> 211 sh = 0. ;
>> 212 x = 1. ;
>> 213 for (G4int k=0; k<=2; k++) {
>> 214 x *= bg2 ;
>> 215 sh += ShellCorrectionVector[k]/x;
>> 216 }
>> 217 }
>> 218 else {
>> 219 sh = 0. ;
>> 220 x = 1. ;
>> 221 for (G4int k=0; k<=2; k++) {
>> 222 x *= bg2lim ;
>> 223 sh += ShellCorrectionVector[k]/x;
>> 224 }
>> 225 sh *= log(tau/taul)/log(taulim/taul) ;
>> 226 }
>> 227
>> 228 ionloss -= delta + sh ;
>> 229 ionloss /= beta2 ;
>> 230
>> 231 // correction of R. Kokoulin
>> 232 G4double E = LowEdgeEnergy+ParticleMass ;
>> 233 G4double epmax = RateMass*E*E/(RateMass*E+ParticleMass) ;
>> 234 G4double apar = log(2.*epmax/electron_mass_c2) ;
>> 235 ionloss += fine_structure_const*(log(2.*E/ParticleMass)-apar/3.)*
>> 236 apar*apar/twopi ;
>> 237
>> 238 ionloss *= Factor*ElectronDensity ;
>> 239 }
>> 240 if ( ionloss <= 0.)
>> 241 ionloss = 0. ;
>> 242 aVector->PutValue(i,ionloss) ;
138 } 243 }
139 EmModel(0)->SetLowEnergyLimit(param->MinKi << 244 theLossTable->insert(aVector);
140 EmModel(0)->SetHighEnergyLimit(elow); << 245 }
>> 246 }
141 247
142 // fluctuation model << 248 void G4MuIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType)
143 if (nullptr == FluctModel()) { << 249 {
144 SetFluctModel(G4EmStandUtil::ModelOfFluc << 250 // Build mean free path tables for the delta ray production process
>> 251 G4double LowEdgeEnergy,Tmax , Value ,sigma ;
>> 252 G4bool isOutRange ;
>> 253 const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable();
>> 254
>> 255 G4int numOfMaterials = theMaterialTable->length();
>> 256
>> 257 if (theMeanFreePathTable) {
>> 258 theMeanFreePathTable->clearAndDestroy();
>> 259 delete theMeanFreePathTable;
>> 260 }
>> 261
>> 262 theMeanFreePathTable = new G4PhysicsTable(numOfMaterials);
>> 263
>> 264 // get electron and particle cuts in kinetic energy
>> 265 DeltaCutInKineticEnergy = theElectron->GetCutsInEnergy() ;
>> 266
>> 267 for (G4int J=0 ; J < numOfMaterials; J++)
>> 268 {
>> 269 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 270 LowerBoundLambda,UpperBoundLambda,NbinLambda);
>> 271
>> 272 const G4Material* material= (*theMaterialTable)[J];
>> 273 const G4ElementVector* theElementVector=
>> 274 material->GetElementVector() ;
>> 275 const G4double* theAtomicNumDensityVector =
>> 276 material->GetAtomicNumDensityVector();
>> 277 const G4int NumberOfElements=
>> 278 material->GetNumberOfElements() ;
>> 279 DeltaCutInKineticEnergyNow = DeltaCutInKineticEnergy[J] ;
>> 280
>> 281 for ( G4int i = 0 ; i < NbinLambda ; i++ )
>> 282 {
>> 283 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 284 sigma = 0. ;
>> 285
>> 286 // check threshold here !
>> 287 G4double Tmax = 2.*electron_mass_c2*LowEdgeEnergy*
>> 288 (LowEdgeEnergy+2.*ParticleMass)/
>> 289 (ParticleMass*ParticleMass+2.*electron_mass_c2*
>> 290 (LowEdgeEnergy+ParticleMass)+
>> 291 electron_mass_c2*electron_mass_c2) ;
>> 292
>> 293 if(Tmax > DeltaCutInKineticEnergyNow)
>> 294 {
>> 295 for (G4int iel=0; iel<NumberOfElements; iel++ )
>> 296 {
>> 297 sigma += theAtomicNumDensityVector[iel]*
>> 298 ComputeMicroscopicCrossSection(aParticleType,
>> 299 LowEdgeEnergy,
>> 300 (*theElementVector)(iel)->GetZ() ) ;
>> 301 }
>> 302 }
>> 303
>> 304 Value = sigma<=0 ? DBL_MAX : 1./sigma ;
>> 305
>> 306 aVector->PutValue(i, Value) ;
145 } 307 }
146 AddEmModel(1, EmModel(0), FluctModel()); <<
147 308
148 // high energy model << 309 theMeanFreePathTable->insert(aVector);
149 if (nullptr == EmModel(1)) { SetEmModel(ne << 310 }
150 EmModel(1)->SetLowEnergyLimit(elow); << 311 }
151 EmModel(1)->SetHighEnergyLimit(emax); << 312
152 AddEmModel(1, EmModel(1), FluctModel()); <<
153 313
154 isInitialised = true; << 314 G4double G4MuIonisation::ComputeMicroscopicCrossSection(
>> 315 const G4ParticleDefinition& aParticleType,
>> 316 G4double KineticEnergy,
>> 317 G4double AtomicNumber)
>> 318 {
>> 319 const G4double xgi[] = {0.06943,0.33001,0.66999,0.93057} ;
>> 320 const G4double wgi[] = {0.17393,0.32607,0.32607,0.17393} ;
>> 321 const G4double ak1 = 4.6 ;
>> 322 const G4int k2 = 2 ;
>> 323 const G4double masspar = 0.5*ParticleMass*ParticleMass/electron_mass_c2 ;
>> 324 G4double TotalEnergy=KineticEnergy + ParticleMass;
>> 325 G4double KnockonMaxEnergy = TotalEnergy/(1.+masspar/TotalEnergy) ;
>> 326
>> 327 G4double TotalCrossSection= 0. ;
>> 328
>> 329 if( KnockonMaxEnergy > DeltaCutInKineticEnergyNow )
>> 330 {
>> 331 G4double aaa = log(DeltaCutInKineticEnergyNow);
>> 332 G4double bbb = log(KnockonMaxEnergy) ;
>> 333 G4int kkk = int((bbb-aaa)/ak1)+k2 ;
>> 334 G4double hhh = (bbb-aaa)/kkk ;
>> 335 G4double step = exp(hhh) ;
>> 336 G4double ymax = 1./KnockonMaxEnergy ;
>> 337
>> 338 for (G4int k=0; k<kkk; k++)
>> 339 {
>> 340 G4double ymin = ymax ;
>> 341 ymax = ymin*step ;
>> 342 G4double hhy = ymax-ymin ;
>> 343 for (G4int i=0; i<4; i++)
>> 344 {
>> 345 G4double y = ymin+hhy*xgi[i];
>> 346 G4double ep = 1./y ;
>> 347 TotalCrossSection += ep*ep*wgi[i]*hhy*
>> 348 ComputeDMicroscopicCrossSection(
>> 349 aParticleType,KineticEnergy,
>> 350 AtomicNumber,ep) ;
>> 351 }
>> 352 }
155 } 353 }
>> 354 return TotalCrossSection ;
156 } 355 }
>> 356
>> 357 G4double G4MuIonisation::ComputeDMicroscopicCrossSection(
>> 358 const G4ParticleDefinition& ParticleType,
>> 359 G4double KineticEnergy, G4double AtomicNumber,
>> 360 G4double KnockonEnergy)
>> 361 // Calculates the differential (D) microscopic cross section
>> 362 // using the cross section formula of R.P. Kokoulin (10/98)
>> 363 {
>> 364 const G4double masspar=0.5*ParticleMass*ParticleMass/electron_mass_c2 ;
>> 365 const G4double alphaprime = fine_structure_const/twopi ;
>> 366 G4double TotalEnergy = KineticEnergy + ParticleMass ;
>> 367 G4double KnockonMaxEnergy = TotalEnergy/(1.+masspar/TotalEnergy) ;
>> 368
>> 369 G4double DCrossSection = 0. ;
>> 370
>> 371 if(KnockonEnergy >= KnockonMaxEnergy) return DCrossSection ;
>> 372
>> 373 G4double v = KnockonEnergy/TotalEnergy ;
>> 374 DCrossSection = twopi_mc2_rcl2*AtomicNumber*
>> 375 (1.-KnockonEnergy/KnockonMaxEnergy+0.5*v*v)/
>> 376 (KnockonEnergy*KnockonEnergy) ;
>> 377 G4double a1 = log(1.+2.*KnockonEnergy/electron_mass_c2) ;
>> 378 G4double a3 = log(4.*TotalEnergy*(TotalEnergy-KnockonEnergy)/
>> 379 (ParticleMass*ParticleMass)) ;
>> 380 DCrossSection *= (1.+alphaprime*a1*(a3-a1)) ;
157 381
158 //....oooOO0OOooo........oooOO0OOooo........oo << 382 return DCrossSection ;
>> 383 }
>> 384
>> 385 G4VParticleChange* G4MuIonisation::PostStepDoIt(
>> 386 const G4Track& trackData,
>> 387 const G4Step& stepData)
>> 388 {
>> 389 const G4DynamicParticle* aParticle ;
>> 390 const G4double alphaprime = fine_structure_const/twopi ;
>> 391 G4Material* aMaterial;
>> 392 G4double KineticEnergy,TotalEnergy,TotalMomentum,
>> 393 betasquare,MaxKineticEnergyTransfer,
>> 394 DeltaKineticEnergy,DeltaTotalMomentum,costheta,sintheta,phi,
>> 395 dirx,diry,dirz,finalKineticEnergy,finalPx,finalPy,finalPz,
>> 396 x,xc,te2,grej,Psquare,Esquare,summass,rate,grejc,finalMomentum ;
>> 397 G4double Charge ;
>> 398
>> 399 aParticleChange.Initialize(trackData) ;
>> 400 aMaterial = trackData.GetMaterial() ;
>> 401 aParticle = trackData.GetDynamicParticle() ;
>> 402 Charge=aParticle->GetDefinition()->GetPDGCharge();
>> 403 KineticEnergy=aParticle->GetKineticEnergy();
>> 404 TotalEnergy=KineticEnergy + ParticleMass ;
>> 405 Psquare=KineticEnergy*(TotalEnergy+ParticleMass) ;
>> 406 Esquare=TotalEnergy*TotalEnergy ;
>> 407 summass = ParticleMass + electron_mass_c2 ;
>> 408 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection() ;
>> 409
>> 410 DeltaCutInKineticEnergyNow = DeltaCutInKineticEnergy[aMaterial->GetIndex()];
>> 411
>> 412 // some kinematics......................
>> 413 betasquare=Psquare/Esquare ;
>> 414 MaxKineticEnergyTransfer = 2.*electron_mass_c2*Psquare
>> 415 /(summass*summass+2.*electron_mass_c2*KineticEnergy);
>> 416
>> 417 // sampling kinetic energy of the delta ray
>> 418 if( MaxKineticEnergyTransfer <= DeltaCutInKineticEnergyNow )
>> 419 {
>> 420 // pathological case (it should not happen ,
>> 421 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 422 }
>> 423 else
>> 424 {
>> 425 // normal case ......................................
>> 426 xc=DeltaCutInKineticEnergyNow/MaxKineticEnergyTransfer ;
>> 427 rate=MaxKineticEnergyTransfer/TotalEnergy ;
>> 428 te2=0.5*rate*rate ;
>> 429
>> 430 // sampling follows ...
>> 431 G4double a0=log(2.*TotalEnergy/ParticleMass) ;
>> 432 grejc=(1.-betasquare*xc+te2*xc*xc)*
>> 433 (1.+ alphaprime*a0*a0) ;
>> 434 do {
>> 435 x=xc/(1.-(1.-xc)*G4UniformRand());
>> 436 G4double twoep = 2.*x*MaxKineticEnergyTransfer ;
>> 437 grej=(1.-x*(betasquare-x*te2))*
>> 438 (1.+alphaprime*log(1.+twoep/electron_mass_c2)*
>> 439 (a0+log((2.*TotalEnergy-twoep)/ParticleMass)-
>> 440 log(1.+twoep/electron_mass_c2)))
>> 441 /grejc ;
>> 442
>> 443 } while( G4UniformRand()>grej );
>> 444 }
>> 445
>> 446 DeltaKineticEnergy = x * MaxKineticEnergyTransfer ;
>> 447
>> 448 if(DeltaKineticEnergy <= 0.)
>> 449 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 450
>> 451 DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
>> 452 2. * electron_mass_c2 )) ;
>> 453 TotalMomentum = sqrt(Psquare) ;
>> 454 costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
>> 455 /(DeltaTotalMomentum * TotalMomentum) ;
>> 456
>> 457 // protection against costheta > 1 or < -1 ---------------
>> 458 if ( costheta < -1. )
>> 459 costheta = -1. ;
>> 460 if ( costheta > +1. )
>> 461 costheta = +1. ;
>> 462
>> 463 // direction of the delta electron ........
>> 464 phi = twopi * G4UniformRand() ;
>> 465 sintheta = sqrt((1.+costheta)*(1.-costheta));
>> 466 dirx = sintheta * cos(phi) ;
>> 467 diry = sintheta * sin(phi) ;
>> 468 dirz = costheta ;
>> 469 G4ThreeVector DeltaDirection(dirx,diry,dirz) ;
>> 470 DeltaDirection.rotateUz(ParticleDirection) ;
>> 471
>> 472 // create G4DynamicParticle object for delta ray
>> 473 G4DynamicParticle *theDeltaRay = new G4DynamicParticle;
>> 474 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
>> 475 theDeltaRay->SetMomentumDirection(
>> 476 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z());
>> 477 theDeltaRay->SetDefinition(G4Electron::Electron());
>> 478 finalKineticEnergy = KineticEnergy - DeltaKineticEnergy ;
>> 479 if (finalKineticEnergy > 0. )
>> 480 {
>> 481 finalPx = TotalMomentum*ParticleDirection.x()
>> 482 - DeltaTotalMomentum*DeltaDirection.x();
>> 483 finalPy = TotalMomentum*ParticleDirection.y()
>> 484 - DeltaTotalMomentum*DeltaDirection.y();
>> 485 finalPz = TotalMomentum*ParticleDirection.z()
>> 486 - DeltaTotalMomentum*DeltaDirection.z();
>> 487 finalMomentum =
>> 488 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
>> 489 finalPx /= finalMomentum ;
>> 490 finalPy /= finalMomentum ;
>> 491 finalPz /= finalMomentum ;
>> 492
>> 493 aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz );
>> 494 }
>> 495 else
>> 496 {
>> 497 finalKineticEnergy = 0. ;
>> 498 aParticleChange.SetStatusChange(fStopButAlive);
>> 499 }
>> 500
>> 501 aParticleChange.SetEnergyChange( finalKineticEnergy );
>> 502 aParticleChange.SetNumberOfSecondaries(1);
>> 503 aParticleChange.AddSecondary( theDeltaRay );
>> 504 aParticleChange.SetLocalEnergyDeposit (0.);
>> 505
>> 506 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 507 }
159 508
160 void G4MuIonisation::ProcessDescription(std::o << 509 void G4MuIonisation::PrintInfoDefinition()
161 { 510 {
162 out << " Muon ionisation"; << 511 G4String comments = "knock-on electron cross sections .\n ";
163 G4VEnergyLossProcess::ProcessDescription(out << 512 comments += " Good description above the mean excitation energy.\n";
>> 513 comments += " delta ray energy sampled from differential Xsection." ;
>> 514
>> 515 G4cout << G4endl << GetProcessName() << ": " << comments
>> 516 << "\n PhysicsTables from " << G4BestUnit(LowerBoundLambda,
>> 517 "Energy")
>> 518 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 519 << " in " << NbinLambda << " bins. \n";
164 } 520 }
165 521
166 //....oooOO0OOooo........oooOO0OOooo........oo <<
167 522