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 #include "G4EnergyLossForExtrapolator.hh"
00051 #include "G4PhysicalConstants.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4PhysicsLogVector.hh"
00054 #include "G4ParticleDefinition.hh"
00055 #include "G4Material.hh"
00056 #include "G4MaterialCutsCouple.hh"
00057 #include "G4Electron.hh"
00058 #include "G4Positron.hh"
00059 #include "G4Proton.hh"
00060 #include "G4MuonPlus.hh"
00061 #include "G4MuonMinus.hh"
00062 #include "G4ParticleTable.hh"
00063 #include "G4LossTableBuilder.hh"
00064 #include "G4MollerBhabhaModel.hh"
00065 #include "G4BetheBlochModel.hh"
00066 #include "G4eBremsstrahlungModel.hh"
00067 #include "G4MuPairProductionModel.hh"
00068 #include "G4MuBremsstrahlungModel.hh"
00069 #include "G4ProductionCuts.hh"
00070 #include "G4LossTableManager.hh"
00071 #include "G4WentzelVIModel.hh"
00072
00073
00074
00075 G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
00076 :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false)
00077 {
00078 currentParticle = 0;
00079 currentMaterial = 0;
00080
00081 linLossLimit = 0.01;
00082 emin = 1.*MeV;
00083 emax = 10.*TeV;
00084 nbins = 70;
00085
00086 nmat = index = 0;
00087 cuts = 0;
00088
00089 mass = charge2 = electronDensity = radLength = bg2 = beta2
00090 = kineticEnergy = tmax = 0;
00091 gam = 1.0;
00092
00093 dedxElectron = dedxPositron = dedxProton = rangeElectron
00094 = rangePositron = rangeProton = invRangeElectron = invRangePositron
00095 = invRangeProton = mscElectron = dedxMuon = rangeMuon = invRangeMuon = 0;
00096 cuts = 0;
00097 electron = positron = proton = muonPlus = muonMinus = 0;
00098 }
00099
00100
00101
00102 G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator()
00103 {
00104 for(G4int i=0; i<nmat; i++) {delete couples[i];}
00105 delete dedxElectron;
00106 delete dedxPositron;
00107 delete dedxProton;
00108 delete dedxMuon;
00109 delete rangeElectron;
00110 delete rangePositron;
00111 delete rangeProton;
00112 delete rangeMuon;
00113 delete invRangeElectron;
00114 delete invRangePositron;
00115 delete invRangeProton;
00116 delete invRangeMuon;
00117 delete mscElectron;
00118 delete cuts;
00119 }
00120
00121
00122
00123 G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
00124 G4double stepLength,
00125 const G4Material* mat,
00126 const G4ParticleDefinition* part)
00127 {
00128 if(!isInitialised) Initialisation();
00129 G4double kinEnergyFinal = kinEnergy;
00130 if(SetupKinematics(part, mat, kinEnergy)) {
00131 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
00132 G4double r = ComputeRange(kinEnergy,part);
00133 if(r <= step) {
00134 kinEnergyFinal = 0.0;
00135 } else if(step < linLossLimit*r) {
00136 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
00137 } else {
00138 G4double r1 = r - step;
00139 kinEnergyFinal = ComputeEnergy(r1,part);
00140 }
00141 }
00142 return kinEnergyFinal;
00143 }
00144
00145
00146
00147 G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
00148 G4double stepLength,
00149 const G4Material* mat,
00150 const G4ParticleDefinition* part)
00151 {
00152 if(!isInitialised) Initialisation();
00153 G4double kinEnergyFinal = kinEnergy;
00154
00155 if(SetupKinematics(part, mat, kinEnergy)) {
00156 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
00157 G4double r = ComputeRange(kinEnergy,part);
00158
00159 if(step < linLossLimit*r) {
00160 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
00161 } else {
00162 G4double r1 = r + step;
00163 kinEnergyFinal = ComputeEnergy(r1,part);
00164 }
00165 }
00166 return kinEnergyFinal;
00167 }
00168
00169
00170
00171 G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
00172 G4double stepLength,
00173 const G4Material* mat,
00174 const G4ParticleDefinition* part)
00175 {
00176 G4double res = stepLength;
00177 if(!isInitialised) Initialisation();
00178 if(SetupKinematics(part, mat, kinEnergy)) {
00179 if(part == electron || part == positron) {
00180 G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
00181 if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
00182 else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
00183 else res = ComputeRange(kinEnergy,part);
00184
00185 } else {
00186 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
00187 }
00188 }
00189 return res;
00190 }
00191
00192
00193
00194 G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
00195 const G4Material* mat,
00196 G4double kinEnergy)
00197 {
00198 if(!part || !mat || kinEnergy < keV) return false;
00199 if(!isInitialised) Initialisation();
00200 G4bool flag = false;
00201 if(part != currentParticle) {
00202 flag = true;
00203 currentParticle = part;
00204 mass = part->GetPDGMass();
00205 G4double q = part->GetPDGCharge()/eplus;
00206 charge2 = q*q;
00207 }
00208 if(mat != currentMaterial) {
00209 G4int i = mat->GetIndex();
00210 if(i >= nmat) {
00211 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
00212 << i << " is out of table - NO extrapolation" << G4endl;
00213 } else {
00214 flag = true;
00215 currentMaterial = mat;
00216 electronDensity = mat->GetElectronDensity();
00217 radLength = mat->GetRadlen();
00218 index = i;
00219 }
00220 }
00221 if(flag || kinEnergy != kineticEnergy) {
00222 kineticEnergy = kinEnergy;
00223 G4double tau = kinEnergy/mass;
00224
00225 gam = tau + 1.0;
00226 bg2 = tau * (tau + 2.0);
00227 beta2 = bg2/(gam*gam);
00228 tmax = kinEnergy;
00229 if(part == electron) tmax *= 0.5;
00230 else if(part != positron) {
00231 G4double r = electron_mass_c2/mass;
00232 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
00233 }
00234 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
00235 }
00236 return true;
00237 }
00238
00239
00240
00241 void G4EnergyLossForExtrapolator::Initialisation()
00242 {
00243 isInitialised = true;
00244 if(verbose>1) {
00245 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
00246 }
00247 currentParticle = 0;
00248 currentMaterial = 0;
00249 kineticEnergy = 0.0;
00250
00251 electron = G4Electron::Electron();
00252 positron = G4Positron::Positron();
00253 proton = G4Proton::Proton();
00254 muonPlus = G4MuonPlus::MuonPlus();
00255 muonMinus= G4MuonMinus::MuonMinus();
00256
00257 currentParticleName = "";
00258
00259 nmat = G4Material::GetNumberOfMaterials();
00260 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00261 cuts = new G4ProductionCuts();
00262
00263 const G4MaterialCutsCouple* couple;
00264 for(G4int i=0; i<nmat; i++) {
00265 couple = new G4MaterialCutsCouple((*mtable)[i],cuts);
00266 couples.push_back(couple);
00267 }
00268
00269 dedxElectron = PrepareTable();
00270 dedxPositron = PrepareTable();
00271 dedxMuon = PrepareTable();
00272 dedxProton = PrepareTable();
00273 rangeElectron = PrepareTable();
00274 rangePositron = PrepareTable();
00275 rangeMuon = PrepareTable();
00276 rangeProton = PrepareTable();
00277 invRangeElectron = PrepareTable();
00278 invRangePositron = PrepareTable();
00279 invRangeMuon = PrepareTable();
00280 invRangeProton = PrepareTable();
00281 mscElectron = PrepareTable();
00282
00283 G4LossTableBuilder builder;
00284
00285 if(verbose>1)
00286 G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl;
00287
00288 ComputeElectronDEDX(electron, dedxElectron);
00289 builder.BuildRangeTable(dedxElectron,rangeElectron);
00290 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
00291
00292 if(verbose>1)
00293 G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl;
00294
00295 ComputeElectronDEDX(positron, dedxPositron);
00296 builder.BuildRangeTable(dedxPositron, rangePositron);
00297 builder.BuildInverseRangeTable(rangePositron, invRangePositron);
00298
00299 if(verbose>1)
00300 G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl;
00301
00302 ComputeMuonDEDX(muonPlus, dedxMuon);
00303 builder.BuildRangeTable(dedxMuon, rangeMuon);
00304 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
00305
00306 if(verbose>1)
00307 G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl;
00308
00309 ComputeProtonDEDX(proton, dedxProton);
00310 builder.BuildRangeTable(dedxProton, rangeProton);
00311 builder.BuildInverseRangeTable(rangeProton, invRangeProton);
00312
00313 ComputeTrasportXS(electron, mscElectron);
00314 }
00315
00316
00317
00318 G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable()
00319 {
00320 G4PhysicsTable* table = new G4PhysicsTable();
00321
00322 for(G4int i=0; i<nmat; i++) {
00323
00324 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
00325 v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
00326 table->push_back(v);
00327 }
00328 return table;
00329 }
00330
00331
00332
00333 const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
00334 {
00335 const G4ParticleDefinition* p = 0;
00336 if(name != currentParticleName) {
00337 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
00338 if(!p) {
00339 G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find "
00340 << name << G4endl;
00341 }
00342 } else {
00343 p = currentParticle;
00344 }
00345 return p;
00346 }
00347
00348
00349
00350 G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy,
00351 const G4ParticleDefinition* part)
00352 {
00353 G4double x = 0.0;
00354 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron);
00355 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
00356 else if(part == muonPlus || part == muonMinus) {
00357 x = ComputeValue(kinEnergy, dedxMuon);
00358 } else {
00359 G4double e = kinEnergy*proton_mass_c2/mass;
00360 x = ComputeValue(e, dedxProton)*charge2;
00361 }
00362 return x;
00363 }
00364
00365
00366
00367 G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy,
00368 const G4ParticleDefinition* part)
00369 {
00370 G4double x = 0.0;
00371 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron);
00372 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
00373 else if(part == muonPlus || part == muonMinus)
00374 x = ComputeValue(kinEnergy, rangeMuon);
00375 else {
00376 G4double massratio = proton_mass_c2/mass;
00377 G4double e = kinEnergy*massratio;
00378 x = ComputeValue(e, rangeProton)/(charge2*massratio);
00379 }
00380 return x;
00381 }
00382
00383
00384
00385 G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
00386 const G4ParticleDefinition* part)
00387 {
00388 G4double x = 0.0;
00389 if(part == electron) x = ComputeValue(range, invRangeElectron);
00390 else if(part == positron) x = ComputeValue(range, invRangePositron);
00391 else if(part == muonPlus || part == muonMinus)
00392 x = ComputeValue(range, invRangeMuon);
00393 else {
00394 G4double massratio = proton_mass_c2/mass;
00395 G4double r = range*massratio*charge2;
00396 x = ComputeValue(r, invRangeProton)/massratio;
00397 }
00398 return x;
00399 }
00400
00401
00402
00403 void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part,
00404 G4PhysicsTable* table)
00405 {
00406 G4DataVector v;
00407 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
00408 G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
00409 ioni->Initialise(part, v);
00410 brem->Initialise(part, v);
00411
00412 mass = electron_mass_c2;
00413 charge2 = 1.0;
00414 currentParticle = part;
00415
00416 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00417 if(0<verbose) {
00418 G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for "
00419 << part->GetParticleName()
00420 << G4endl;
00421 }
00422 for(G4int i=0; i<nmat; i++) {
00423
00424 const G4Material* mat = (*mtable)[i];
00425 if(1<verbose)
00426 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
00427 const G4MaterialCutsCouple* couple = couples[i];
00428 G4PhysicsVector* aVector = (*table)[i];
00429
00430 for(G4int j=0; j<=nbins; j++) {
00431
00432 G4double e = aVector->Energy(j);
00433 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e);
00434 if(1<verbose) {
00435 G4cout << "j= " << j
00436 << " e(MeV)= " << e/MeV
00437 << " dedx(Mev/cm)= " << dedx*cm/MeV
00438 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
00439 }
00440 aVector->PutValue(j,dedx);
00441 }
00442 }
00443 delete ioni;
00444 delete brem;
00445 }
00446
00447
00448
00449 void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
00450 G4PhysicsTable* table)
00451 {
00452 G4DataVector v;
00453 G4BetheBlochModel* ioni = new G4BetheBlochModel();
00454 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
00455 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
00456 ioni->Initialise(part, v);
00457 pair->Initialise(part, v);
00458 brem->Initialise(part, v);
00459
00460 mass = part->GetPDGMass();
00461 charge2 = 1.0;
00462 currentParticle = part;
00463
00464 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00465
00466 if(0<verbose) {
00467 G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName()
00468 << G4endl;
00469 }
00470
00471 for(G4int i=0; i<nmat; i++) {
00472
00473 const G4Material* mat = (*mtable)[i];
00474 if(1<verbose)
00475 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
00476 const G4MaterialCutsCouple* couple = couples[i];
00477 G4PhysicsVector* aVector = (*table)[i];
00478 for(G4int j=0; j<=nbins; j++) {
00479
00480 G4double e = aVector->Energy(j);
00481 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
00482 pair->ComputeDEDX(couple,part,e,e) +
00483 brem->ComputeDEDX(couple,part,e,e);
00484 aVector->PutValue(j,dedx);
00485 if(1<verbose) {
00486 G4cout << "j= " << j
00487 << " e(MeV)= " << e/MeV
00488 << " dedx(Mev/cm)= " << dedx*cm/MeV
00489 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
00490 }
00491 }
00492 }
00493 delete ioni;
00494 }
00495
00496
00497
00498 void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
00499 G4PhysicsTable* table)
00500 {
00501 G4DataVector v;
00502 G4BetheBlochModel* ioni = new G4BetheBlochModel();
00503 ioni->Initialise(part, v);
00504
00505 mass = part->GetPDGMass();
00506 charge2 = 1.0;
00507 currentParticle = part;
00508
00509 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00510
00511 if(0<verbose) {
00512 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
00513 << G4endl;
00514 }
00515
00516 for(G4int i=0; i<nmat; i++) {
00517
00518 const G4Material* mat = (*mtable)[i];
00519 if(1<verbose)
00520 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
00521 const G4MaterialCutsCouple* couple = couples[i];
00522 G4PhysicsVector* aVector = (*table)[i];
00523 for(G4int j=0; j<=nbins; j++) {
00524
00525 G4double e = aVector->Energy(j);
00526 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
00527 aVector->PutValue(j,dedx);
00528 if(1<verbose) {
00529 G4cout << "j= " << j
00530 << " e(MeV)= " << e/MeV
00531 << " dedx(Mev/cm)= " << dedx*cm/MeV
00532 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl;
00533 }
00534 }
00535 }
00536 delete ioni;
00537 }
00538
00539
00540
00541 void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
00542 G4PhysicsTable* table)
00543 {
00544 G4DataVector v;
00545 G4WentzelVIModel* msc = new G4WentzelVIModel();
00546 msc->SetPolarAngleLimit(CLHEP::pi);
00547 msc->Initialise(part, v);
00548
00549 mass = part->GetPDGMass();
00550 charge2 = 1.0;
00551 currentParticle = part;
00552
00553 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
00554
00555 if(0<verbose) {
00556 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
00557 << G4endl;
00558 }
00559
00560 for(G4int i=0; i<nmat; i++) {
00561
00562 const G4Material* mat = (*mtable)[i];
00563 msc->SetCurrentCouple(couples[i]);
00564 if(1<verbose)
00565 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
00566 G4PhysicsVector* aVector = (*table)[i];
00567 for(G4int j=0; j<=nbins; j++) {
00568
00569 G4double e = aVector->Energy(j);
00570 G4double xs = msc->CrossSectionPerVolume(mat,part,e);
00571 aVector->PutValue(j,xs);
00572 if(1<verbose) {
00573 G4cout << "j= " << j << " e(MeV)= " << e/MeV
00574 << " xs(1/mm)= " << xs*mm << G4endl;
00575 }
00576 }
00577 }
00578 delete msc;
00579 }
00580
00581
00582