#include <G4WentzelVIRelModel.hh>
Inheritance diagram for G4WentzelVIRelModel:
Public Member Functions | |
G4WentzelVIRelModel (const G4String &nam="WentzelVIUni") | |
virtual | ~G4WentzelVIRelModel () |
virtual void | Initialise (const G4ParticleDefinition *, const G4DataVector &) |
void | StartTracking (G4Track *) |
virtual G4double | ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX) |
virtual G4ThreeVector & | SampleScattering (const G4ThreeVector &, G4double safety) |
virtual G4double | ComputeTruePathLengthLimit (const G4Track &track, G4double ¤tMinimalStep) |
virtual G4double | ComputeGeomPathLength (G4double truePathLength) |
virtual G4double | ComputeTrueStepLength (G4double geomStepLength) |
Definition at line 69 of file G4WentzelVIRelModel.hh.
G4WentzelVIRelModel::G4WentzelVIRelModel | ( | const G4String & | nam = "WentzelVIUni" |
) |
Definition at line 71 of file G4WentzelVIRelModel.cc.
References G4Pow::GetInstance(), G4NistManager::Instance(), and G4LossTableManager::Instance().
00071 : 00072 G4VMscModel(nam), 00073 numlimit(0.1), 00074 currentCouple(0), 00075 cosThetaMin(1.0), 00076 inside(false), 00077 singleScatteringMode(false) 00078 { 00079 invsqrt12 = 1./sqrt(12.); 00080 tlimitminfix = 1.e-6*mm; 00081 lowEnergyLimit = 1.0*eV; 00082 particle = 0; 00083 nelments = 5; 00084 xsecn.resize(nelments); 00085 prob.resize(nelments); 00086 theManager = G4LossTableManager::Instance(); 00087 fNistManager = G4NistManager::Instance(); 00088 fG4pow = G4Pow::GetInstance(); 00089 wokvi = new G4WentzelVIRelXSection(); 00090 00091 preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange = xtsec = 0; 00092 currentMaterialIndex = 0; 00093 cosThetaMax = cosTetMaxNuc = 1.0; 00094 00095 fParticleChange = 0; 00096 currentCuts = 0; 00097 currentMaterial = 0; 00098 }
G4WentzelVIRelModel::~G4WentzelVIRelModel | ( | ) | [virtual] |
G4double G4WentzelVIRelModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | KineticEnergy, | |||
G4double | AtomicNumber, | |||
G4double | AtomicWeight = 0. , |
|||
G4double | cut = DBL_MAX , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 131 of file G4WentzelVIRelModel.cc.
References G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom(), G4VEmModel::CurrentCouple(), FatalException, G4Exception(), G4lrint(), G4WentzelVIRelXSection::SetupKinematic(), and G4WentzelVIRelXSection::SetupTarget().
00136 { 00137 G4double cross = 0.0; 00138 if(p != particle) { SetupParticle(p); } 00139 if(kinEnergy < lowEnergyLimit) { return cross; } 00140 if(!CurrentCouple()) { 00141 G4Exception("G4WentzelVIRelModel::ComputeCrossSectionPerAtom", "em0011", 00142 FatalException, " G4MaterialCutsCouple is not defined"); 00143 return 0.0; 00144 } 00145 DefineMaterial(CurrentCouple()); 00146 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial); 00147 if(cosTetMaxNuc < 1.0) { 00148 cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cutEnergy); 00149 cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc); 00150 /* 00151 if(p->GetParticleName() == "e-") 00152 G4cout << "G4WentzelVIRelModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy 00153 << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn 00154 << " " << particle->GetParticleName() << G4endl; 00155 */ 00156 } 00157 return cross; 00158 }
Reimplemented from G4VMscModel.
Definition at line 266 of file G4WentzelVIRelModel.cc.
References DBL_MAX, G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4WentzelVIRelXSection::SetupKinematic().
00267 { 00268 tPathLength = truelength; 00269 zPathLength = tPathLength; 00270 00271 if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) { 00272 G4double tau = tPathLength/lambdaeff; 00273 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength 00274 // << " Leff= " << lambdaeff << " tau= " << tau << G4endl; 00275 // small step 00276 if(tau < numlimit) { 00277 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0); 00278 00279 // medium step 00280 } else { 00281 G4double e1 = 0.0; 00282 if(currentRange > tPathLength) { 00283 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple); 00284 } 00285 e1 = 0.5*(e1 + preKinEnergy); 00286 cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial); 00287 lambdaeff = GetTransportMeanFreePath(particle,e1); 00288 zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff)); 00289 } 00290 } else { lambdaeff = DBL_MAX; } 00291 //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl; 00292 return zPathLength; 00293 }
G4double G4WentzelVIRelModel::ComputeTruePathLengthLimit | ( | const G4Track & | track, | |
G4double & | currentMinimalStep | |||
) | [virtual] |
Reimplemented from G4VMscModel.
Definition at line 170 of file G4WentzelVIRelModel.cc.
References G4VMscModel::ComputeGeomLimit(), G4VMscModel::ComputeSafety(), G4VMscModel::ConvertTrueToGeom(), G4VMscModel::facgeom, G4VMscModel::facrange, G4VMscModel::facsafety, fGeomBoundary, fUseDistanceToBoundary, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4Step::GetPreStepPoint(), G4ProductionCuts::GetProductionCut(), G4MaterialCutsCouple::GetProductionCuts(), G4Material::GetRadlen(), G4VMscModel::GetRange(), G4Track::GetStep(), G4VMscModel::GetTransportMeanFreePath(), G4WentzelVIRelXSection::SetupKinematic(), G4InuclParticleNames::sp, and G4VMscModel::steppingAlgorithm.
00173 { 00174 G4double tlimit = currentMinimalStep; 00175 const G4DynamicParticle* dp = track.GetDynamicParticle(); 00176 G4StepPoint* sp = track.GetStep()->GetPreStepPoint(); 00177 G4StepStatus stepStatus = sp->GetStepStatus(); 00178 singleScatteringMode = false; 00179 //G4cout << "G4WentzelVIRelModel::ComputeTruePathLengthLimit stepStatus= " 00180 // << stepStatus << G4endl; 00181 00182 00183 // initialisation for each step, lambda may be computed from scratch 00184 preKinEnergy = dp->GetKineticEnergy(); 00185 DefineMaterial(track.GetMaterialCutsCouple()); 00186 lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy); 00187 currentRange = GetRange(particle,preKinEnergy,currentCouple); 00188 cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial); 00189 00190 // extra check for abnormal situation 00191 // this check needed to run MSC with eIoni and eBrem inactivated 00192 if(tlimit > currentRange) { tlimit = currentRange; } 00193 00194 // stop here if small range particle 00195 if(inside || tlimit < tlimitminfix) { 00196 return ConvertTrueToGeom(tlimit, currentMinimalStep); 00197 } 00198 00199 // pre step 00200 G4double presafety = sp->GetSafety(); 00201 // far from geometry boundary 00202 if(currentRange < presafety) { 00203 inside = true; 00204 return ConvertTrueToGeom(tlimit, currentMinimalStep); 00205 } 00206 00207 // compute presafety again if presafety <= 0 and no boundary 00208 // i.e. when it is needed for optimization purposes 00209 if(stepStatus != fGeomBoundary && presafety < tlimitminfix) { 00210 presafety = ComputeSafety(sp->GetPosition(), tlimit); 00211 if(currentRange < presafety) { 00212 inside = true; 00213 return ConvertTrueToGeom(tlimit, currentMinimalStep); 00214 } 00215 } 00216 /* 00217 G4cout << "e(MeV)= " << preKinEnergy/MeV 00218 << " " << particle->GetParticleName() 00219 << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm 00220 << " R(mm)= " <<currentRange/mm 00221 << " L0(mm^-1)= " << lambdaeff*mm 00222 <<G4endl; 00223 */ 00224 00225 // natural limit for high energy 00226 G4double rlimit = std::max(facrange*currentRange, 00227 0.7*(1.0 - cosTetMaxNuc)*lambdaeff); 00228 00229 // low-energy e- 00230 if(cosThetaMax > cosTetMaxNuc) { 00231 rlimit = std::min(rlimit, facsafety*presafety); 00232 } 00233 00234 // cut correction 00235 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1); 00236 //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " << presafety 00237 // << " 1-cosThetaMax= " <<1-cosThetaMax << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 00238 // << G4endl; 00239 if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); } 00240 00241 if(rlimit < tlimit) { tlimit = rlimit; } 00242 00243 tlimit = std::max(tlimit, tlimitminfix); 00244 00245 // step limit in infinite media 00246 tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom); 00247 00248 //compute geomlimit and force few steps within a volume 00249 if (steppingAlgorithm == fUseDistanceToBoundary && stepStatus == fGeomBoundary) 00250 { 00251 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange); 00252 tlimit = std::min(tlimit, geomlimit/facgeom); 00253 } 00254 00255 /* 00256 G4cout << particle->GetParticleName() << " e= " << preKinEnergy 00257 << " L0= " << lambdaeff << " R= " << currentRange 00258 << "tlimit= " << tlimit 00259 << " currentMinimalStep= " << currentMinimalStep << G4endl; 00260 */ 00261 return ConvertTrueToGeom(tlimit, currentMinimalStep); 00262 }
Reimplemented from G4VMscModel.
Definition at line 297 of file G4WentzelVIRelModel.cc.
References DBL_MAX, G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4WentzelVIRelXSection::SetupKinematic().
00298 { 00299 // initialisation of single scattering x-section 00300 xtsec = 0.0; 00301 cosThetaMin = cosTetMaxNuc; 00302 00303 //G4cout << "Step= " << geomStepLength << " Lambda= " << lambdaeff 00304 // << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl; 00305 // pathalogical case 00306 if(lambdaeff == DBL_MAX) { 00307 singleScatteringMode = true; 00308 zPathLength = geomStepLength; 00309 tPathLength = geomStepLength; 00310 cosThetaMin = 1.0; 00311 00312 // normal case 00313 } else { 00314 00315 // small step use only single scattering 00316 const G4double singleScatLimit = 1.0e-7; 00317 if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) { 00318 singleScatteringMode = true; 00319 cosThetaMin = 1.0; 00320 zPathLength = geomStepLength; 00321 tPathLength = geomStepLength; 00322 00323 // step defined by transportation 00324 } else if(geomStepLength != zPathLength) { 00325 00326 // step defined by transportation 00327 zPathLength = geomStepLength; 00328 G4double tau = geomStepLength/lambdaeff; 00329 tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); 00330 00331 // energy correction for a big step 00332 if(tau > numlimit) { 00333 G4double e1 = 0.0; 00334 if(currentRange > tPathLength) { 00335 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple); 00336 } 00337 e1 = 0.5*(e1 + preKinEnergy); 00338 cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial); 00339 lambdaeff = GetTransportMeanFreePath(particle,e1); 00340 tau = zPathLength/lambdaeff; 00341 00342 if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 00343 else { tPathLength = currentRange; } 00344 } 00345 } 00346 } 00347 00348 // check of step length 00349 // define threshold angle between single and multiple scattering 00350 if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; } 00351 00352 // recompute transport cross section - do not change energy 00353 // anymore - cannot be applied for big steps 00354 if(cosThetaMin > cosTetMaxNuc) { 00355 00356 // new computation 00357 G4double cross = ComputeXSectionPerVolume(); 00358 //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec << G4endl; 00359 if(cross <= 0.0) { 00360 singleScatteringMode = true; 00361 tPathLength = zPathLength; 00362 lambdaeff = DBL_MAX; 00363 } else if(xtsec > 0.0) { 00364 00365 lambdaeff = 1./cross; 00366 G4double tau = zPathLength*cross; 00367 if(tau < numlimit) { tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); } 00368 else if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 00369 else { tPathLength = currentRange; } 00370 00371 if(tPathLength > currentRange) { tPathLength = currentRange; } 00372 } 00373 } 00374 00375 /* 00376 G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength 00377 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl; 00378 G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin 00379 << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 00380 << " e(MeV)= " << preKinEnergy/MeV << " " << singleScatteringMode << G4endl; 00381 */ 00382 return tPathLength; 00383 }
void G4WentzelVIRelModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 109 of file G4WentzelVIRelModel.cc.
References G4VMscModel::GetParticleChangeForMSC(), G4WentzelVIRelXSection::Initialise(), and G4VEmModel::PolarAngleLimit().
00111 { 00112 // reset parameters 00113 SetupParticle(p); 00114 currentRange = 0.0; 00115 00116 cosThetaMax = cos(PolarAngleLimit()); 00117 wokvi->Initialise(p, cosThetaMax); 00118 /* 00119 G4cout << "G4WentzelVIRelModel: " << particle->GetParticleName() 00120 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax 00121 << G4endl; 00122 */ 00123 currentCuts = &cuts; 00124 00125 // set values of some data members 00126 fParticleChange = GetParticleChangeForMSC(p); 00127 }
G4ThreeVector & G4WentzelVIRelModel::SampleScattering | ( | const G4ThreeVector & | , | |
G4double | safety | |||
) | [virtual] |
Reimplemented from G4VMscModel.
Definition at line 388 of file G4WentzelVIRelModel.cc.
References DBL_MAX, G4VMscModel::fDisplacement, G4lrint(), G4UniformRand, G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4VMscModel::latDisplasment, G4ParticleChangeForMSC::ProposeMomentumDirection(), G4WentzelVIRelXSection::SampleSingleScattering(), and G4WentzelVIRelXSection::SetupTarget().
00390 { 00391 fDisplacement.set(0.0,0.0,0.0); 00392 //G4cout << "!##! G4WentzelVIRelModel::SampleScattering for " 00393 // << particle->GetParticleName() << G4endl; 00394 00395 // ignore scattering for zero step length and energy below the limit 00396 if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0) 00397 { return fDisplacement; } 00398 00399 G4double invlambda = 0.0; 00400 if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; } 00401 00402 // use average kinetic energy over the step 00403 G4double cut = (*currentCuts)[currentMaterialIndex]; 00404 /* 00405 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV 00406 << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec 00407 << " xmsc= " << tPathLength*invlambda << " safety= " << safety << G4endl; 00408 */ 00409 00410 // step limit due msc 00411 G4double x0 = tPathLength; 00412 // const G4double thinlimit = 0.5; 00413 const G4double thinlimit = 0.1; 00414 G4int nMscSteps = 1; 00415 // large scattering angle case - two step approach 00416 if(tPathLength*invlambda > thinlimit && safety > tlimitminfix) { 00417 x0 *= 0.5; 00418 nMscSteps = 2; 00419 } 00420 00421 // step limit due to single scattering 00422 G4double x1 = 2*tPathLength; 00423 if(0.0 < xtsec) { x1 = -log(G4UniformRand())/xtsec; } 00424 00425 const G4ElementVector* theElementVector = 00426 currentMaterial->GetElementVector(); 00427 G4int nelm = currentMaterial->GetNumberOfElements(); 00428 00429 // geometry 00430 G4double sint, cost, phi; 00431 G4ThreeVector temp(0.0,0.0,1.0); 00432 00433 // current position and direction relative to the end point 00434 // because of magnetic field geometry is computed relatively to the 00435 // end point of the step 00436 G4ThreeVector dir(0.0,0.0,1.0); 00437 fDisplacement.set(0.0,0.0,-zPathLength); 00438 G4double mscfac = zPathLength/tPathLength; 00439 00440 // start a loop 00441 G4double x2 = x0; 00442 G4double step, z; 00443 G4bool singleScat; 00444 /* 00445 G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2 00446 << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode 00447 << " xtsec= " << xtsec << G4endl; 00448 */ 00449 do { 00450 00451 // single scattering case 00452 if(x1 < x2) { 00453 step = x1; 00454 singleScat = true; 00455 } else { 00456 step = x2; 00457 singleScat = false; 00458 } 00459 00460 // new position 00461 fDisplacement += step*mscfac*dir; 00462 00463 if(singleScat) { 00464 00465 // select element 00466 G4int i = 0; 00467 if(nelm > 1) { 00468 G4double qsec = G4UniformRand()*xtsec; 00469 for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } } 00470 } 00471 G4double cosTetM = 00472 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut); 00473 //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " " << prob[i] << G4endl; 00474 temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]); 00475 00476 // direction is changed 00477 temp.rotateUz(dir); 00478 dir = temp; 00479 00480 // new proposed step length 00481 x1 = -log(G4UniformRand())/xtsec; 00482 x2 -= step; 00483 if(x2 <= 0.0) { --nMscSteps; } 00484 00485 // multiple scattering 00486 } else { 00487 --nMscSteps; 00488 x1 -= step; 00489 x2 = x0; 00490 00491 // for multiple scattering x0 should be used as a step size 00492 if(!singleScatteringMode) { 00493 G4double z0 = x0*invlambda; 00494 00495 // correction to keep first moment 00496 00497 // sample z in interval 0 - 1 00498 if(z0 > 5.0) { z = G4UniformRand(); } 00499 else { 00500 G4double zzz = 0.0; 00501 if(z0 > 0.01) { zzz = exp(-1.0/z0); } 00502 z = -z0*log(1.0 - (1.0 - zzz)*G4UniformRand()); 00503 // /(1.0 - (1.0/z0 + 1.0)*zzz); 00504 } 00505 00506 cost = 1.0 - 2.0*z/*factCM*/; 00507 if(cost > 1.0) { cost = 1.0; } 00508 else if(cost < -1.0) { cost =-1.0; } 00509 sint = sqrt((1.0 - cost)*(1.0 + cost)); 00510 phi = twopi*G4UniformRand(); 00511 G4double vx1 = sint*cos(phi); 00512 G4double vy1 = sint*sin(phi); 00513 00514 // change direction 00515 temp.set(vx1,vy1,cost); 00516 temp.rotateUz(dir); 00517 dir = temp; 00518 00519 // lateral displacement 00520 if (latDisplasment && x0 > tlimitminfix) { 00521 G4double rms = invsqrt12*sqrt(2*z0); 00522 G4double r = x0*mscfac; 00523 G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0)); 00524 G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0)); 00525 G4double dz; 00526 G4double d = r*r - dx*dx - dy*dy; 00527 if(d >= 0.0) { dz = sqrt(d) - r; } 00528 else { dx = dy = dz = 0.0; } 00529 00530 // change position 00531 temp.set(dx,dy,dz); 00532 temp.rotateUz(dir); 00533 fDisplacement += temp; 00534 } 00535 } 00536 } 00537 } while (0 < nMscSteps); 00538 00539 dir.rotateUz(oldDirection); 00540 00541 //G4cout << "G4WentzelVIRelModel sampling of scattering is done" << G4endl; 00542 // end of sampling ------------------------------- 00543 00544 fParticleChange->ProposeMomentumDirection(dir); 00545 00546 // lateral displacement 00547 fDisplacement.rotateUz(oldDirection); 00548 00549 /* 00550 G4cout << " r(mm)= " << fDisplacement.mag() 00551 << " safety= " << safety 00552 << " trueStep(mm)= " << tPathLength 00553 << " geomStep(mm)= " << zPathLength 00554 << " x= " << fDisplacement.x() 00555 << " y= " << fDisplacement.y() 00556 << " z= " << fDisplacement.z() 00557 << G4endl; 00558 */ 00559 00560 //G4cout<< "G4WentzelVIRelModel::SampleScattering end NewDir= " << dir<< G4endl; 00561 return fDisplacement; 00562 }
void G4WentzelVIRelModel::StartTracking | ( | G4Track * | ) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 162 of file G4WentzelVIRelModel.cc.
References G4DynamicParticle::GetDefinition(), and G4Track::GetDynamicParticle().
00163 { 00164 SetupParticle(track->GetDynamicParticle()->GetDefinition()); 00165 inside = false; 00166 }