Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4VCrossSectionDataSet.cc 89024 2015-03-18 08:17:25Z gcosmo $ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class file 30 // GEANT4 Class file 29 // 31 // >> 32 // 30 // File name: G4VCrossSectionDataSet 33 // File name: G4VCrossSectionDataSet 31 // 34 // 32 // Author F.W. Jones, TRIUMF, 20-JAN-97 35 // Author F.W. Jones, TRIUMF, 20-JAN-97 33 // 36 // 34 // Modifications: 37 // Modifications: 35 // 23.01.2009 V.Ivanchenko move constructor an 38 // 23.01.2009 V.Ivanchenko move constructor and destructor to source 36 // 12.08.2011 G.Folger, V.Ivanchenko, T.Koi, D 39 // 12.08.2011 G.Folger, V.Ivanchenko, T.Koi, D.Wright redesign the class 37 // 40 // 38 41 39 #include "G4VCrossSectionDataSet.hh" 42 #include "G4VCrossSectionDataSet.hh" 40 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 41 #include "G4CrossSectionDataSetRegistry.hh" 44 #include "G4CrossSectionDataSetRegistry.hh" >> 45 #include "G4DynamicParticle.hh" 42 #include "G4Material.hh" 46 #include "G4Material.hh" 43 #include "G4Element.hh" 47 #include "G4Element.hh" 44 #include "G4Isotope.hh" 48 #include "G4Isotope.hh" 45 #include "G4NistManager.hh" 49 #include "G4NistManager.hh" >> 50 #include "G4HadronicException.hh" >> 51 #include "G4HadTmpUtil.hh" 46 #include "Randomize.hh" 52 #include "Randomize.hh" 47 #include "G4HadronicParameters.hh" << 48 53 49 G4VCrossSectionDataSet::G4VCrossSectionDataSet 54 G4VCrossSectionDataSet::G4VCrossSectionDataSet(const G4String& nam) : 50 name(nam), << 55 verboseLevel(0),minKinEnergy(0.0),maxKinEnergy(100*TeV),name(nam) 51 maxKinEnergy(G4HadronicParameters::Instance( << 52 { 56 { 53 registry = G4CrossSectionDataSetRegistry::In 57 registry = G4CrossSectionDataSetRegistry::Instance(); 54 registry->Register(this); 58 registry->Register(this); 55 } 59 } 56 60 57 G4VCrossSectionDataSet::~G4VCrossSectionDataSe 61 G4VCrossSectionDataSet::~G4VCrossSectionDataSet() 58 { 62 { 59 registry->DeRegister(this); 63 registry->DeRegister(this); 60 } 64 } 61 65 62 G4bool 66 G4bool 63 G4VCrossSectionDataSet::IsElementApplicable(co 67 G4VCrossSectionDataSet::IsElementApplicable(const G4DynamicParticle*, 64 G4int, 68 G4int, 65 const G4Material*) 69 const G4Material*) 66 { 70 { 67 return false; 71 return false; 68 } 72 } 69 73 70 G4bool 74 G4bool 71 G4VCrossSectionDataSet::IsIsoApplicable(const 75 G4VCrossSectionDataSet::IsIsoApplicable(const G4DynamicParticle*, 72 G4int, 76 G4int, G4int, 73 const 77 const G4Element*, 74 const 78 const G4Material*) 75 { 79 { 76 return false; 80 return false; 77 } 81 } 78 82 79 G4double 83 G4double 80 G4VCrossSectionDataSet::ComputeCrossSection(co 84 G4VCrossSectionDataSet::ComputeCrossSection(const G4DynamicParticle* part, 81 co << 85 const G4Element* elm, 82 co << 86 const G4Material* mat) 83 { 87 { 84 G4int Z = elm->GetZasInt(); << 88 G4int Z = G4lrint(elm->GetZ()); 85 89 86 if (IsElementApplicable(part, Z, mat)) { 90 if (IsElementApplicable(part, Z, mat)) { 87 return GetElementCrossSection(part, Z, mat 91 return GetElementCrossSection(part, Z, mat); 88 } 92 } 89 93 90 // isotope-wise cross section making sum ove 94 // isotope-wise cross section making sum over available 91 // isotope cross sections, which may be inco 95 // isotope cross sections, which may be incomplete, so 92 // the result is corrected 96 // the result is corrected 93 std::size_t nIso = elm->GetNumberOfIsotopes( << 97 G4int nIso = elm->GetNumberOfIsotopes(); 94 G4double fact = 0.0; 98 G4double fact = 0.0; 95 G4double xsec = 0.0; 99 G4double xsec = 0.0; >> 100 G4Isotope* iso = 0; 96 101 97 // user-defined isotope abundances << 102 if (0 < nIso) { 98 const G4IsotopeVector* isoVector = elm->GetI << 103 99 const G4double* abundVector = elm->GetRelati << 104 // user-defined isotope abundances 100 for (std::size_t j=0; j<nIso; ++j) { << 105 G4IsotopeVector* isoVector = elm->GetIsotopeVector(); 101 const G4Isotope* iso = (*isoVector)[j]; << 106 G4double* abundVector = elm->GetRelativeAbundanceVector(); 102 G4int A = iso->GetN(); << 107 103 if(abundVector[j] > 0.0 && IsIsoApplicable << 108 for (G4int j = 0; j<nIso; ++j) { 104 fact += abundVector[j]; << 109 iso = (*isoVector)[j]; 105 xsec += abundVector[j]*GetIsoCrossSectio << 110 G4int A = iso->GetN(); >> 111 if(abundVector[j] > 0.0 && IsIsoApplicable(part, Z, A, elm, mat)) { >> 112 fact += abundVector[j]; >> 113 xsec += abundVector[j]*GetIsoCrossSection(part, Z, A, iso, elm, mat); >> 114 } 106 } 115 } 107 } << 108 return (fact > 0.0) ? xsec/fact : 0.0; << 109 } << 110 116 111 G4double << 117 } else { 112 G4VCrossSectionDataSet::ComputeCrossSectionPer << 118 113 G4double kinEnergy, G4 << 119 // natural isotope abundances 114 const G4ParticleDefini << 120 G4NistManager* nist = G4NistManager::Instance(); 115 const G4Element* elm, << 121 G4int n0 = nist->GetNistFirstIsotopeN(Z); 116 { << 122 G4int nn = nist->GetNumberOfNistIsotopes(Z); 117 G4int Z = elm->GetZasInt(); << 123 for (G4int A = n0; A < n0+nn; ++A) { 118 std::size_t nIso = elm->GetNumberOfIsotopes( << 124 G4double abund = nist->GetIsotopeAbundance(Z, A); 119 G4double xsec = 0.0; << 125 if(abund > 0.0 && IsIsoApplicable(part, Z, A, elm, mat)) { 120 const G4IsotopeVector* isoVector = elm->GetI << 126 fact += abund; 121 const G4double* abundVector = elm->GetRelati << 127 xsec += abund*GetIsoCrossSection(part, Z, A, iso, elm, mat); 122 for (std::size_t j=0; j<nIso; ++j) { << 128 } 123 const G4Isotope* iso = (*isoVector)[j]; << 129 } 124 G4int A = iso->GetN(); << 125 xsec += abundVector[j]* << 126 ComputeIsoCrossSection(kinEnergy, loge, << 127 } 130 } >> 131 if(fact > 0.0) { xsec /= fact; } 128 return xsec; 132 return xsec; 129 } 133 } 130 134 131 G4double 135 G4double 132 G4VCrossSectionDataSet::GetElementCrossSection 136 G4VCrossSectionDataSet::GetElementCrossSection(const G4DynamicParticle* dynPart, 133 G4int Z, 137 G4int Z, 134 const G4Material* mat) 138 const G4Material* mat) 135 { 139 { 136 G4ExceptionDescription ed; << 140 G4cout << "G4VCrossSectionDataSet::GetCrossSection per element ERROR: " 137 ed << "GetElementCrossSection is not impleme << 141 << " there is no cross section for " 138 << "Particle: " << dynPart->GetDefinition << 142 << dynPart->GetDefinition()->GetParticleName() 139 << " Ekin(MeV)= " << dynPart->GetKineti << 143 << " E(MeV)= " << dynPart->GetKineticEnergy()/MeV; 140 if(nullptr != mat) { ed << " material: " << << 144 if(mat) { G4cout << " inside " << mat->GetName(); } 141 ed << " target Z= " << Z << G4endl; << 145 G4cout << " for Z= " << Z << G4endl; 142 G4Exception("G4VCrossSectionDataSet::GetElem << 146 throw G4HadronicException(__FILE__, __LINE__, 143 FatalException, ed); << 147 "G4VCrossSectionDataSet::GetElementCrossSection is absent"); 144 return 0.0; 148 return 0.0; 145 } 149 } 146 150 147 G4double 151 G4double 148 G4VCrossSectionDataSet::GetIsoCrossSection(con 152 G4VCrossSectionDataSet::GetIsoCrossSection(const G4DynamicParticle* dynPart, 149 G4int Z, G4int A, 153 G4int Z, G4int A, 150 const G4Isotope*, 154 const G4Isotope*, 151 const G4Element* elm, 155 const G4Element* elm, 152 const G4Material* mat) 156 const G4Material* mat) 153 { 157 { 154 G4ExceptionDescription ed; << 158 G4cout << "G4VCrossSectionDataSet::GetCrossSection per isotope ERROR: " 155 ed << "GetIsoCrossSection is not implemented << 159 << " there is no cross section for " 156 << "Particle: " << dynPart->GetDefinition << 160 << dynPart->GetDefinition()->GetParticleName() 157 << " Ekin(MeV)= " << dynPart->GetKineti << 161 << " E(MeV)= " << dynPart->GetKineticEnergy()/MeV; 158 if(nullptr != mat) { ed << " material: " << << 162 if(mat) { G4cout << " inside " << mat->GetName(); } 159 if(nullptr != elm) { ed << " element: " << e << 163 if(elm) { G4cout << " for " << elm->GetName(); } 160 ed << " target Z= " << Z << " A= " << A << G << 164 G4cout << " Z= " << Z << " A= " << A << G4endl; 161 G4Exception("G4VCrossSectionDataSet::GetIsoC << 165 throw G4HadronicException(__FILE__, __LINE__, 162 FatalException, ed); << 166 "G4VCrossSectionDataSet::GetIsoCrossSection is absent"); 163 return 0.0; << 164 } << 165 << 166 G4double << 167 G4VCrossSectionDataSet::ComputeIsoCrossSection << 168 const G4ParticleDefinition* p << 169 G4int Z, G4int A, << 170 const G4Isotope*, << 171 const G4Element* elm, << 172 const G4Material* mat) << 173 { << 174 G4ExceptionDescription ed; << 175 ed << "GetIsoCrossSection is not implemented << 176 << "Particle: " << pd->GetParticleName() << 177 << " Ekin(MeV)= " << kinEnergy/CLHEP::M << 178 if(nullptr != mat) { ed << " material: " << << 179 if(nullptr != elm) { ed << " element: " << e << 180 ed << " target Z= " << Z << " A= " << A << G << 181 G4Exception("G4VCrossSectionDataSet::GetIsoC << 182 FatalException, ed); << 183 return 0.0; 167 return 0.0; 184 } 168 } 185 169 186 const G4Isotope* << 170 G4Isotope* 187 G4VCrossSectionDataSet::SelectIsotope(const G4 << 171 G4VCrossSectionDataSet::SelectIsotope(const G4Element* anElement, G4double) 188 G4double << 189 { 172 { 190 G4int nIso = (G4int)anElement->GetNumberOfIs << 173 G4int nIso = anElement->GetNumberOfIsotopes(); 191 const G4Isotope* iso = anElement->GetIsotope << 174 G4IsotopeVector* isoVector = anElement->GetIsotopeVector(); >> 175 G4Isotope* iso = (*isoVector)[0]; 192 176 193 // more than 1 isotope 177 // more than 1 isotope 194 if(1 < nIso) { 178 if(1 < nIso) { 195 const G4double* abundVector = anElement->G << 179 G4double* abundVector = anElement->GetRelativeAbundanceVector(); 196 G4double sum = 0.0; 180 G4double sum = 0.0; 197 G4double q = G4UniformRand(); 181 G4double q = G4UniformRand(); 198 for (G4int j=0; j<nIso; ++j) { << 182 for (G4int j = 0; j<nIso; ++j) { 199 sum += abundVector[j]; 183 sum += abundVector[j]; 200 if(q <= sum) { 184 if(q <= sum) { 201 iso = anElement->GetIsotope(j); << 185 iso = (*isoVector)[j]; 202 break; 186 break; 203 } 187 } 204 } 188 } 205 } 189 } 206 return iso; 190 return iso; 207 } 191 } 208 192 209 void G4VCrossSectionDataSet::BuildPhysicsTable 193 void G4VCrossSectionDataSet::BuildPhysicsTable(const G4ParticleDefinition&) 210 {} 194 {} 211 195 212 void G4VCrossSectionDataSet::DumpPhysicsTable( 196 void G4VCrossSectionDataSet::DumpPhysicsTable(const G4ParticleDefinition&) 213 {} 197 {} 214 198 215 void G4VCrossSectionDataSet::CrossSectionDescr 199 void G4VCrossSectionDataSet::CrossSectionDescription(std::ostream& outFile) const 216 { 200 { 217 outFile << "The description for this cross s 201 outFile << "The description for this cross section data set has not been written yet.\n"; 218 } 202 } 219 203