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