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 #include "G4PenelopeBremsstrahlungModel.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4PenelopeBremsstrahlungFS.hh"
00043 #include "G4PenelopeBremsstrahlungAngular.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4MaterialCutsCouple.hh"
00046 #include "G4ProductionCutsTable.hh"
00047 #include "G4DynamicParticle.hh"
00048 #include "G4Gamma.hh"
00049 #include "G4Electron.hh"
00050 #include "G4Positron.hh"
00051 #include "G4PenelopeOscillatorManager.hh"
00052 #include "G4PenelopeCrossSection.hh"
00053 #include "G4PhysicsFreeVector.hh"
00054 #include "G4PhysicsLogVector.hh"
00055 #include "G4PhysicsTable.hh"
00056
00057
00058
00059 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
00060 const G4String& nam)
00061 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),energyGrid(0),
00062 XSTableElectron(0),XSTablePositron(0)
00063
00064 {
00065 fIntrinsicLowEnergyLimit = 100.0*eV;
00066 fIntrinsicHighEnergyLimit = 100.0*GeV;
00067 nBins = 200;
00068
00069 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00070
00071 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
00072
00073 verboseLevel= 0;
00074
00075
00076
00077
00078
00079
00080
00081
00082 fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS();
00083 fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
00084
00085
00086 SetDeexcitationFlag(true);
00087
00088 }
00089
00090
00091
00092 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
00093 {
00094 ClearTables();
00095 if (fPenelopeFSHelper)
00096 delete fPenelopeFSHelper;
00097 if (fPenelopeAngular)
00098 delete fPenelopeAngular;
00099 }
00100
00101
00102
00103 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition*,
00104 const G4DataVector&)
00105 {
00106 if (verboseLevel > 3)
00107 G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
00108
00109
00110 ClearTables();
00111
00112
00113 if (fPenelopeAngular)
00114 fPenelopeAngular->Initialize();
00115
00116
00117
00118 nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
00119 nBins = std::max(nBins,(size_t)100);
00120 energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
00121 HighEnergyLimit(),
00122 nBins-1);
00123
00124
00125 XSTableElectron = new
00126 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
00127 XSTablePositron = new
00128 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
00129
00130 if (verboseLevel > 2) {
00131 G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
00132 << "Energy range: "
00133 << LowEnergyLimit() / keV << " keV - "
00134 << HighEnergyLimit() / GeV << " GeV."
00135 << G4endl;
00136 }
00137
00138 if(isInitialised) return;
00139 fParticleChange = GetParticleChangeForLoss();
00140 isInitialised = true;
00141 }
00142
00143
00144
00145 G4double G4PenelopeBremsstrahlungModel::CrossSectionPerVolume(const G4Material* material,
00146 const G4ParticleDefinition* theParticle,
00147 G4double energy,
00148 G4double cutEnergy,
00149 G4double)
00150 {
00151
00152 if (verboseLevel > 3)
00153 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
00154
00155 SetupForMaterial(theParticle, material, energy);
00156
00157 G4double crossPerMolecule = 0.;
00158
00159 G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
00160 cutEnergy);
00161
00162 if (theXS)
00163 crossPerMolecule = theXS->GetHardCrossSection(energy);
00164
00165 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00166 G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
00167
00168 if (verboseLevel > 3)
00169 G4cout << "Material " << material->GetName() << " has " << atPerMol <<
00170 "atoms per molecule" << G4endl;
00171
00172 G4double moleculeDensity = 0.;
00173 if (atPerMol)
00174 moleculeDensity = atomDensity/atPerMol;
00175
00176 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
00177
00178 if (verboseLevel > 2)
00179 {
00180 G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
00181 G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
00182 energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
00183 }
00184
00185 return crossPerVolume;
00186 }
00187
00188
00189
00190
00191
00192
00193 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00194 G4double,
00195 G4double,
00196 G4double,
00197 G4double,
00198 G4double)
00199 {
00200 G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
00201 G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
00202 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
00203 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
00204 return 0;
00205 }
00206
00207
00208
00209 G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
00210 const G4ParticleDefinition* theParticle,
00211 G4double kineticEnergy,
00212 G4double cutEnergy)
00213 {
00214 if (verboseLevel > 3)
00215 G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
00216
00217 G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
00218 cutEnergy);
00219 G4double sPowerPerMolecule = 0.0;
00220 if (theXS)
00221 sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
00222
00223
00224 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
00225 G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
00226
00227 G4double moleculeDensity = 0.;
00228 if (atPerMol)
00229 moleculeDensity = atomDensity/atPerMol;
00230
00231 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
00232
00233 if (verboseLevel > 2)
00234 {
00235 G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
00236 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
00237 kineticEnergy/keV << " keV = " <<
00238 sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
00239 }
00240 return sPowerPerVolume;
00241 }
00242
00243
00244
00245 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
00246 const G4MaterialCutsCouple* couple,
00247 const G4DynamicParticle* aDynamicParticle,
00248 G4double cutG,
00249 G4double)
00250 {
00251 if (verboseLevel > 3)
00252 G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
00253
00254 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00255 const G4Material* material = couple->GetMaterial();
00256
00257 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
00258 {
00259 fParticleChange->SetProposedKineticEnergy(0.);
00260 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
00261 return ;
00262 }
00263
00264 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
00265
00266 G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
00267
00268
00269 if (kineticEnergy < cutG)
00270 return;
00271
00272 if (verboseLevel > 3)
00273 G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " <<
00274 "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
00275
00276
00277 G4double gammaEnergy =
00278 fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
00279
00280 if (verboseLevel > 3)
00281 G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
00282
00283
00284
00285 G4ThreeVector gammaDirection1 =
00286 fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
00287
00288 if (verboseLevel > 3)
00289 G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
00290
00291 G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
00292 if (residualPrimaryEnergy < 0)
00293 {
00294
00295 gammaEnergy += residualPrimaryEnergy;
00296 residualPrimaryEnergy = 0.0;
00297 }
00298
00299
00300 G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
00301 particleDirection1 = particleDirection1.unit();
00302
00303
00304 if (residualPrimaryEnergy > 0.)
00305 {
00306 fParticleChange->ProposeMomentumDirection(particleDirection1);
00307 fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
00308 }
00309 else
00310 {
00311 fParticleChange->SetProposedKineticEnergy(0.);
00312 }
00313
00314
00315 G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
00316 gammaDirection1,
00317 gammaEnergy);
00318 fvect->push_back(theGamma);
00319
00320 if (verboseLevel > 1)
00321 {
00322 G4cout << "-----------------------------------------------------------" << G4endl;
00323 G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
00324 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
00325 G4cout << "-----------------------------------------------------------" << G4endl;
00326 G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
00327 G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
00328 G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
00329 << " keV" << G4endl;
00330 G4cout << "-----------------------------------------------------------" << G4endl;
00331 }
00332
00333 if (verboseLevel > 0)
00334 {
00335 G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
00336 if (energyDiff > 0.05*keV)
00337 G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: "
00338 <<
00339 (residualPrimaryEnergy+gammaEnergy)/keV <<
00340 " keV (final) vs. " <<
00341 kineticEnergy/keV << " keV (initial)" << G4endl;
00342 }
00343 return;
00344 }
00345
00346
00347
00348 void G4PenelopeBremsstrahlungModel::ClearTables()
00349 {
00350 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>::iterator i;
00351 if (XSTableElectron)
00352 {
00353 for (i=XSTableElectron->begin(); i != XSTableElectron->end(); i++)
00354 {
00355 G4PenelopeCrossSection* tab = i->second;
00356 delete tab;
00357 }
00358 delete XSTableElectron;
00359 XSTableElectron = 0;
00360 }
00361 if (XSTablePositron)
00362 {
00363 for (i=XSTablePositron->begin(); i != XSTablePositron->end(); i++)
00364 {
00365 G4PenelopeCrossSection* tab = i->second;
00366 delete tab;
00367 }
00368 delete XSTablePositron;
00369 XSTablePositron = 0;
00370 }
00371
00372 if (energyGrid)
00373 delete energyGrid;
00374
00375 if (fPenelopeFSHelper)
00376 fPenelopeFSHelper->ClearTables();
00377
00378 if (verboseLevel > 2)
00379 G4cout << "G4PenelopeBremsstrahlungModel: cleared tables" << G4endl;
00380
00381 return;
00382 }
00383
00384
00385
00386 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
00387 const G4MaterialCutsCouple*)
00388 {
00389 return fIntrinsicLowEnergyLimit;
00390 }
00391
00392
00393
00394 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
00395 {
00396
00397
00398
00399
00400
00401 if (verboseLevel > 2)
00402 {
00403 G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
00404 G4cout << "for e+/e- in " << mat->GetName() << G4endl;
00405 }
00406
00407
00408 if (energyGrid->GetVectorLength() != nBins)
00409 {
00410 G4ExceptionDescription ed;
00411 ed << "Energy Grid looks not initialized" << G4endl;
00412 ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
00413 G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
00414 "em2016",FatalException,ed);
00415 }
00416
00417 G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
00418 G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
00419
00420 G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
00421
00422
00423
00424 for (size_t bin=0;bin<nBins;bin++)
00425 {
00426 G4double energy = energyGrid->GetLowEdgeEnergy(bin);
00427 G4double XH0=0, XH1=0, XH2=0;
00428 G4double XS0=0, XS1=0, XS2=0;
00429
00430
00431 G4double fact = fPenelopeFSHelper->GetEffectiveZSquared(mat)*
00432 ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
00433 (energy*(energy+2.0*electron_mass_c2)));
00434
00435 G4double restrictedCut = cut/energy;
00436
00437
00438
00439 size_t nBinsX = fPenelopeFSHelper->GetNBinsX();
00440 G4double* tempData = new G4double[nBinsX];
00441 G4double logene = std::log(energy);
00442 for (size_t ix=0;ix<nBinsX;ix++)
00443 {
00444
00445 G4double val = (*table)[ix]->Value(logene);
00446 tempData[ix] = std::exp(val);
00447 }
00448
00449 G4double XH0A = 0.;
00450 if (restrictedCut <= 1)
00451 XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
00452 fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
00453 G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
00454 restrictedCut,0);
00455 G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
00456 restrictedCut,1);
00457 G4double XH1A=0, XH2A=0;
00458 if (restrictedCut <=1)
00459 {
00460 XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
00461 XS1A;
00462 XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
00463 XS2A;
00464 }
00465 delete[] tempData;
00466
00467 XH0 = XH0A*fact;
00468 XS1 = XS1A*fact*energy;
00469 XH1 = XH1A*fact*energy;
00470 XS2 = XS2A*fact*energy*energy;
00471 XH2 = XH2A*fact*energy*energy;
00472
00473 XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
00474
00475
00476 G4double posCorrection = GetPositronXSCorrection(mat,energy);
00477 XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
00478 XH1*posCorrection,
00479 XH2*posCorrection,
00480 XS0,
00481 XS1*posCorrection,
00482 XS2*posCorrection);
00483 }
00484
00485
00486 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00487 XSTableElectron->insert(std::make_pair(theKey,XSEntry));
00488 XSTablePositron->insert(std::make_pair(theKey,XSEntry2));
00489
00490 return;
00491 }
00492
00493
00494
00495
00496 G4PenelopeCrossSection*
00497 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
00498 const G4Material* mat,
00499 G4double cut)
00500 {
00501 if (part != G4Electron::Electron() && part != G4Positron::Positron())
00502 {
00503 G4ExceptionDescription ed;
00504 ed << "Invalid particle: " << part->GetParticleName() << G4endl;
00505 G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00506 "em0001",FatalException,ed);
00507 return NULL;
00508 }
00509
00510 if (part == G4Electron::Electron())
00511 {
00512 if (!XSTableElectron)
00513 {
00514 G4String excep = "The Cross Section Table for e- was not initialized correctly!";
00515 G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00516 "em2013",FatalException,excep);
00517 return NULL;
00518 }
00519 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00520 if (XSTableElectron->count(theKey))
00521 return XSTableElectron->find(theKey)->second;
00522 else
00523 {
00524 BuildXSTable(mat,cut);
00525 if (XSTableElectron->count(theKey))
00526 return XSTableElectron->find(theKey)->second;
00527 else
00528 {
00529 G4ExceptionDescription ed;
00530 ed << "Unable to build e- table for " << mat->GetName() << G4endl;
00531 G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00532 "em2009",FatalException,ed);
00533 }
00534 }
00535 }
00536 if (part == G4Positron::Positron())
00537 {
00538 if (!XSTablePositron)
00539 {
00540 G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
00541 G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00542 "em2013",FatalException,excep);
00543 return NULL;
00544 }
00545 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
00546 if (XSTablePositron->count(theKey))
00547 return XSTablePositron->find(theKey)->second;
00548 else
00549 {
00550 BuildXSTable(mat,cut);
00551 if (XSTablePositron->count(theKey))
00552 return XSTablePositron->find(theKey)->second;
00553 else
00554 {
00555 G4ExceptionDescription ed;
00556 ed << "Unable to build e+ table for " << mat->GetName() << G4endl;
00557 G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
00558 "em2009",FatalException,ed);
00559 }
00560 }
00561 }
00562 return NULL;
00563 }
00564
00565
00566
00567
00568
00569 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
00570 G4double energy)
00571 {
00572
00573
00574
00575
00576
00577 G4double t=std::log(1.0+1e6*energy/
00578 (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
00579 G4double corr = 1.0-std::exp(-t*(1.2359e-1-t*(6.1274e-2-t*
00580 (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
00581 (7.0568e-5-t*
00582 1.8080e-6)))))));
00583 return corr;
00584 }