G4CrossSectionDataStore.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4CrossSectionDataStore
00034 //
00035 // Modifications:
00036 // 23.01.2009 V.Ivanchenko add destruction of data sets
00037 // 29.04.2010 G.Folger     modifictaions for integer A & Z
00038 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTable
00039 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D.Wright redesign the class
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;  //ALB 14-Aug-2012 Coverity fix.
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     // element wise cross section
00117     elmCrossSection = dataSetList[i]->GetElementCrossSection(part, Z, mat);
00118 
00119   } else {
00120     // isotope wise cross section
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       //ALB 14-Aug-2012 Coverity fix.  return elmCrossSection;
00130     }
00131     // user-defined isotope abundances        
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   //G4cout << "xsec(barn)= " <<  xsec/barn << G4endl;
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   // this methods is called after the check that dataSetList[idx] 
00156   // depend on isotopes, so for this DataSet only isotopes are checked
00157 
00158   // isotope-wise cross section does exist
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     // seach for other dataSet
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   //return dataSetList[idx]->ComputeCrossSection(part, elm, mat);
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   // zero cross section case
00223   if(0.0 >= cross) { return anElement; }
00224 
00225   // select element from a compound 
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     // element-wise cross section
00244     // isotope cross section is not computed
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     // isotope abundances        
00255     G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
00256     iso = (*isoVector)[0];
00257 
00258     // more than 1 isotope
00259     if(1 < nIso) { 
00260       iso = dataSetList[i]->SelectIsotope(anElement, part->GetKineticEnergy());
00261     }
00262 
00263   } else {
00264 
00265     //----------------------------------------------------------------
00266     // isotope-wise cross section
00267     // isotope cross section is computed
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     // user-defined isotope abundances        
00281     G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
00282     iso = (*isoVector)[0];
00283 
00284     // more than 1 isotope
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   // Print out all cross section data sets used and the energies at
00330   // which they apply
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   // Write cross section data set info to html physics list
00357   // documentation page
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 }

Generated on Mon May 27 17:47:58 2013 for Geant4 by  doxygen 1.4.7