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