00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 #include "G4WentzelVIModel.hh"
00059 #include "G4PhysicalConstants.hh"
00060 #include "G4SystemOfUnits.hh"
00061 #include "Randomize.hh"
00062 #include "G4ParticleChangeForMSC.hh"
00063 #include "G4PhysicsTableHelper.hh"
00064 #include "G4ElementVector.hh"
00065 #include "G4ProductionCutsTable.hh"
00066 #include "G4LossTableManager.hh"
00067 #include "G4Pow.hh"
00068 #include "G4NistManager.hh"
00069
00070
00071
00072 using namespace std;
00073
00074 G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
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 }
00102
00103
00104
00105 G4WentzelVIModel::~G4WentzelVIModel()
00106 {
00107 delete wokvi;
00108 }
00109
00110
00111
00112 void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
00113 const G4DataVector& cuts)
00114 {
00115
00116 SetupParticle(p);
00117 currentRange = 0.0;
00118
00119 cosThetaMax = cos(PolarAngleLimit());
00120 wokvi->Initialise(p, cosThetaMax);
00121
00122
00123
00124
00125
00126 currentCuts = &cuts;
00127
00128
00129 fParticleChange = GetParticleChangeForMSC(p);
00130 }
00131
00132
00133
00134 G4double G4WentzelVIModel::ComputeCrossSectionPerAtom(
00135 const G4ParticleDefinition* p,
00136 G4double kinEnergy,
00137 G4double Z, G4double,
00138 G4double cutEnergy, G4double)
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
00155
00156
00157
00158
00159 }
00160 return cross;
00161 }
00162
00163
00164
00165 void G4WentzelVIModel::StartTracking(G4Track* track)
00166 {
00167 SetupParticle(track->GetDynamicParticle()->GetDefinition());
00168 inside = false;
00169 }
00170
00171
00172
00173 G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
00174 const G4Track& track,
00175 G4double& currentMinimalStep)
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
00183
00184
00185
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
00193
00194
00195
00196
00197 if(tlimit > currentRange) { tlimit = currentRange; }
00198
00199
00200 if(inside || tlimit < tlimitminfix) {
00201 return ConvertTrueToGeom(tlimit, currentMinimalStep);
00202 }
00203
00204
00205 G4double presafety = sp->GetSafety();
00206
00207 if(currentRange < presafety) {
00208 inside = true;
00209 return ConvertTrueToGeom(tlimit, currentMinimalStep);
00210 }
00211
00212
00213
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
00223
00224
00225
00226
00227
00228
00229
00230
00231 G4double rlimit = std::max(facrange*currentRange,
00232 0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
00233
00234
00235 if(cosThetaMax > cosTetMaxNuc) {
00236 rlimit = std::min(rlimit, facsafety*presafety);
00237 }
00238
00239
00240 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
00241
00242
00243
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
00251 tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
00252
00253
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
00263
00264
00265
00266
00267 return ConvertTrueToGeom(tlimit, currentMinimalStep);
00268 }
00269
00270
00271
00272 G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
00273 {
00274 tPathLength = truelength;
00275 zPathLength = tPathLength;
00276
00277 if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
00278 G4double tau = tPathLength/lambdaeff;
00279
00280
00281
00282 if(tau < numlimit) {
00283 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
00284
00285
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
00298 return zPathLength;
00299 }
00300
00301
00302
00303 G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
00304 {
00305
00306 xtsec = 0.0;
00307 cosThetaMin = cosTetMaxNuc;
00308
00309
00310
00311
00312
00313
00314 if(lambdaeff == DBL_MAX) {
00315 singleScatteringMode = true;
00316 zPathLength = geomStepLength;
00317 tPathLength = geomStepLength;
00318 cosThetaMin = 1.0;
00319
00320
00321 } else {
00322
00323
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
00332 } else if(geomStepLength != zPathLength) {
00333
00334
00335 zPathLength = geomStepLength;
00336 G4double tau = geomStepLength/lambdaeff;
00337 tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
00338
00339
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
00357
00358 if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; }
00359
00360
00361
00362 if(cosThetaMin > cosTetMaxNuc) {
00363
00364
00365 G4double cross = ComputeXSectionPerVolume();
00366
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
00386
00387
00388
00389
00390
00391
00392 return tPathLength;
00393 }
00394
00395
00396
00397 G4ThreeVector&
00398 G4WentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection,
00399 G4double safety)
00400 {
00401 fDisplacement.set(0.0,0.0,0.0);
00402
00403
00404
00405
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
00413 G4double cut = (*currentCuts)[currentMaterialIndex];
00414
00415
00416
00417
00418
00419
00420
00421
00422 G4double x0 = tPathLength;
00423
00424 const G4double thinlimit = 0.1;
00425 G4int nMscSteps = 1;
00426
00427 if(!singleScatteringMode && tPathLength*invlambda > thinlimit
00428 && safety > tlimitminfix) {
00429 x0 *= 0.5;
00430 nMscSteps = 2;
00431 }
00432
00433
00434 G4double x1 = 2*tPathLength;
00435 if(0.0 < xtsec) { x1 = -log(G4UniformRand())/xtsec; }
00436
00437
00438 if(singleScatteringMode && x1 > tPathLength)
00439 { return fDisplacement; }
00440
00441 const G4ElementVector* theElementVector =
00442 currentMaterial->GetElementVector();
00443 G4int nelm = currentMaterial->GetNumberOfElements();
00444
00445
00446 G4double sint, cost, phi;
00447 G4ThreeVector temp(0.0,0.0,1.0);
00448
00449
00450
00451
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
00458 G4double x2 = x0;
00459 G4double step, z;
00460 G4bool singleScat;
00461
00462
00463
00464
00465
00466 do {
00467
00468
00469 if(singleScatteringMode || x1 <= x2) {
00470 step = x1;
00471 singleScat = true;
00472 } else {
00473 step = x2;
00474 singleScat = false;
00475 }
00476
00477
00478 fDisplacement += step*mscfac*dir;
00479
00480 if(singleScat) {
00481
00482
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
00491 temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
00492
00493
00494 temp.rotateUz(dir);
00495 dir = temp;
00496
00497
00498 x2 -= step;
00499 x1 = -log(G4UniformRand())/xtsec;
00500 if(singleScatteringMode && x1 > x2) { break; }
00501
00502
00503 } else {
00504 --nMscSteps;
00505 x1 -= step;
00506 x2 = x0;
00507
00508
00509 G4double z0 = x0*invlambda;
00510
00511
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;
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
00526 temp.set(vx1,vy1,cost);
00527 temp.rotateUz(dir);
00528 dir = temp;
00529
00530
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
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
00552
00553
00554 fParticleChange->ProposeMomentumDirection(dir);
00555
00556
00557 fDisplacement.rotateUz(oldDirection);
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 return fDisplacement;
00572 }
00573
00574
00575
00576 G4double G4WentzelVIModel::ComputeXSectionPerVolume()
00577 {
00578
00579 const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
00580 const G4double* theAtomNumDensityVector =
00581 currentMaterial->GetVecNbOfAtomsPerVolume();
00582 G4int nelm = currentMaterial->GetNumberOfElements();
00583 if(nelm > nelments) {
00584 nelments = nelm;
00585 xsecn.resize(nelm);
00586 prob.resize(nelm);
00587 }
00588 G4double cut = (*currentCuts)[currentMaterialIndex];
00589
00590
00591
00592 xtsec = 0.0;
00593 if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
00594
00595
00596 G4double xs = 0.0;
00597 for (G4int i=0; i<nelm; ++i) {
00598 G4double costm =
00599 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
00600 G4double density = theAtomNumDensityVector[i];
00601
00602 G4double esec = 0.0;
00603 if(costm < cosThetaMin) {
00604
00605
00606 if(1.0 > cosThetaMin) {
00607 xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
00608 }
00609
00610 G4double nucsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
00611 esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
00612 nucsec += esec;
00613 if(nucsec > 0.0) { esec /= nucsec; }
00614 xtsec += nucsec*density;
00615 }
00616 xsecn[i] = xtsec;
00617 prob[i] = esec;
00618
00619
00620
00621 }
00622
00623
00624
00625 return xs;
00626 }
00627
00628