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 #include "G4CrossSectionDataStore.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4HadronicException.hh"
00045 #include "G4HadTmpUtil.hh"
00046 #include "Randomize.hh"
00047 #include "G4Nucleus.hh"
00048
00049 #include "G4DynamicParticle.hh"
00050 #include "G4Isotope.hh"
00051 #include "G4Element.hh"
00052 #include "G4Material.hh"
00053 #include "G4NistManager.hh"
00054 #include <iostream>
00055
00056
00057 G4CrossSectionDataStore::G4CrossSectionDataStore() :
00058 nDataSetList(0), verboseLevel(0)
00059 {
00060 nist = G4NistManager::Instance();
00061 currentMaterial = elmMaterial = 0;
00062 currentElement = 0;
00063 matParticle = elmParticle = 0;
00064 matKinEnergy = elmKinEnergy = matCrossSection = elmCrossSection = 0.0;
00065 }
00066
00067 G4CrossSectionDataStore::~G4CrossSectionDataStore()
00068 {}
00069
00070 G4double
00071 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
00072 const G4Material* mat)
00073 {
00074 if(mat == currentMaterial && part->GetDefinition() == matParticle
00075 && part->GetKineticEnergy() == matKinEnergy)
00076 { return matCrossSection; }
00077
00078 currentMaterial = mat;
00079 matParticle = part->GetDefinition();
00080 matKinEnergy = part->GetKineticEnergy();
00081 matCrossSection = 0;
00082
00083 G4int nElements = mat->GetNumberOfElements();
00084 const G4double* nAtomsPerVolume = mat->GetVecNbOfAtomsPerVolume();
00085
00086 if(G4int(xsecelm.size()) < nElements) { xsecelm.resize(nElements); }
00087
00088 for(G4int i=0; i<nElements; ++i) {
00089 matCrossSection += nAtomsPerVolume[i] *
00090 GetCrossSection(part, (*mat->GetElementVector())[i], mat);
00091 xsecelm[i] = matCrossSection;
00092 }
00093 return matCrossSection;
00094 }
00095
00096 G4double
00097 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
00098 const G4Element* elm,
00099 const G4Material* mat)
00100 {
00101 if(mat == elmMaterial && elm == currentElement &&
00102 part->GetDefinition() == elmParticle &&
00103 part->GetKineticEnergy() == elmKinEnergy)
00104 { return elmCrossSection; }
00105
00106 elmMaterial = mat;
00107 currentElement = elm;
00108 elmParticle = part->GetDefinition();
00109 elmKinEnergy = part->GetKineticEnergy();
00110 elmCrossSection = 0.0;
00111
00112 G4int i = nDataSetList-1;
00113 G4int Z = G4lrint(elm->GetZ());
00114 if (dataSetList[i]->IsElementApplicable(part, Z, mat)) {
00115
00116
00117 elmCrossSection = dataSetList[i]->GetElementCrossSection(part, Z, mat);
00118
00119 } else {
00120
00121 G4int nIso = elm->GetNumberOfIsotopes();
00122 G4Isotope* iso = 0;
00123
00124 if (0 >= nIso) {
00125 G4cout << " Element " << elm->GetName() << " Z= " << Z
00126 << " has no isotopes " << G4endl;
00127 throw G4HadronicException(__FILE__, __LINE__,
00128 " Isotope vector is not defined");
00129
00130 }
00131
00132 G4IsotopeVector* isoVector = elm->GetIsotopeVector();
00133 G4double* abundVector = elm->GetRelativeAbundanceVector();
00134
00135 for (G4int j = 0; j<nIso; ++j) {
00136 if(abundVector[j] > 0.0) {
00137 iso = (*isoVector)[j];
00138 elmCrossSection += abundVector[j]*
00139 GetIsoCrossSection(part, Z, iso->GetN(), iso, elm, mat, i);
00140 }
00141 }
00142 }
00143
00144 return elmCrossSection;
00145 }
00146
00147 G4double
00148 G4CrossSectionDataStore::GetIsoCrossSection(const G4DynamicParticle* part,
00149 G4int Z, G4int A,
00150 const G4Isotope* iso,
00151 const G4Element* elm,
00152 const G4Material* mat,
00153 G4int idx)
00154 {
00155
00156
00157
00158
00159 if(dataSetList[idx]->IsIsoApplicable(part, Z, A, elm, mat) ) {
00160 return dataSetList[idx]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
00161
00162 } else {
00163
00164 for (G4int j = idx-1; j >= 0; --j) {
00165 if (dataSetList[j]->IsElementApplicable(part, Z, mat)) {
00166 return dataSetList[j]->GetElementCrossSection(part, Z, mat);
00167 } else if (dataSetList[j]->IsIsoApplicable(part, Z, A, elm, mat)) {
00168 return dataSetList[j]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
00169 }
00170 }
00171 }
00172 G4cout << "G4CrossSectionDataStore::GetCrossSection ERROR: "
00173 << " no isotope cross section found"
00174 << G4endl;
00175 G4cout << " for " << part->GetDefinition()->GetParticleName()
00176 << " off Element " << elm->GetName()
00177 << " in " << mat->GetName()
00178 << " Z= " << Z << " A= " << A
00179 << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
00180 throw G4HadronicException(__FILE__, __LINE__,
00181 " no applicable data set found for the isotope");
00182 return 0.0;
00183
00184 }
00185
00186 G4double
00187 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* part,
00188 G4int Z, G4int A,
00189 const G4Isotope* iso,
00190 const G4Element* elm,
00191 const G4Material* mat)
00192 {
00193 for (G4int i = nDataSetList-1; i >= 0; --i) {
00194 if (dataSetList[i]->IsIsoApplicable(part, Z, A, elm, mat) ) {
00195 return dataSetList[i]->GetIsoCrossSection(part, Z, A, iso, elm, mat);
00196 }
00197 }
00198 G4cout << "G4CrossSectionDataStore::GetCrossSection ERROR: "
00199 << " no isotope cross section found"
00200 << G4endl;
00201 G4cout << " for " << part->GetDefinition()->GetParticleName()
00202 << " off Element " << elm->GetName()
00203 << " in " << mat->GetName()
00204 << " Z= " << Z << " A= " << A
00205 << " E(MeV)= " << part->GetKineticEnergy()/MeV << G4endl;
00206 throw G4HadronicException(__FILE__, __LINE__,
00207 " no applicable data set found for the isotope");
00208 return 0.0;
00209 }
00210
00211 G4Element*
00212 G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* part,
00213 const G4Material* mat,
00214 G4Nucleus& target)
00215 {
00216 G4int nElements = mat->GetNumberOfElements();
00217 const G4ElementVector* theElementVector = mat->GetElementVector();
00218 G4Element* anElement = (*theElementVector)[0];
00219
00220 G4double cross = GetCrossSection(part, mat);
00221
00222
00223 if(0.0 >= cross) { return anElement; }
00224
00225
00226 if(1 < nElements) {
00227 cross *= G4UniformRand();
00228 for(G4int i=0; i<nElements; ++i) {
00229 if(cross <= xsecelm[i]) {
00230 anElement = (*theElementVector)[i];
00231 break;
00232 }
00233 }
00234 }
00235
00236 G4int Z = G4lrint(anElement->GetZ());
00237 G4Isotope* iso = 0;
00238
00239 G4int i = nDataSetList-1;
00240 if (dataSetList[i]->IsElementApplicable(part, Z, mat)) {
00241
00242
00243
00244
00245
00246 G4int nIso = anElement->GetNumberOfIsotopes();
00247 if (0 >= nIso) {
00248 G4cout << " Element " << anElement->GetName() << " Z= " << Z
00249 << " has no isotopes " << G4endl;
00250 throw G4HadronicException(__FILE__, __LINE__,
00251 " Isotope vector is not defined");
00252 return anElement;
00253 }
00254
00255 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
00256 iso = (*isoVector)[0];
00257
00258
00259 if(1 < nIso) {
00260 iso = dataSetList[i]->SelectIsotope(anElement, part->GetKineticEnergy());
00261 }
00262
00263 } else {
00264
00265
00266
00267
00268
00269 G4int nIso = anElement->GetNumberOfIsotopes();
00270 cross = 0.0;
00271
00272 if (0 >= nIso) {
00273 G4cout << " Element " << anElement->GetName() << " Z= " << Z
00274 << " has no isotopes " << G4endl;
00275 throw G4HadronicException(__FILE__, __LINE__,
00276 " Isotope vector is not defined");
00277 return anElement;
00278 }
00279
00280
00281 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
00282 iso = (*isoVector)[0];
00283
00284
00285 if(1 < nIso) {
00286 G4double* abundVector = anElement->GetRelativeAbundanceVector();
00287 if(G4int(xseciso.size()) < nIso) { xseciso.resize(nIso); }
00288
00289 for (G4int j = 0; j<nIso; ++j) {
00290 G4double xsec = 0.0;
00291 if(abundVector[j] > 0.0) {
00292 iso = (*isoVector)[j];
00293 xsec = abundVector[j]*
00294 GetIsoCrossSection(part, Z, iso->GetN(), iso, anElement, mat, i);
00295 }
00296 cross += xsec;
00297 xseciso[j] = cross;
00298 }
00299 cross *= G4UniformRand();
00300 for (G4int j = 0; j<nIso; ++j) {
00301 if(cross <= xseciso[j]) {
00302 iso = (*isoVector)[j];
00303 break;
00304 }
00305 }
00306 }
00307 }
00308 target.SetIsotope(iso);
00309 return anElement;
00310 }
00311
00312 void
00313 G4CrossSectionDataStore::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
00314 {
00315 if (nDataSetList == 0)
00316 {
00317 throw G4HadronicException(__FILE__, __LINE__,
00318 "G4CrossSectionDataStore: no data sets registered");
00319 return;
00320 }
00321 for (G4int i=0; i<nDataSetList; ++i) {
00322 dataSetList[i]->BuildPhysicsTable(aParticleType);
00323 }
00324 }
00325
00326 void
00327 G4CrossSectionDataStore::DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
00328 {
00329
00330
00331
00332 if (nDataSetList == 0) {
00333 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
00334 << " no data sets registered" << G4endl;
00335 return;
00336 }
00337
00338 for (G4int i = nDataSetList-1; i >= 0; --i) {
00339 G4double e1 = dataSetList[i]->GetMinKinEnergy();
00340 G4double e2 = dataSetList[i]->GetMaxKinEnergy();
00341 if (i < nDataSetList-1) { G4cout << " "; }
00342 G4cout << std::setw(25) << dataSetList[i]->GetName() << ": Emin(GeV)= "
00343 << std::setw(4) << e1/GeV << " Emax(GeV)= "
00344 << e2/GeV << G4endl;
00345 if (dataSetList[i]->GetName() == "G4CrossSectionPairGG") {
00346 dataSetList[i]->DumpPhysicsTable(aParticleType);
00347 }
00348 }
00349
00350 G4cout << G4endl;
00351 }
00352
00353 void G4CrossSectionDataStore::DumpHtml(const G4ParticleDefinition&,
00354 std::ofstream& outFile)
00355 {
00356
00357
00358
00359 G4double ehi = 0;
00360 G4double elo = 0;
00361 for (G4int i = nDataSetList-1; i > 0; i--) {
00362 elo = dataSetList[i]->GetMinKinEnergy()/GeV;
00363 ehi = dataSetList[i]->GetMaxKinEnergy()/GeV;
00364 outFile << " <li><b><a href=\"" << dataSetList[i]->GetName() << ".html\"> "
00365 << dataSetList[i]->GetName() << "</a> from "
00366 << elo << " GeV to " << ehi << " GeV </b></li>\n";
00367 }
00368
00369 G4double defaultHi = dataSetList[0]->GetMaxKinEnergy()/GeV;
00370 if (ehi < defaultHi) {
00371 outFile << " <li><b><a href=\"" << dataSetList[0]->GetName() << ".html\"> "
00372 << dataSetList[0]->GetName() << "</a> from "
00373 << ehi << " GeV to " << defaultHi << " GeV </b></li>\n";
00374 }
00375 }