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
00037
00038 #include "G4PenelopeOscillatorManager.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4PenelopeIonisationCrossSection.hh"
00041 #include "G4PenelopeIonisationXSHandler.hh"
00042 #include "G4PenelopeCrossSection.hh"
00043 #include "G4Electron.hh"
00044 #include "G4AtomicTransitionManager.hh"
00045
00046 G4PenelopeIonisationCrossSection::G4PenelopeIonisationCrossSection() :
00047 G4VhShellCrossSection("Penelope"),shellIDTable(0),
00048 theCrossSectionHandler(0)
00049 {
00050 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
00051 nMaxLevels = 9;
00052
00053
00054
00055
00056
00057 verboseLevel = 0;
00058
00059 fLowEnergyLimit = 10.0*eV;
00060 fHighEnergyLimit = 100.0*GeV;
00061
00062 transitionManager = G4AtomicTransitionManager::Instance();
00063 }
00064
00065
00066
00067 G4PenelopeIonisationCrossSection::~G4PenelopeIonisationCrossSection()
00068 {
00069 if (theCrossSectionHandler)
00070 delete theCrossSectionHandler;
00071 }
00072
00073
00074
00075 G4double G4PenelopeIonisationCrossSection::CrossSection(G4int Z,
00076 G4AtomicShellEnumerator shell,
00077 G4double incidentEnergy,
00078 G4double ,
00079 const G4Material* material)
00080 {
00081 if (verboseLevel > 1)
00082 G4cout << "Entering in method G4PenelopeIonisationCrossSection::CrossSection()" << G4endl;
00083
00084 G4double cross = 0.;
00085
00086
00087 if (!material)
00088 {
00089
00090 G4ExceptionDescription ed;
00091 ed << "The method has been called with a NULL G4Material pointer" << G4endl;
00092 G4Exception("G4PenelopeIonisationCrossSection::CrossSection()","em2042",
00093 FatalException,ed);
00094 return cross;
00095 }
00096
00097 if (!theCrossSectionHandler)
00098 theCrossSectionHandler = new G4PenelopeIonisationXSHandler();
00099
00100 G4int nmax = std::min(nMaxLevels,transitionManager->NumberOfShells(Z));
00101
00102 if(G4int(shell) < nmax &&
00103 incidentEnergy >= fLowEnergyLimit && incidentEnergy <= fHighEnergyLimit)
00104 {
00105
00106
00107
00108
00109 G4int index = FindShellIDIndex(material,Z,shell);
00110
00111
00112 if (index < 0)
00113 return cross;
00114
00115 G4PenelopeCrossSection* theXS =
00116 theCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),
00117 material,
00118 0.);
00119
00120
00121 G4PenelopeOscillator* theOsc = oscManager->GetOscillatorIonisation(material,index);
00122 if (theOsc->GetParentZ() != Z || theOsc->GetShellFlag()-1 != G4int(shell))
00123 {
00124
00125 G4ExceptionDescription ed;
00126 ed << "There is something wrong here: it looks like the index is wrong" << G4endl;
00127 ed << "Requested: shell " << G4int(shell) << " and Z = " << Z << G4endl;
00128 ed << "Retrieved: " << theOsc->GetShellFlag()-1 << " and Z = " << theOsc->GetParentZ() << G4endl;
00129 G4Exception("G4PenelopeIonisationCrossSection::CrossSection()","em2043",
00130 JustWarning,ed);
00131 return cross;
00132 }
00133
00134
00135 G4double crossPerMolecule = (theXS) ? theXS->GetShellCrossSection(index,incidentEnergy) : 0.;
00136
00137
00138
00139 G4double atomsPerMolec = oscManager->GetNumberOfZAtomsPerMolecule(material,Z);
00140 if (atomsPerMolec)
00141 cross = crossPerMolecule/atomsPerMolec;
00142
00143
00144 if (verboseLevel > 0)
00145 {
00146 G4cout << "Cross section of shell " << G4int(shell) << " and Z= " << Z;
00147 G4cout << " of material: " << material->GetName() << " and energy = " << incidentEnergy/keV << " keV" << G4endl;
00148 G4cout << "--> " << cross/barn << " barn" << G4endl;
00149 G4cout << "Shell binding energy: " << theOsc->GetIonisationEnergy()/eV << " eV;" ;
00150 G4cout << " resonance energy: " << theOsc->GetResonanceEnergy()/eV << "eV" << G4endl;
00151 if (verboseLevel > 2)
00152 {
00153 G4cout << "Cross section per molecule: " << crossPerMolecule/barn << " barn" << G4endl;
00154 G4cout << "Atoms " << Z << " per molecule: " << atomsPerMolec << G4endl;
00155 }
00156 }
00157 }
00158
00159 return cross;
00160 }
00161
00162
00163 std::vector<G4double>
00164 G4PenelopeIonisationCrossSection::GetCrossSection(G4int Z,
00165 G4double kinEnergy,
00166 G4double, G4double,
00167 const G4Material* mat)
00168 {
00169 G4int nmax = std::min(nMaxLevels,transitionManager->NumberOfShells(Z));
00170 std::vector<G4double> vec(nmax,0.0);
00171 for(G4int i=0; i<nmax; ++i) {
00172 vec[i] = CrossSection(Z, G4AtomicShellEnumerator(i), kinEnergy,0.,mat);
00173 }
00174 return vec;
00175 }
00176
00177
00178
00179 std::vector<G4double>
00180 G4PenelopeIonisationCrossSection::Probabilities(G4int Z,
00181 G4double kinEnergy,
00182 G4double,
00183 G4double,
00184 const G4Material* mat)
00185 {
00186 std::vector<G4double> vec = GetCrossSection(Z, kinEnergy,0,0,mat);
00187 size_t n = vec.size();
00188 size_t i=0;
00189 G4double sum = 0.0;
00190 for(i=0; i<n; ++i) { sum += vec[i]; }
00191 if(sum > 0.0) {
00192 sum = 1.0/sum;
00193 for(i=0; i<n; ++i) { vec[i] = vec[i]*sum; }
00194 }
00195 return vec;
00196 }
00197
00198
00199
00200 G4int G4PenelopeIonisationCrossSection::FindShellIDIndex(const G4Material* mat,
00201 G4int Z,
00202 G4AtomicShellEnumerator shell)
00203 {
00204 if (verboseLevel > 1)
00205 G4cout << "Entering in method G4PenelopeIonisationCrossSection::FindShellIDIndex()" << G4endl;
00206
00207 if (!shellIDTable)
00208 shellIDTable = new std::map< std::pair<const G4Material*,G4int>, G4DataVector*>;
00209
00210 std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
00211 G4int result = -1;
00212 G4int ishell = G4int(shell);
00213
00214 if (shellIDTable->count(theKey))
00215 {
00216 if (verboseLevel > 2)
00217 G4cout << "FindShellIDIndex: Table already built for " << mat->GetName() << G4endl;
00218 G4DataVector* theVec = shellIDTable->find(theKey)->second;
00219
00220 if (ishell>=0 && ishell < (G4int) theVec->size())
00221 result = (G4int) (*theVec)[ishell];
00222 else
00223 {
00224 G4ExceptionDescription ed;
00225 ed << "Shell ID: " << ishell << " not available for material " << mat->GetName() << " and Z = " <<
00226 Z << G4endl;
00227 G4Exception("G4PenelopeIonisationCrossSection::FindShellIDIndex()","em2041",JustWarning,
00228 ed);
00229 return -1;
00230 }
00231 }
00232 else
00233 {
00234 if (verboseLevel > 2)
00235 G4cout << "FindShellIDIndex: Table to be built for " << mat->GetName() << G4endl;
00236
00237 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
00238 size_t numberOfOscillators = theTable->size();
00239
00240
00241
00242 G4DataVector* dat = new G4DataVector(nMaxLevels,-1);
00243 for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
00244 {
00245 G4PenelopeOscillator* theOsc = (*theTable)[iosc];
00246
00247 if (theOsc->GetParentZ() == Z)
00248 {
00249
00250 G4int shFlag = theOsc->GetShellFlag();
00251
00252 if (shFlag < 30)
00253 (*dat)[shFlag-1] = (G4double) iosc;
00254 if ((shFlag-1) == ishell)
00255 result = (G4int) iosc;
00256 }
00257 }
00258 shellIDTable->insert(std::make_pair(theKey,dat));
00259 }
00260
00261 if (verboseLevel > 1)
00262 G4cout << "Leaving method G4PenelopeIonisationCrossSection::FindShellIDIndex() with index = " << result << G4endl;
00263
00264 return result;
00265
00266 }