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 << 6 // * the Geant4 Collaboration. It is provided << 7 // * conditions of the Geant4 Software License << 8 // * LICENSE and available at http://cern.ch/ << 9 // * include a list of copyright holders. << 10 // * << 11 // * Neither the authors of this software syst << 12 // * institutes,nor the agencies providing fin << 13 // * work make any representation or warran << 14 // * regarding this software system or assum << 15 // * use. Please see the license in the file << 16 // * for the full disclaimer and the limitatio << 17 // * << 18 // * This code implementation is the result << 19 // * technical work of the GEANT4 collaboratio << 20 // * By using, copying, modifying or distri << 21 // * any work based on the software) you ag << 22 // * use in resulting scientific publicati << 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* << 25 // 7 // 26 // ------------------------------------------- << 8 // $Id: G4eIonisation.cc,v 1.10 2000/05/23 14:42:21 urban Exp $ >> 9 // GEANT4 tag $Name: geant4-03-00 $ 27 // 10 // 28 // GEANT4 Class file << 11 // >> 12 // ------------------------------------------------------------- >> 13 // GEANT 4 class implementation file 29 // 14 // >> 15 // For information related to this code contact: >> 16 // CERN, IT Division, ASD group >> 17 // History: based on object model of >> 18 // 2nd December 1995, G.Cosmo >> 19 // ---------- G4eIonisation physics process ----------- >> 20 // by Laszlo Urban, 20 March 1997 >> 21 // ************************************************************** >> 22 // It is the first implementation of the NEW IONISATION PROCESS. >> 23 // It calculates the ionisation of e+/e-. >> 24 // ************************************************************** 30 // 25 // 31 // File name: G4eIonisation << 26 // 07-04-98: remove 'tracking cut' of the ionizing particle, MMa 32 // << 27 // 04-09-98: new methods SetBining() PrintInfo() 33 // Author: Laszlo Urban << 28 // 07-09-98: Cleanup 34 // << 29 // 02/02/99: correction inDoIt , L.Urban 35 // Creation date: 20.03.1997 << 30 // 10/02/00 modifications , new e.m. structure, L.Urban 36 // << 31 // -------------------------------------------------------------- 37 // Modified by Michel Maire and Vladimir Ivanc << 32 38 // << 39 // ------------------------------------------- << 40 // << 41 //....oooOO0OOooo........oooOO0OOooo........oo << 42 //....oooOO0OOooo........oooOO0OOooo........oo << 43 << 44 #include "G4eIonisation.hh" 33 #include "G4eIonisation.hh" 45 #include "G4Electron.hh" << 34 #include "G4EnergyLossTables.hh" 46 #include "G4MollerBhabhaModel.hh" << 35 #include "G4ios.hh" 47 #include "G4EmParameters.hh" << 36 #include "G4UnitsTable.hh" 48 #include "G4EmStandUtil.hh" << 37 >> 38 G4double G4eIonisation::LowerBoundLambda = 1.*keV ; >> 39 G4double G4eIonisation::UpperBoundLambda = 100.*TeV ; >> 40 G4int G4eIonisation::NbinLambda = 100 ; >> 41 >> 42 // constructor and destructor >> 43 >> 44 G4eIonisation::G4eIonisation(const G4String& processName) >> 45 : G4VeEnergyLoss(processName), >> 46 theMeanFreePathTable(NULL) >> 47 { } 49 48 50 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 50 52 G4eIonisation::G4eIonisation(const G4String& n << 51 G4eIonisation::~G4eIonisation() 53 : G4VEnergyLossProcess(name), << 54 theElectron(G4Electron::Electron()), << 55 isElectron(true), << 56 isInitialised(false) << 57 { 52 { 58 SetProcessSubType(fIonisation); << 53 if (theMeanFreePathTable) { 59 SetSecondaryParticle(theElectron); << 54 theMeanFreePathTable->clearAndDestroy(); >> 55 delete theMeanFreePathTable; >> 56 } >> 57 60 } 58 } 61 59 62 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 61 64 G4eIonisation::~G4eIonisation() = default; << 62 void G4eIonisation::BuildPhysicsTable(const G4ParticleDefinition& aParticleType) >> 63 // just call BuildLossTable+BuildLambdaTable >> 64 { >> 65 // get bining from EnergyLoss >> 66 LowestKineticEnergy = GetLowerBoundEloss() ; >> 67 HighestKineticEnergy = GetUpperBoundEloss() ; >> 68 TotBin = GetNbinEloss() ; >> 69 >> 70 BuildLossTable(aParticleType) ; >> 71 >> 72 if(&aParticleType==G4Electron::Electron()) >> 73 { >> 74 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable ; >> 75 CounterOfElectronProcess++; >> 76 } >> 77 else >> 78 { >> 79 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable ; >> 80 CounterOfPositronProcess++; >> 81 } >> 82 >> 83 BuildLambdaTable(aParticleType) ; >> 84 >> 85 BuildDEDXTable(aParticleType) ; >> 86 >> 87 if(&aParticleType==G4Electron::Electron()) >> 88 PrintInfoDefinition(); >> 89 } 65 90 66 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 67 92 68 G4double G4eIonisation::MinPrimaryEnergy(const << 93 void G4eIonisation::BuildLossTable(const G4ParticleDefinition& aParticleType) 69 const G4Material*, << 70 G4double cut) << 71 { 94 { 72 G4double x = cut; << 95 // Build tables for the ionization energy loss 73 if(isElectron) x += cut; << 96 // the tables are built for *MATERIALS* 74 return x; << 97 >> 98 const G4double twoln10 = 2.*log(10.); >> 99 const G4double Factor = twopi_mc2_rcl2; >> 100 >> 101 static const G4double Tl = 0.2*keV ; >> 102 >> 103 G4double LowEdgeEnergy, ionloss; >> 104 >> 105 // material properties >> 106 G4double ElectronDensity,Eexc,Eexcm2,Cden,Mden,Aden,X0den,X1den ; >> 107 // some local variables >> 108 G4double tau,Tmax,gamma,gamma2,bg2,beta2,d,d2,d3,d4,delta,x,y ; >> 109 >> 110 ParticleMass = aParticleType.GetPDGMass(); >> 111 G4double* ParticleCutInKineticEnergy = aParticleType.GetEnergyCuts() ; >> 112 >> 113 // create table >> 114 >> 115 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); >> 116 G4int numOfMaterials = theMaterialTable->length(); >> 117 >> 118 if (theLossTable) { theLossTable->clearAndDestroy(); >> 119 delete theLossTable; >> 120 } >> 121 >> 122 theLossTable = new G4PhysicsTable(numOfMaterials); >> 123 >> 124 // loop for materials >> 125 >> 126 for (G4int J=0; J<numOfMaterials; J++) >> 127 { >> 128 // create physics vector and fill it >> 129 >> 130 G4PhysicsLogVector* aVector = new G4PhysicsLogVector( >> 131 LowestKineticEnergy, HighestKineticEnergy, TotBin); >> 132 >> 133 // get material parameters needed for the energy loss calculation >> 134 const G4Material* material= (*theMaterialTable)[J]; >> 135 >> 136 ElectronDensity = material->GetElectronDensity(); >> 137 Eexc = material->GetIonisation()->GetMeanExcitationEnergy(); >> 138 Eexc /= ParticleMass; Eexcm2 = Eexc*Eexc; >> 139 Cden = material->GetIonisation()->GetCdensity(); >> 140 Mden = material->GetIonisation()->GetMdensity(); >> 141 Aden = material->GetIonisation()->GetAdensity(); >> 142 X0den = material->GetIonisation()->GetX0density(); >> 143 X1den = material->GetIonisation()->GetX1density(); >> 144 >> 145 // for the lowenergy extrapolation >> 146 G4double Zeff = material->GetTotNbOfElectPerVolume()/ >> 147 material->GetTotNbOfAtomsPerVolume() ; >> 148 G4double Th = 0.25*sqrt(Zeff)*keV ; >> 149 >> 150 G4double Tsav ; >> 151 >> 152 // now comes the loop for the kinetic energy values >> 153 >> 154 for (G4int i = 0 ; i < TotBin ; i++) >> 155 { >> 156 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ; >> 157 // low energy ? >> 158 if(LowEdgeEnergy < Th) >> 159 { >> 160 Tsav = LowEdgeEnergy ; >> 161 LowEdgeEnergy = Th ; >> 162 } >> 163 else >> 164 Tsav = 0. ; >> 165 >> 166 tau = LowEdgeEnergy/ParticleMass ; >> 167 >> 168 // Seltzer-Berger formula >> 169 gamma = tau + 1.; gamma2 = gamma*gamma; >> 170 bg2 = tau*(tau+2.); >> 171 beta2 = bg2/gamma2; >> 172 >> 173 // electron >> 174 if (&aParticleType==G4Electron::Electron()) >> 175 { >> 176 Tmax = LowEdgeEnergy/2.; >> 177 d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass; >> 178 ionloss = log(2.*(tau+2.)/Eexcm2)-1.-beta2 >> 179 + log((tau-d)*d)+tau/(tau-d) >> 180 + (0.5*d*d+(2.*tau+1.)*log(1.-d/tau))/gamma2; >> 181 } >> 182 else //positron >> 183 { >> 184 Tmax = LowEdgeEnergy ; >> 185 d = G4std::min(ParticleCutInKineticEnergy[J], Tmax)/ParticleMass; >> 186 d2=d*d/2.; d3=d*d*d/3.; d4=d*d*d*d/4.; >> 187 y=1./(1.+gamma); >> 188 ionloss = log(2.*(tau+2.)/Eexcm2)+log(tau*d) >> 189 - beta2*(tau+2.*d-y*(3.*d2+y*(d-d3+y*(d2-tau*d3+d4))))/tau; >> 190 } >> 191 >> 192 //density correction >> 193 x = log(bg2)/twoln10; >> 194 if (x < X0den) delta = 0.; >> 195 else { delta = twoln10*x - Cden; >> 196 if (x < X1den) delta += Aden*pow((X1den-x),Mden); >> 197 } >> 198 >> 199 //now you can compute the total ionization loss >> 200 ionloss -= delta ; >> 201 ionloss *= Factor*ElectronDensity/beta2 ; >> 202 if (ionloss <= 0.) ionloss = 0.; >> 203 >> 204 // low energy ? >> 205 if(Tsav > 0.) >> 206 { >> 207 if(Tsav >= Tl) >> 208 ionloss *= sqrt(LowEdgeEnergy/Tsav) ; >> 209 else >> 210 ionloss *= sqrt(LowEdgeEnergy*Tsav)/Tl ; >> 211 } >> 212 >> 213 aVector->PutValue(i,ionloss) ; >> 214 } >> 215 theLossTable->insert(aVector); >> 216 } 75 } 217 } 76 218 77 //....oooOO0OOooo........oooOO0OOooo........oo 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 220 79 G4bool G4eIonisation::IsApplicable(const G4Par << 221 void G4eIonisation::BuildLambdaTable(const G4ParticleDefinition& aParticleType) 80 { 222 { 81 return (&p == theElectron || &p == G4Positro << 223 // Build mean free path tables for the delta ray production process >> 224 // tables are built for MATERIALS >> 225 >> 226 G4double LowEdgeEnergy, Value, SIGMA; >> 227 >> 228 //create table >> 229 const G4MaterialTable* theMaterialTable=G4Material::GetMaterialTable(); >> 230 G4int numOfMaterials = theMaterialTable->length(); >> 231 >> 232 if (theMeanFreePathTable) { theMeanFreePathTable->clearAndDestroy(); >> 233 delete theMeanFreePathTable; >> 234 } >> 235 >> 236 theMeanFreePathTable = new G4PhysicsTable(numOfMaterials); >> 237 >> 238 // get electron cuts in kinetic energy >> 239 // The electron cuts needed in the case of the positron , too! >> 240 // This is the reason why SetCut has to be called for electron first !!!!!!! >> 241 >> 242 if((G4Electron::Electron()->GetCutsInEnergy() == 0) && >> 243 ( &aParticleType == G4Positron::Positron())) >> 244 { >> 245 G4cout << " The ELECTRON energy cuts needed to compute energy loss/mean free path " >> 246 " for POSITRON , too. " << G4endl; >> 247 G4Exception(" Call SetCut for e- first !!!!!!") ; >> 248 } >> 249 >> 250 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ; >> 251 >> 252 // loop for materials >> 253 >> 254 for (G4int J=0 ; J < numOfMaterials; J++) >> 255 { >> 256 //create physics vector then fill it .... >> 257 >> 258 G4PhysicsLogVector* aVector = new G4PhysicsLogVector( >> 259 LowerBoundLambda, UpperBoundLambda, NbinLambda); >> 260 >> 261 // compute the (macroscopic) cross section first >> 262 >> 263 const G4Material* material= (*theMaterialTable)[J]; >> 264 const >> 265 G4ElementVector* theElementVector = material->GetElementVector(); >> 266 const >> 267 G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector(); >> 268 const >> 269 G4int NumberOfElements = material->GetNumberOfElements() ; >> 270 >> 271 // get the electron kinetic energy cut for the actual material, >> 272 // it will be used in ComputeMicroscopicCrossSection >> 273 // (--> it will be the same for all the elements in this material ) >> 274 G4double DeltaThreshold = DeltaCutInKineticEnergy[J] ; >> 275 >> 276 for (G4int i = 0 ; i < NbinLambda ; i++) >> 277 { >> 278 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ; >> 279 SIGMA = 0.; >> 280 for (G4int iel=0; iel<NumberOfElements; iel++ ) >> 281 { >> 282 SIGMA += theAtomicNumDensityVector[iel]* >> 283 ComputeMicroscopicCrossSection( aParticleType, >> 284 LowEdgeEnergy, >> 285 (*theElementVector)(iel)->GetZ(), >> 286 DeltaThreshold); >> 287 } >> 288 >> 289 // mean free path = 1./macroscopic cross section >> 290 Value = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX; >> 291 aVector->PutValue(i, Value) ; >> 292 } >> 293 theMeanFreePathTable->insert(aVector); >> 294 } 82 } 295 } 83 296 84 //....oooOO0OOooo........oooOO0OOooo........oo 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 298 86 void G4eIonisation::InitialiseEnergyLossProces << 299 G4double G4eIonisation::ComputeMicroscopicCrossSection( 87 const G4ParticleDefinition* part, << 300 const G4ParticleDefinition& aParticleType, 88 const G4ParticleDefinition*) << 301 G4double KineticEnergy, >> 302 G4double AtomicNumber , >> 303 G4double DeltaThreshold) >> 304 { >> 305 // calculates the microscopic cross section >> 306 //(it is called for elements , AtomicNumber = Z ) >> 307 >> 308 G4double MaxKineticEnergyTransfer, TotalCrossSection(0.); >> 309 >> 310 ParticleMass = aParticleType.GetPDGMass(); >> 311 G4double TotalEnergy = KineticEnergy + ParticleMass; >> 312 >> 313 G4double betasquare = KineticEnergy*(TotalEnergy+ParticleMass) >> 314 /(TotalEnergy*TotalEnergy); >> 315 G4double gamma = TotalEnergy/ParticleMass, gamma2 = gamma*gamma; >> 316 G4double x=DeltaThreshold/KineticEnergy, x2 = x*x; >> 317 >> 318 if (&aParticleType==G4Electron::Electron()) >> 319 MaxKineticEnergyTransfer = 0.5*KineticEnergy; >> 320 else MaxKineticEnergyTransfer = KineticEnergy; >> 321 >> 322 // now you can calculate the total cross section >> 323 >> 324 if (MaxKineticEnergyTransfer > DeltaThreshold) >> 325 { >> 326 if (&aParticleType==G4Electron::Electron()) //Moller (e-e-) scattering >> 327 { >> 328 TotalCrossSection = (gamma-1.)*(gamma-1.)*(0.5-x)/gamma2 + 1./x >> 329 - 1./(1.-x)-(2.*gamma-1.)*log((1.-x)/x)/gamma2; >> 330 TotalCrossSection /= betasquare; >> 331 } >> 332 else //Bhabha (e+e-) scattering >> 333 { >> 334 G4double y=1./(1.+gamma), y2 =y*y, y12=1.-2.*y; >> 335 G4double b1=2.-y2, b2=y12*(3.+y2), b4=y12*y12*y12, b3=b4+y12*y12; >> 336 TotalCrossSection = (1./x-1.)/betasquare+b1*log(x)+b2*(1.-x) >> 337 - b3*(1.-x2)/2.+b4*(1.-x2*x)/3.; >> 338 } >> 339 TotalCrossSection *= (twopi_mc2_rcl2*AtomicNumber/KineticEnergy); >> 340 } >> 341 return TotalCrossSection ; >> 342 } >> 343 >> 344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 345 >> 346 G4VParticleChange* G4eIonisation::PostStepDoIt( const G4Track& trackData, >> 347 const G4Step& stepData) 89 { 348 { 90 if(!isInitialised) { << 349 aParticleChange.Initialize(trackData) ; 91 if(part != theElectron) { isElectron = fal << 350 92 if (nullptr == EmModel(0)) { SetEmModel(ne << 351 G4Material* aMaterial = trackData.GetMaterial() ; 93 G4EmParameters* param = G4EmParameters::In << 352 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle() ; 94 EmModel(0)->SetLowEnergyLimit(param->MinKi << 353 95 EmModel(0)->SetHighEnergyLimit(param->MaxK << 354 ParticleMass = aParticle->GetDefinition()->GetPDGMass(); 96 if (nullptr == FluctModel()) { << 355 G4double KineticEnergy = aParticle->GetKineticEnergy(); 97 SetFluctModel(G4EmStandUtil::ModelOfFluc << 356 G4double TotalEnergy = KineticEnergy + ParticleMass; >> 357 G4double Psquare = KineticEnergy*(TotalEnergy+ParticleMass); >> 358 G4double TotalMomentum = sqrt(Psquare); >> 359 G4double Esquare=TotalEnergy*TotalEnergy; >> 360 G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection(); >> 361 >> 362 // get kinetic energy cut for the electron >> 363 G4double* DeltaCutInKineticEnergy = G4Electron::Electron()->GetCutsInEnergy() ; >> 364 G4double DeltaThreshold = DeltaCutInKineticEnergy[aMaterial->GetIndex()]; >> 365 >> 366 // some kinematics >> 367 G4double MaxKineticEnergyTransfer; >> 368 if (Charge < 0.) MaxKineticEnergyTransfer = 0.5*KineticEnergy; >> 369 else MaxKineticEnergyTransfer = KineticEnergy; >> 370 >> 371 // sampling kinetic energy of the delta ray >> 372 >> 373 if (MaxKineticEnergyTransfer <= DeltaThreshold) // pathological case (should not happen, >> 374 // there is no change at all) >> 375 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); >> 376 >> 377 >> 378 // normal case >> 379 G4double cc,y,y2,c2,b0,b1,b2,b3,b4,x,x1,grej,grejc; >> 380 >> 381 G4double tau = KineticEnergy/ParticleMass; >> 382 G4double gamma = tau+1., gamma2=gamma*gamma; >> 383 G4double xc = DeltaThreshold/KineticEnergy, xc1=1.-xc; >> 384 >> 385 if (Charge < 0.) // Moller (e-e-) scattering >> 386 { >> 387 b1=4./(9.*gamma2-10.*gamma+5.); >> 388 b2=tau*tau*b1; b3=(2.*gamma2+2.*gamma-1.)*b1; >> 389 cc=1.-2.*xc; >> 390 do { >> 391 x = xc/(1.-cc*G4UniformRand()); x1 = 1.-x; >> 392 grej = b2*x*x-b3*x/x1+b1*gamma2/(x1*x1); >> 393 } while (G4UniformRand()>grej) ; 98 } 394 } 99 AddEmModel(1, EmModel(), FluctModel()); << 395 else // Bhabha (e+e-) scattering 100 isInitialised = true; << 396 { 101 } << 397 y=1./(gamma+1.); y2=y*y; cc=1.-2.*y; >> 398 b1=2.-y2; b2=cc*(3.+y2); >> 399 c2=cc*cc; b4=c2*cc; b3=c2+b4; >> 400 b0=gamma2/(gamma2-1.); >> 401 grejc=(((b4*xc-b3)*xc+b2)*xc-b1)*xc+b0; >> 402 do { >> 403 x = xc/(1.-xc1*G4UniformRand()); >> 404 grej = ((((b4*x-b3)*x+b2)*x-b1)*x+b0)/grejc; >> 405 } while (G4UniformRand()>grej); >> 406 } >> 407 >> 408 G4double DeltaKineticEnergy = x * KineticEnergy; >> 409 >> 410 // protection :do not produce a secondary with 0. kinetic energy ! >> 411 if (DeltaKineticEnergy <= 0.) >> 412 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); >> 413 >> 414 G4double DeltaTotalMomentum = sqrt(DeltaKineticEnergy * (DeltaKineticEnergy + >> 415 2. * electron_mass_c2 )); >> 416 >> 417 G4double costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2) >> 418 /(DeltaTotalMomentum * TotalMomentum); >> 419 >> 420 if (costheta < -1.) costheta = -1.; >> 421 if (costheta > +1.) costheta = +1.; >> 422 >> 423 // direction of the delta electron >> 424 >> 425 G4double phi = twopi * G4UniformRand(); >> 426 G4double sintheta = sqrt((1.+costheta)*(1.-costheta)); >> 427 G4double dirx = sintheta * cos(phi), diry = sintheta * sin(phi), dirz = costheta; >> 428 >> 429 G4ThreeVector DeltaDirection(dirx,diry,dirz); >> 430 DeltaDirection.rotateUz(ParticleDirection); >> 431 >> 432 // create G4DynamicParticle object for delta ray >> 433 >> 434 G4DynamicParticle* theDeltaRay = new G4DynamicParticle; >> 435 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy ); >> 436 theDeltaRay->SetMomentumDirection( >> 437 DeltaDirection.x(),DeltaDirection.y(),DeltaDirection.z()); >> 438 theDeltaRay->SetDefinition(G4Electron::Electron()); >> 439 >> 440 // fill aParticleChange >> 441 // changed energy and momentum of the actual particle >> 442 G4double finalKineticEnergy = KineticEnergy - DeltaKineticEnergy; >> 443 >> 444 G4double Edep = 0. ; >> 445 >> 446 if (finalKineticEnergy > MinKineticEnergy) >> 447 { >> 448 G4double finalPx = TotalMomentum*ParticleDirection.x() >> 449 - DeltaTotalMomentum*DeltaDirection.x(); >> 450 G4double finalPy = TotalMomentum*ParticleDirection.y() >> 451 - DeltaTotalMomentum*DeltaDirection.y(); >> 452 G4double finalPz = TotalMomentum*ParticleDirection.z() >> 453 - DeltaTotalMomentum*DeltaDirection.z(); >> 454 G4double finalMomentum = >> 455 sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ; >> 456 finalPx /= finalMomentum ; >> 457 finalPy /= finalMomentum ; >> 458 finalPz /= finalMomentum ; >> 459 >> 460 aParticleChange.SetMomentumChange( finalPx,finalPy,finalPz ); >> 461 } >> 462 else >> 463 { >> 464 finalKineticEnergy = 0.; >> 465 Edep = finalKineticEnergy ; >> 466 if (Charge < 0.) aParticleChange.SetStatusChange(fStopAndKill); >> 467 else aParticleChange.SetStatusChange(fStopButAlive); >> 468 } >> 469 >> 470 aParticleChange.SetEnergyChange( finalKineticEnergy ); >> 471 aParticleChange.SetNumberOfSecondaries(1); >> 472 aParticleChange.AddSecondary( theDeltaRay ); >> 473 aParticleChange.SetLocalEnergyDeposit (Edep); >> 474 >> 475 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); 102 } 476 } 103 477 104 //....oooOO0OOooo........oooOO0OOooo........oo 478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 105 479 106 void G4eIonisation::ProcessDescription(std::os << 480 void G4eIonisation::PrintInfoDefinition() 107 { 481 { 108 out << " Ionisation"; << 482 G4String comments = "delta cross sections from Moller+Bhabha. "; 109 G4VEnergyLossProcess::ProcessDescription(out << 483 comments += "Good description from 1 KeV to 100 GeV.\n"; 110 } << 484 comments += " delta ray energy sampled from differential Xsection."; >> 485 >> 486 G4cout << G4endl << GetProcessName() << ": " << comments >> 487 << "\n PhysicsTables from " << G4BestUnit(LowerBoundLambda,"Energy") >> 488 << " to " << G4BestUnit(UpperBoundLambda,"Energy") >> 489 << " in " << NbinLambda << " bins. \n"; >> 490 } 111 491 112 //....oooOO0OOooo........oooOO0OOooo........oo << 492 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 113 493