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