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 #include "G4MuBremsstrahlungModel.hh"
00065 #include "G4PhysicalConstants.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4Gamma.hh"
00068 #include "G4MuonMinus.hh"
00069 #include "G4MuonPlus.hh"
00070 #include "Randomize.hh"
00071 #include "G4Material.hh"
00072 #include "G4Element.hh"
00073 #include "G4ElementVector.hh"
00074 #include "G4ProductionCutsTable.hh"
00075 #include "G4ParticleChangeForLoss.hh"
00076
00077
00078
00079
00080 using namespace std;
00081
00082 G4MuBremsstrahlungModel::G4MuBremsstrahlungModel(const G4ParticleDefinition* p,
00083 const G4String& nam)
00084 : G4VEmModel(nam),
00085 particle(0),
00086 sqrte(sqrt(exp(1.))),
00087 bh(202.4),
00088 bh1(446.),
00089 btf(183.),
00090 btf1(1429.),
00091 fParticleChange(0),
00092 lowestKinEnergy(1.0*GeV),
00093 minThreshold(1.0*keV)
00094 {
00095 theGamma = G4Gamma::Gamma();
00096 nist = G4NistManager::Instance();
00097
00098 mass = rmass = cc = coeff = 1.0;
00099
00100 fDN[0] = 0.0;
00101 for(G4int i=1; i<93; ++i) {
00102 G4double dn = 1.54*nist->GetA27(i);
00103 fDN[i] = dn;
00104 if(1 < i) {
00105 fDN[i] /= std::pow(dn, 1./G4double(i));
00106 }
00107 }
00108
00109 if(p) { SetParticle(p); }
00110 }
00111
00112
00113
00114 G4MuBremsstrahlungModel::~G4MuBremsstrahlungModel()
00115 {
00116 size_t n = partialSumSigma.size();
00117 if(n > 0) {
00118 for(size_t i=0; i<n; i++) {
00119 delete partialSumSigma[i];
00120 }
00121 }
00122 }
00123
00124
00125
00126 G4double G4MuBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
00127 const G4MaterialCutsCouple*)
00128 {
00129 return minThreshold;
00130 }
00131
00132
00133
00134 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
00135 const G4DataVector& cuts)
00136 {
00137 if(p) { SetParticle(p); }
00138
00139
00140 G4double fixedEnergy = 0.5*HighEnergyLimit();
00141
00142 const G4ProductionCutsTable* theCoupleTable=
00143 G4ProductionCutsTable::GetProductionCutsTable();
00144 if(theCoupleTable) {
00145 G4int numOfCouples = theCoupleTable->GetTableSize();
00146
00147 G4int nn = partialSumSigma.size();
00148 G4int nc = cuts.size();
00149
00150
00151 if(nn == numOfCouples) { return; }
00152
00153
00154 if(nn > 0) {
00155 for (G4int ii=0; ii<nn; ii++){
00156 G4DataVector* a = partialSumSigma[ii];
00157 if ( a ) { delete a; }
00158 }
00159 partialSumSigma.clear();
00160 }
00161
00162 if (numOfCouples>0) {
00163 for (G4int i=0; i<numOfCouples; i++) {
00164 G4double cute = DBL_MAX;
00165
00166
00167 if(i < nc) { cute = cuts[i]; }
00168
00169 const G4MaterialCutsCouple* couple =
00170 theCoupleTable->GetMaterialCutsCouple(i);
00171 const G4Material* material = couple->GetMaterial();
00172 G4DataVector* dv = ComputePartialSumSigma(material,fixedEnergy,cute);
00173 partialSumSigma.push_back(dv);
00174 }
00175 }
00176 }
00177
00178
00179 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00180 }
00181
00182
00183
00184 G4double G4MuBremsstrahlungModel::ComputeDEDXPerVolume(
00185 const G4Material* material,
00186 const G4ParticleDefinition*,
00187 G4double kineticEnergy,
00188 G4double cutEnergy)
00189 {
00190 G4double dedx = 0.0;
00191 if (kineticEnergy <= lowestKinEnergy) return dedx;
00192
00193 G4double tmax = kineticEnergy;
00194 G4double cut = std::min(cutEnergy,tmax);
00195 if(cut < minThreshold) cut = minThreshold;
00196
00197 const G4ElementVector* theElementVector = material->GetElementVector();
00198 const G4double* theAtomicNumDensityVector =
00199 material->GetAtomicNumDensityVector();
00200
00201
00202 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00203
00204 G4double loss =
00205 ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
00206
00207 dedx += loss*theAtomicNumDensityVector[i];
00208 }
00209
00210 if(dedx < 0.) dedx = 0.;
00211 return dedx;
00212 }
00213
00214
00215
00216 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
00217 G4double tkin, G4double cut)
00218 {
00219 G4double totalEnergy = mass + tkin;
00220 G4double ak1 = 0.05;
00221 G4int k2=5;
00222 G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
00223 G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
00224 G4double loss = 0.;
00225
00226 G4double vcut = cut/totalEnergy;
00227 G4double vmax = tkin/totalEnergy;
00228
00229 G4double aaa = 0.;
00230 G4double bbb = vcut;
00231 if(vcut>vmax) bbb=vmax ;
00232 G4int kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
00233 if (kkk < 1) { kkk = 1; }
00234 G4double hhh=(bbb-aaa)/float(kkk) ;
00235
00236 G4double aa = aaa;
00237 for(G4int l=0; l<kkk; l++)
00238 {
00239 for(G4int i=0; i<6; i++)
00240 {
00241 G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
00242 loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00243 }
00244 aa += hhh;
00245 }
00246
00247 loss *=hhh*totalEnergy ;
00248
00249 return loss;
00250 }
00251
00252
00253
00254 G4double G4MuBremsstrahlungModel::ComputeMicroscopicCrossSection(
00255 G4double tkin,
00256 G4double Z,
00257 G4double cut)
00258 {
00259 G4double totalEnergy = tkin + mass;
00260 G4double ak1 = 2.3;
00261 G4int k2 = 4;
00262 G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
00263 G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
00264 G4double cross = 0.;
00265
00266 if(cut >= tkin) return cross;
00267
00268 G4double vcut = cut/totalEnergy;
00269 G4double vmax = tkin/totalEnergy;
00270
00271 G4double aaa = log(vcut);
00272 G4double bbb = log(vmax);
00273 G4int kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
00274 G4double hhh = (bbb-aaa)/G4double(kkk);
00275
00276 G4double aa = aaa;
00277
00278 for(G4int l=0; l<kkk; l++)
00279 {
00280 for(G4int i=0; i<6; i++)
00281 {
00282 G4double ep = exp(aa + xgi[i]*hhh)*totalEnergy;
00283 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00284 }
00285 aa += hhh;
00286 }
00287
00288 cross *=hhh;
00289
00290
00291
00292 return cross;
00293 }
00294
00295
00296
00297 G4double G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection(
00298 G4double tkin,
00299 G4double Z,
00300 G4double gammaEnergy)
00301
00302 {
00303 G4double dxsection = 0.;
00304
00305 if( gammaEnergy > tkin) return dxsection ;
00306
00307 G4double E = tkin + mass ;
00308 G4double v = gammaEnergy/E ;
00309 G4double delta = 0.5*mass*mass*v/(E-gammaEnergy) ;
00310 G4double rab0=delta*sqrte ;
00311
00312 G4int iz = G4int(Z);
00313 if(iz < 1) iz = 1;
00314 else if(iz > 92) iz = 92;
00315
00316 G4double z13 = 1.0/nist->GetZ13(iz);
00317 G4double dnstar = fDN[iz];
00318
00319 G4double b,b1;
00320
00321 if(1 == iz) {
00322 b = bh;
00323 b1 = bh1;
00324 } else {
00325 b = btf;
00326 b1 = btf1;
00327 }
00328
00329
00330 G4double rab1=b*z13;
00331 G4double fn=log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
00332 (mass+delta*(dnstar*sqrte-2.))) ;
00333 if(fn <0.) fn = 0. ;
00334
00335 G4double epmax1=E/(1.+0.5*mass*rmass/E) ;
00336 G4double fe=0.;
00337 if(gammaEnergy<epmax1)
00338 {
00339 G4double rab2=b1*z13*z13 ;
00340 fe=log(rab2*mass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
00341 (electron_mass_c2+rab0*rab2))) ;
00342 if(fe<0.) fe=0. ;
00343 }
00344
00345 dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn*Z + fe)/gammaEnergy;
00346
00347 return dxsection;
00348 }
00349
00350
00351
00352 G4double G4MuBremsstrahlungModel::ComputeCrossSectionPerAtom(
00353 const G4ParticleDefinition*,
00354 G4double kineticEnergy,
00355 G4double Z, G4double,
00356 G4double cutEnergy,
00357 G4double maxEnergy)
00358 {
00359 G4double cross = 0.0;
00360 if (kineticEnergy <= lowestKinEnergy) return cross;
00361 G4double tmax = std::min(maxEnergy, kineticEnergy);
00362 G4double cut = std::min(cutEnergy, kineticEnergy);
00363 if(cut < minThreshold) cut = minThreshold;
00364 if (cut >= tmax) return cross;
00365
00366 cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
00367 if(tmax < kineticEnergy) {
00368 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
00369 }
00370 return cross;
00371 }
00372
00373
00374
00375 G4DataVector* G4MuBremsstrahlungModel::ComputePartialSumSigma(
00376 const G4Material* material,
00377 G4double kineticEnergy,
00378 G4double cut)
00379
00380
00381
00382
00383 {
00384 G4int nElements = material->GetNumberOfElements();
00385 const G4ElementVector* theElementVector = material->GetElementVector();
00386 const G4double* theAtomNumDensityVector =
00387 material->GetAtomicNumDensityVector();
00388
00389 G4DataVector* dv = new G4DataVector();
00390
00391 G4double cross = 0.0;
00392
00393 for (G4int i=0; i<nElements; i++ ) {
00394 cross += theAtomNumDensityVector[i]
00395 * ComputeMicroscopicCrossSection(kineticEnergy,
00396 (*theElementVector)[i]->GetZ(), cut);
00397 dv->push_back(cross);
00398 }
00399 return dv;
00400 }
00401
00402
00403
00404 void G4MuBremsstrahlungModel::SampleSecondaries(
00405 std::vector<G4DynamicParticle*>* vdp,
00406 const G4MaterialCutsCouple* couple,
00407 const G4DynamicParticle* dp,
00408 G4double minEnergy,
00409 G4double maxEnergy)
00410 {
00411 G4double kineticEnergy = dp->GetKineticEnergy();
00412
00413 G4double tmax = std::min(kineticEnergy, maxEnergy);
00414 G4double tmin = std::min(kineticEnergy, minEnergy);
00415 if(tmin < minThreshold) tmin = minThreshold;
00416 if(tmin >= tmax) return;
00417
00418
00419
00420 G4ParticleMomentum partDirection = dp->GetMomentumDirection();
00421
00422
00423 const G4Element* anElement = SelectRandomAtom(couple);
00424 G4double Z = anElement->GetZ();
00425
00426 G4double totalEnergy = kineticEnergy + mass;
00427 G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
00428
00429 G4double func1 = tmin*
00430 ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
00431
00432 G4double lnepksi, epksi;
00433 G4double func2;
00434
00435 G4double xmin = log(tmin/MeV);
00436 G4double xmax = log(kineticEnergy/tmin);
00437
00438 do {
00439 lnepksi = xmin + G4UniformRand()*xmax;
00440 epksi = MeV*exp(lnepksi);
00441 func2 = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
00442
00443 } while(func2 < func1*G4UniformRand());
00444
00445 G4double gEnergy = epksi;
00446
00447
00448
00449 G4double gam = totalEnergy/mass;
00450 G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
00451 G4double rmax2= rmax*rmax;
00452 G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
00453
00454 G4double theta = sqrt(x/(1.0 - x))/gam;
00455 G4double sint = sin(theta);
00456 G4double phi = twopi * G4UniformRand() ;
00457 G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
00458
00459 G4ThreeVector gDirection(dirx, diry, dirz);
00460 gDirection.rotateUz(partDirection);
00461
00462 partDirection *= totalMomentum;
00463 partDirection -= gEnergy*gDirection;
00464 partDirection = partDirection.unit();
00465
00466
00467 kineticEnergy -= gEnergy;
00468 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00469 fParticleChange->SetProposedMomentumDirection(partDirection);
00470
00471
00472 G4DynamicParticle* aGamma =
00473 new G4DynamicParticle(theGamma,gDirection,gEnergy);
00474 vdp->push_back(aGamma);
00475 }
00476
00477
00478
00479 const G4Element* G4MuBremsstrahlungModel::SelectRandomAtom(
00480 const G4MaterialCutsCouple* couple) const
00481 {
00482
00483
00484 const G4Material* material = couple->GetMaterial();
00485 G4int nElements = material->GetNumberOfElements();
00486 const G4ElementVector* theElementVector = material->GetElementVector();
00487 if(1 == nElements) { return (*theElementVector)[0]; }
00488 else if(1 > nElements) { return 0; }
00489
00490 G4DataVector* dv = partialSumSigma[couple->GetIndex()];
00491 G4double rval = G4UniformRand()*((*dv)[nElements-1]);
00492 for (G4int i=0; i<nElements; i++) {
00493 if (rval <= (*dv)[i]) { return (*theElementVector)[i]; }
00494 }
00495 return (*theElementVector)[nElements-1];
00496 }
00497
00498