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