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
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 #include "G4MuPairProductionModel.hh"
00073 #include "G4PhysicalConstants.hh"
00074 #include "G4SystemOfUnits.hh"
00075 #include "G4Electron.hh"
00076 #include "G4Positron.hh"
00077 #include "G4MuonMinus.hh"
00078 #include "G4MuonPlus.hh"
00079 #include "Randomize.hh"
00080 #include "G4Material.hh"
00081 #include "G4Element.hh"
00082 #include "G4ElementVector.hh"
00083 #include "G4ProductionCutsTable.hh"
00084 #include "G4ParticleChangeForLoss.hh"
00085 #include "G4ParticleChangeForGamma.hh"
00086
00087
00088
00089
00090
00091 G4double G4MuPairProductionModel::zdat[]={1., 4., 13., 29., 92.};
00092 G4double G4MuPairProductionModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
00093 G4double G4MuPairProductionModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
00094 1.e9, 1.e10};
00095 G4double G4MuPairProductionModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
00096 0.5917, 0.7628, 0.8983, 0.9801 };
00097 G4double G4MuPairProductionModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
00098 0.1813, 0.1569, 0.1112, 0.0506 };
00099
00100
00101
00102 using namespace std;
00103
00104 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
00105 const G4String& nam)
00106 : G4VEmModel(nam),
00107 particle(0),
00108 factorForCross(4.*fine_structure_const*fine_structure_const
00109 *classic_electr_radius*classic_electr_radius/(3.*pi)),
00110 sqrte(sqrt(exp(1.))),
00111 currentZ(0),
00112 fParticleChange(0),
00113 minPairEnergy(4.*electron_mass_c2),
00114 lowestKinEnergy(GeV),
00115 nzdat(5),
00116 ntdat(8),
00117 nbiny(1000),
00118 nmaxElements(0),
00119 ymin(-5.),
00120 ymax(0.),
00121 dy((ymax-ymin)/nbiny),
00122 samplingTablesAreFilled(false)
00123 {
00124 SetLowEnergyLimit(minPairEnergy);
00125 nist = G4NistManager::Instance();
00126
00127 theElectron = G4Electron::Electron();
00128 thePositron = G4Positron::Positron();
00129
00130 particleMass = lnZ = z13 = z23 = 0;
00131
00132 for(size_t i=0; i<1001; ++i) { ya[i] = 0.0; }
00133
00134 if(p) { SetParticle(p); }
00135 }
00136
00137
00138
00139 G4MuPairProductionModel::~G4MuPairProductionModel()
00140 {}
00141
00142
00143
00144 G4double G4MuPairProductionModel::MinEnergyCut(const G4ParticleDefinition*,
00145 const G4MaterialCutsCouple* )
00146 {
00147 return minPairEnergy;
00148 }
00149
00150
00151
00152 G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
00153 G4double kineticEnergy)
00154 {
00155 G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
00156 return maxPairEnergy;
00157 }
00158
00159
00160
00161 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
00162 const G4DataVector&)
00163 {
00164 if (!samplingTablesAreFilled) {
00165 if(p) { SetParticle(p); }
00166 MakeSamplingTables();
00167 }
00168 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00169 }
00170
00171
00172
00173 G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
00174 const G4Material* material,
00175 const G4ParticleDefinition*,
00176 G4double kineticEnergy,
00177 G4double cutEnergy)
00178 {
00179 G4double dedx = 0.0;
00180 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
00181 { return dedx; }
00182
00183 const G4ElementVector* theElementVector = material->GetElementVector();
00184 const G4double* theAtomicNumDensityVector =
00185 material->GetAtomicNumDensityVector();
00186
00187
00188 for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
00189 G4double Z = (*theElementVector)[i]->GetZ();
00190 SetCurrentElement(Z);
00191 G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
00192 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
00193 dedx += loss*theAtomicNumDensityVector[i];
00194 }
00195 if (dedx < 0.) { dedx = 0.; }
00196 return dedx;
00197 }
00198
00199
00200
00201 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
00202 G4double tkin,
00203 G4double cutEnergy,
00204 G4double tmax)
00205 {
00206 SetCurrentElement(Z);
00207 G4double loss = 0.0;
00208
00209 G4double cut = std::min(cutEnergy,tmax);
00210 if(cut <= minPairEnergy) { return loss; }
00211
00212
00213
00214 G4double ak1=6.9;
00215 G4double ak2=1.0;
00216 G4double aaa = log(minPairEnergy);
00217 G4double bbb = log(cut);
00218 G4int kkk = (G4int)((bbb-aaa)/ak1+ak2);
00219 if (kkk > 8) kkk = 8;
00220 else if (kkk < 1) { kkk = 1; }
00221 G4double hhh = (bbb-aaa)/(G4double)kkk;
00222 G4double x = aaa;
00223
00224 for (G4int l=0 ; l<kkk; l++)
00225 {
00226
00227 for (G4int ll=0; ll<8; ll++)
00228 {
00229 G4double ep = exp(x+xgi[ll]*hhh);
00230 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00231 }
00232 x += hhh;
00233 }
00234 loss *= hhh;
00235 if (loss < 0.) loss = 0.;
00236 return loss;
00237 }
00238
00239
00240
00241 G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
00242 G4double tkin,
00243 G4double Z,
00244 G4double cut)
00245 {
00246 G4double cross = 0.;
00247 SetCurrentElement(Z);
00248 G4double tmax = MaxSecondaryEnergy(particle, tkin);
00249 if (tmax <= cut) { return cross; }
00250
00251 G4double ak1=6.9 ;
00252 G4double ak2=1.0 ;
00253 G4double aaa = log(cut);
00254 G4double bbb = log(tmax);
00255 G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
00256 if(kkk > 8) { kkk = 8; }
00257 G4double hhh = (bbb-aaa)/G4double(kkk);
00258 G4double x = aaa;
00259
00260 for(G4int l=0; l<kkk; ++l)
00261 {
00262 for(G4int i=0; i<8; ++i)
00263 {
00264 G4double ep = exp(x + xgi[i]*hhh);
00265 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00266 }
00267 x += hhh;
00268 }
00269
00270 cross *= hhh;
00271 if(cross < 0.0) { cross = 0.0; }
00272 return cross;
00273 }
00274
00275 G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
00276 G4double tkin,
00277 G4double Z,
00278 G4double pairEnergy)
00279
00280
00281
00282 {
00283 G4double bbbtf= 183. ;
00284 G4double bbbh = 202.4 ;
00285 G4double g1tf = 1.95e-5 ;
00286 G4double g2tf = 5.3e-5 ;
00287 G4double g1h = 4.4e-5 ;
00288 G4double g2h = 4.8e-5 ;
00289
00290 G4double totalEnergy = tkin + particleMass;
00291 G4double residEnergy = totalEnergy - pairEnergy;
00292 G4double massratio = particleMass/electron_mass_c2 ;
00293 G4double massratio2 = massratio*massratio ;
00294 G4double cross = 0.;
00295
00296 SetCurrentElement(Z);
00297
00298 G4double c3 = 0.75*sqrte*particleMass;
00299 if (residEnergy <= c3*z13) { return cross; }
00300
00301 G4double c7 = 4.*electron_mass_c2;
00302 G4double c8 = 6.*particleMass*particleMass;
00303 G4double alf = c7/pairEnergy;
00304 G4double a3 = 1. - alf;
00305 if (a3 <= 0.) { return cross; }
00306
00307
00308 G4double bbb,g1,g2;
00309 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
00310 else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
00311
00312 G4double zeta = 0;
00313 G4double zeta1 = 0.073*log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26;
00314 if ( zeta1 > 0.)
00315 {
00316 G4double zeta2 = 0.058*log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14;
00317 zeta = zeta1/zeta2 ;
00318 }
00319
00320 G4double z2 = Z*(Z+zeta);
00321 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
00322 G4double a0 = totalEnergy*residEnergy;
00323 G4double a1 = pairEnergy*pairEnergy/a0;
00324 G4double bet = 0.5*a1;
00325 G4double xi0 = 0.25*massratio2*a1;
00326 G4double del = c8/a0;
00327
00328 G4double rta3 = sqrt(a3);
00329 G4double tmnexp = alf/(1. + rta3) + del*rta3;
00330 if(tmnexp >= 1.0) return cross;
00331
00332 G4double tmn = log(tmnexp);
00333 G4double sum = 0.;
00334
00335
00336 for (G4int i=0; i<8; ++i)
00337 {
00338 G4double a4 = exp(tmn*xgi[i]);
00339 G4double a5 = a4*(2.-a4) ;
00340 G4double a6 = 1.-a5 ;
00341 G4double a7 = 1.+a6 ;
00342 G4double a9 = 3.+a6 ;
00343 G4double xi = xi0*a5 ;
00344 G4double xii = 1./xi ;
00345 G4double xi1 = 1.+xi ;
00346 G4double screen = screen0*xi1/a5 ;
00347 G4double yeu = 5.-a6+4.*bet*a7 ;
00348 G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ;
00349 G4double ye1 = 1.+yeu/yed ;
00350 G4double ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
00351 G4double cre = 0.5*log(1.+2.25*z23*xi1*ye1/massratio2) ;
00352 G4double be;
00353
00354 if (xi <= 1.e3) be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
00355 else be = (3.-a6+a1*a7)/(2.*xi);
00356
00357 G4double fe = (ale-cre)*be;
00358 if ( fe < 0.) fe = 0. ;
00359
00360 G4double ymu = 4.+a6 +3.*bet*a7 ;
00361 G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ;
00362 G4double ym1 = 1.+ymu/ymd ;
00363 G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
00364 G4double a10,bm;
00365 if ( xi >= 1.e-3)
00366 {
00367 a10 = (1.+a1)*a5 ;
00368 bm = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
00369 } else {
00370 bm = (5.-a6+bet*a9)*(xi/2.);
00371 }
00372
00373 G4double fm = alm_crm*bm;
00374 if ( fm < 0.) fm = 0. ;
00375
00376 sum += wgi[i]*a4*(fe+fm/massratio2);
00377 }
00378
00379 cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
00380
00381 return cross;
00382 }
00383
00384
00385
00386 G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
00387 const G4ParticleDefinition*,
00388 G4double kineticEnergy,
00389 G4double Z, G4double,
00390 G4double cutEnergy,
00391 G4double maxEnergy)
00392 {
00393 G4double cross = 0.0;
00394 if (kineticEnergy <= lowestKinEnergy) { return cross; }
00395
00396 SetCurrentElement(Z);
00397
00398 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
00399 G4double tmax = std::min(maxEnergy, maxPairEnergy);
00400 G4double cut = std::max(cutEnergy, minPairEnergy);
00401 if (cut >= tmax) return cross;
00402
00403 cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
00404 if(tmax < kineticEnergy) {
00405 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
00406 }
00407 return cross;
00408 }
00409
00410
00411
00412 void G4MuPairProductionModel::MakeSamplingTables()
00413 {
00414 for (G4int iz=0; iz<nzdat; ++iz)
00415 {
00416 G4double Z = zdat[iz];
00417 SetCurrentElement(Z);
00418
00419 for (G4int it=0; it<ntdat; ++it) {
00420
00421 G4double kineticEnergy = tdat[it];
00422 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
00423
00424
00425 G4double xSec = 0.0 ;
00426
00427 if(maxPairEnergy > minPairEnergy) {
00428
00429 G4double y = ymin - 0.5*dy ;
00430 G4double yy = ymin - dy ;
00431 G4double x = exp(y);
00432 G4double fac = exp(dy);
00433 G4double dx = exp(yy)*(fac - 1.0);
00434
00435 G4double c = log(maxPairEnergy/minPairEnergy);
00436
00437 for (G4int i=0 ; i<nbiny; ++i) {
00438 y += dy ;
00439 if(c > 0.0) {
00440 x *= fac;
00441 dx*= fac;
00442 G4double ep = minPairEnergy*exp(c*x) ;
00443 xSec +=
00444 ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
00445 }
00446 ya[i] = y;
00447 proba[iz][it][i] = xSec;
00448 }
00449
00450 } else {
00451 for (G4int i=0 ; i<nbiny; ++i) {
00452 proba[iz][it][i] = xSec;
00453 }
00454 }
00455
00456 ya[nbiny]=ymax;
00457 proba[iz][it][nbiny] = xSec;
00458
00459 }
00460 }
00461 samplingTablesAreFilled = true;
00462 }
00463
00464
00465
00466 void
00467 G4MuPairProductionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
00468 const G4MaterialCutsCouple* couple,
00469 const G4DynamicParticle* aDynamicParticle,
00470 G4double tmin,
00471 G4double tmax)
00472 {
00473 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00474 G4double totalEnergy = kineticEnergy + particleMass;
00475 G4double totalMomentum =
00476 sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
00477
00478 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
00479
00480 G4int it;
00481 for(it=1; it<ntdat; ++it) { if(kineticEnergy <= tdat[it]) { break; } }
00482 if(it == ntdat) { --it; }
00483 G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
00484
00485
00486 const G4Element* anElement =
00487 SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
00488 SetCurrentElement(anElement->GetZ());
00489
00490
00491 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
00492 G4double maxEnergy = std::min(tmax, maxPairEnergy);
00493 G4double minEnergy = std::max(tmin, minPairEnergy);
00494
00495 if(minEnergy >= maxEnergy) { return; }
00496
00497
00498
00499
00500 G4double logmaxmin = log(maxPairEnergy/minPairEnergy);
00501
00502
00503 G4int iymin = 0;
00504 G4int iymax = nbiny-1;
00505 if( minEnergy > minPairEnergy)
00506 {
00507 G4double xc = log(minEnergy/minPairEnergy)/logmaxmin;
00508 iymin = (G4int)((log(xc) - ymin)/dy);
00509 if(iymin >= nbiny) iymin = nbiny-1;
00510 else if(iymin < 0) iymin = 0;
00511 xc = log(maxEnergy/minPairEnergy)/logmaxmin;
00512 iymax = (G4int)((log(xc) - ymin)/dy) + 1;
00513 if(iymax >= nbiny) iymax = nbiny-1;
00514 else if(iymax < 0) iymax = 0;
00515 }
00516
00517
00518 G4int iz, iy;
00519
00520 for(iz=1; iz<nzdat; ++iz) { if(currentZ <= zdat[iz]) { break; } }
00521 if(iz == nzdat) { --iz; }
00522
00523 G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
00524
00525 G4double pmin = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymin,currentZ);
00526 G4double pmax = InterpolatedIntegralCrossSection(dt,dz,iz,it,iymax,currentZ);
00527
00528 G4double p = pmin+G4UniformRand()*(pmax - pmin);
00529
00530
00531 G4double p1 = pmin;
00532 G4double p2 = pmin;
00533 for(iy=iymin+1; iy<=iymax; ++iy) {
00534 p1 = p2;
00535 p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
00536 if(p <= p2) { break; }
00537 }
00538
00539
00540 G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
00541
00542 G4double PairEnergy = minPairEnergy*exp( exp(y)*logmaxmin );
00543
00544 if(PairEnergy < minEnergy) { PairEnergy = minEnergy; }
00545 if(PairEnergy > maxEnergy) { PairEnergy = maxEnergy; }
00546
00547
00548 G4double rmax =
00549 (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
00550 *sqrt(1.-minPairEnergy/PairEnergy);
00551 G4double r = rmax * (-1.+2.*G4UniformRand()) ;
00552
00553
00554 G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
00555 G4double PositronEnergy = PairEnergy - ElectronEnergy;
00556
00557
00558
00559
00560 G4double gam = totalEnergy/particleMass;
00561 G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
00562 G4double gmax2= gmax*gmax;
00563 G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
00564
00565 G4double theta = sqrt(x/(1.0 - x))/gam;
00566 G4double sint = sin(theta);
00567 G4double phi = twopi * G4UniformRand() ;
00568 G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
00569
00570 G4ThreeVector gDirection(dirx, diry, dirz);
00571 gDirection.rotateUz(partDirection);
00572
00573
00574
00575
00576 G4DynamicParticle* aParticle1 =
00577 new G4DynamicParticle(theElectron, gDirection,
00578 ElectronEnergy - electron_mass_c2);
00579
00580
00581 G4DynamicParticle* aParticle2 =
00582 new G4DynamicParticle(thePositron, gDirection,
00583 PositronEnergy - electron_mass_c2);
00584
00585
00586 kineticEnergy -= (ElectronEnergy + PositronEnergy);
00587 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00588
00589 partDirection *= totalMomentum;
00590 partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
00591 partDirection = partDirection.unit();
00592 fParticleChange->SetProposedMomentumDirection(partDirection);
00593
00594
00595 vdp->push_back(aParticle1);
00596 vdp->push_back(aParticle2);
00597 }
00598
00599
00600
00601 const G4Element* G4MuPairProductionModel::SelectRandomAtom(
00602 G4double kinEnergy, G4double dt, G4int it,
00603 const G4MaterialCutsCouple* couple, G4double tmin)
00604 {
00605
00606
00607 const G4Material* material = couple->GetMaterial();
00608 size_t nElements = material->GetNumberOfElements();
00609 const G4ElementVector* theElementVector = material->GetElementVector();
00610 if (nElements == 1) { return (*theElementVector)[0]; }
00611
00612 if(nElements > nmaxElements) {
00613 nmaxElements = nElements;
00614 partialSum.resize(nmaxElements);
00615 }
00616
00617 const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
00618
00619 G4double sum = 0.0;
00620 G4double dl;
00621
00622 size_t i;
00623 for (i=0; i<nElements; ++i) {
00624 G4double Z = ((*theElementVector)[i])->GetZ();
00625 SetCurrentElement(Z);
00626 G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
00627 G4double minEnergy = std::max(tmin, minPairEnergy);
00628 dl = 0.0;
00629 if(minEnergy < maxPairEnergy) {
00630
00631 G4int iz;
00632 for(iz=1; iz<nzdat; ++iz) {if(Z <= zdat[iz]) { break; } }
00633 if(iz == nzdat) { --iz; }
00634 G4double dz = log(Z/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
00635
00636 G4double sigcut;
00637 if(minEnergy <= minPairEnergy)
00638 sigcut = 0.;
00639 else
00640 {
00641 G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
00642 G4int iy = (G4int)((log(xc) - ymin)/dy);
00643 if(iy < 0) { iy = 0; }
00644 if(iy >= nbiny) { iy = nbiny-1; }
00645 sigcut = InterpolatedIntegralCrossSection(dt,dz,iz,it,iy, Z);
00646 }
00647
00648 G4double sigtot = InterpolatedIntegralCrossSection(dt,dz,iz,it,nbiny,Z);
00649 dl = (sigtot - sigcut)*theAtomNumDensityVector[i];
00650 }
00651
00652 if(dl < 0.0) { dl = 0.0; }
00653 sum += dl;
00654 partialSum[i] = sum;
00655 }
00656
00657 G4double rval = G4UniformRand()*sum;
00658 for (i=0; i<nElements; ++i) {
00659 if(rval<=partialSum[i]) { return (*theElementVector)[i]; }
00660 }
00661
00662 return (*theElementVector)[nElements - 1];
00663
00664 }
00665
00666
00667
00668