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 #include "G4LivermoreIonisationModel.hh"
00055 #include "G4PhysicalConstants.hh"
00056 #include "G4SystemOfUnits.hh"
00057 #include "G4ParticleDefinition.hh"
00058 #include "G4MaterialCutsCouple.hh"
00059 #include "G4ProductionCutsTable.hh"
00060 #include "G4DynamicParticle.hh"
00061 #include "G4Element.hh"
00062 #include "G4ParticleChangeForLoss.hh"
00063 #include "G4Electron.hh"
00064 #include "G4CrossSectionHandler.hh"
00065 #include "G4VEMDataSet.hh"
00066 #include "G4eIonisationCrossSectionHandler.hh"
00067 #include "G4eIonisationSpectrum.hh"
00068 #include "G4VEnergySpectrum.hh"
00069 #include "G4SemiLogInterpolation.hh"
00070 #include "G4AtomicTransitionManager.hh"
00071 #include "G4AtomicShell.hh"
00072
00073
00074
00075
00076 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
00077 const G4String& nam) :
00078 G4VEmModel(nam), fParticleChange(0),
00079 isInitialised(false),
00080 crossSectionHandler(0), energySpectrum(0)
00081 {
00082 fIntrinsicLowEnergyLimit = 10.0*eV;
00083 fIntrinsicHighEnergyLimit = 100.0*GeV;
00084
00085 verboseLevel = 0;
00086
00087 transitionManager = G4AtomicTransitionManager::Instance();
00088 }
00089
00090
00091
00092 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
00093 {
00094 delete energySpectrum;
00095 delete crossSectionHandler;
00096 }
00097
00098
00099
00100 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
00101 const G4DataVector& cuts)
00102 {
00103
00104 if (particle != G4Electron::Electron())
00105 {
00106 G4Exception("G4LivermoreIonisationModel::Initialise",
00107 "em0002",FatalException,"Livermore Ionisation Model is applicable only to electrons");
00108 }
00109
00110
00111 if (energySpectrum)
00112 {
00113 delete energySpectrum;
00114 energySpectrum = 0;
00115 }
00116 energySpectrum = new G4eIonisationSpectrum();
00117 if (verboseLevel > 3)
00118 G4cout << "G4VEnergySpectrum is initialized" << G4endl;
00119
00120
00121 if (crossSectionHandler)
00122 {
00123 delete crossSectionHandler;
00124 crossSectionHandler = 0;
00125 }
00126
00127 const size_t nbins = 20;
00128 G4double emin = LowEnergyLimit();
00129 G4double emax = HighEnergyLimit();
00130 G4int ndec = G4int(std::log10(emax/emin) + 0.5);
00131 if(ndec <= 0) { ndec = 1; }
00132
00133 G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
00134 crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
00135 emin,emax,nbins*ndec);
00136 crossSectionHandler->Clear();
00137 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
00138
00139 G4VEMDataSet* emdata =
00140 crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
00141
00142
00143 delete emdata;
00144
00145 if (verboseLevel > 0)
00146 {
00147 G4cout << "Livermore Ionisation model is initialized " << G4endl
00148 << "Energy range: "
00149 << LowEnergyLimit() / keV << " keV - "
00150 << HighEnergyLimit() / GeV << " GeV"
00151 << G4endl;
00152 }
00153
00154 if (verboseLevel > 3)
00155 {
00156 G4cout << "Cross section data: " << G4endl;
00157 crossSectionHandler->PrintData();
00158 G4cout << "Parameters: " << G4endl;
00159 energySpectrum->PrintData();
00160 }
00161
00162 if(isInitialised) { return; }
00163 fParticleChange = GetParticleChangeForLoss();
00164 isInitialised = true;
00165 }
00166
00167
00168
00169 G4double
00170 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00171 G4double energy,
00172 G4double Z, G4double,
00173 G4double cutEnergy,
00174 G4double)
00175 {
00176 G4int iZ = (G4int) Z;
00177 if (!crossSectionHandler)
00178 {
00179 G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
00180 "em1007",FatalException,"The cross section handler is not correctly initialized");
00181 return 0;
00182 }
00183
00184
00185 G4double cs =
00186 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
00187
00188 if (verboseLevel > 1)
00189 {
00190 G4cout << "G4LivermoreIonisationModel " << G4endl;
00191 G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
00192 energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
00193 }
00194 return cs;
00195 }
00196
00197
00198
00199
00200 G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
00201 const G4ParticleDefinition* ,
00202 G4double kineticEnergy,
00203 G4double cutEnergy)
00204 {
00205 G4double sPower = 0.0;
00206
00207 const G4ElementVector* theElementVector = material->GetElementVector();
00208 size_t NumberOfElements = material->GetNumberOfElements() ;
00209 const G4double* theAtomicNumDensityVector =
00210 material->GetAtomicNumDensityVector();
00211
00212
00213 for (size_t iel=0; iel<NumberOfElements; iel++ )
00214 {
00215 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
00216 G4int nShells = transitionManager->NumberOfShells(iZ);
00217 for (G4int n=0; n<nShells; n++)
00218 {
00219 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
00220 kineticEnergy, n);
00221 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
00222 sPower += e * cs * theAtomicNumDensityVector[iel];
00223 }
00224 G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
00225 sPower += esp * theAtomicNumDensityVector[iel];
00226 }
00227
00228 if (verboseLevel > 2)
00229 {
00230 G4cout << "G4LivermoreIonisationModel " << G4endl;
00231 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
00232 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
00233 }
00234
00235 return sPower;
00236 }
00237
00238
00239
00240 void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00241 const G4MaterialCutsCouple* couple,
00242 const G4DynamicParticle* aDynamicParticle,
00243 G4double cutE,
00244 G4double maxE)
00245 {
00246
00247 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
00248
00249 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
00250 {
00251 fParticleChange->SetProposedKineticEnergy(0.);
00252 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
00253 return;
00254 }
00255
00256
00257 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
00258 G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
00259 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
00260 G4double bindingEnergy = shell->BindingEnergy();
00261
00262
00263 G4double energyMax =
00264 std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
00265 G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
00266 kineticEnergy, shellIndex);
00267
00268 if (energyDelta == 0.)
00269 { return; }
00270
00271
00272 G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
00273 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
00274
00275
00276 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
00277 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
00278
00279 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
00280 / (deltaMom * primaryMom);
00281 if (cost > 1.) { cost = 1.; }
00282 G4double sint = std::sqrt((1. - cost)*(1. + cost));
00283 G4double phi = twopi * G4UniformRand();
00284 G4double dirx = sint * std::cos(phi);
00285 G4double diry = sint * std::sin(phi);
00286 G4double dirz = cost;
00287
00288
00289 G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
00290 G4ThreeVector deltaDir(dirx,diry,dirz);
00291 deltaDir.rotateUz(primaryDirection);
00292
00293 dirx = deltaDir.x();
00294 diry = deltaDir.y();
00295 dirz = deltaDir.z();
00296
00297
00298
00299 cost = 2.0*G4UniformRand() - 1.0;
00300 sint = std::sqrt(1. - cost*cost);
00301 phi = twopi * G4UniformRand();
00302 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
00303 / deltaMom;
00304 dirx += del* sint * std::cos(phi);
00305 diry += del* sint * std::sin(phi);
00306 dirz += del* cost;
00307
00308
00309 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
00310 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
00311 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
00312
00313
00314 G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
00315 theDeltaRay->SetKineticEnergy(energyDelta);
00316 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
00317 dirx *= norm;
00318 diry *= norm;
00319 dirz *= norm;
00320 theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
00321 theDeltaRay->SetDefinition(G4Electron::Electron());
00322 fvect->push_back(theDeltaRay);
00323
00324
00325 G4double theEnergyDeposit = bindingEnergy;
00326
00327
00328
00329 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
00330 if(finalKinEnergy < 0.0)
00331 {
00332 theEnergyDeposit += finalKinEnergy;
00333 finalKinEnergy = 0.0;
00334 }
00335 else
00336 {
00337 G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
00338 finalPx *= normLocal;
00339 finalPy *= normLocal;
00340 finalPz *= normLocal;
00341 fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
00342 }
00343 fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
00344
00345 if (theEnergyDeposit < 0)
00346 {
00347 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
00348 << theEnergyDeposit/eV << " eV" << G4endl;
00349 theEnergyDeposit = 0.0;
00350 }
00351
00352
00353 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
00354
00355 if (verboseLevel > 1)
00356 {
00357 G4cout << "-----------------------------------------------------------" << G4endl;
00358 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
00359 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
00360 G4cout << "-----------------------------------------------------------" << G4endl;
00361 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
00362 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
00363 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
00364 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
00365 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
00366 << " keV" << G4endl;
00367 G4cout << "-----------------------------------------------------------" << G4endl;
00368 }
00369 return;
00370 }
00371
00372
00373