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 // Class Description 26 // Class Description 27 // Cross-section data set for a high precision 27 // Cross-section data set for a high precision (based on JENDL_HE evaluated data 28 // libraries) description of elastic scatterin 28 // libraries) description of elastic scattering 20 MeV ~ 3 GeV; 29 // Class Description - End 29 // Class Description - End 30 30 31 // 15-Nov-06 First Implementation is done by T 31 // 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS) 32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 33 // 33 // 34 #include "G4ParticleHPJENDLHEData.hh" 34 #include "G4ParticleHPJENDLHEData.hh" 35 35 36 #include "G4ElementTable.hh" 36 #include "G4ElementTable.hh" 37 #include "G4ParticleHPData.hh" 37 #include "G4ParticleHPData.hh" 38 #include "G4PhysicsFreeVector.hh" 38 #include "G4PhysicsFreeVector.hh" 39 #include "G4Pow.hh" 39 #include "G4Pow.hh" 40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 41 41 42 G4bool G4ParticleHPJENDLHEData::IsApplicable(c 42 G4bool G4ParticleHPJENDLHEData::IsApplicable(const G4DynamicParticle* aP, const G4Element* anE) 43 { 43 { 44 G4bool result = true; 44 G4bool result = true; 45 G4double eKin = aP->GetKineticEnergy(); 45 G4double eKin = aP->GetKineticEnergy(); 46 // if(eKin>20*MeV||aP->GetDefinition()!=G4Ne 46 // if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false; 47 if (eKin < 20 * MeV || 3 * GeV < eKin || aP- 47 if (eKin < 20 * MeV || 3 * GeV < eKin || aP->GetDefinition() != G4Neutron::Neutron()) { 48 result = false; 48 result = false; 49 } 49 } 50 // Element Check 50 // Element Check 51 else if (!(vElement[anE->GetIndex()])) 51 else if (!(vElement[anE->GetIndex()])) 52 result = false; 52 result = false; 53 53 54 return result; 54 return result; 55 } 55 } 56 56 57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa 57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData() 58 { 58 { 59 for (auto& itZ : mIsotope) { 59 for (auto& itZ : mIsotope) { 60 std::map<G4int, G4PhysicsVector*>* pointer 60 std::map<G4int, G4PhysicsVector*>* pointer_map = itZ.second; 61 if (pointer_map != nullptr) { 61 if (pointer_map != nullptr) { 62 for (auto& itA : *pointer_map) { 62 for (auto& itA : *pointer_map) { 63 G4PhysicsVector* pointerPhysicsVector 63 G4PhysicsVector* pointerPhysicsVector = itA.second; 64 if (pointerPhysicsVector != nullptr) { 64 if (pointerPhysicsVector != nullptr) { 65 delete pointerPhysicsVector; 65 delete pointerPhysicsVector; 66 itA.second = NULL; 66 itA.second = NULL; 67 } 67 } 68 } 68 } 69 delete pointer_map; 69 delete pointer_map; 70 itZ.second = NULL; 70 itZ.second = NULL; 71 } 71 } 72 } 72 } 73 mIsotope.clear(); 73 mIsotope.clear(); 74 } 74 } 75 75 76 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa << 76 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData(G4String reaction, G4ParticleDefinition* pd) 77 : G4VCrossSectionDataSet("JENDLHE" + reactio 77 : G4VCrossSectionDataSet("JENDLHE" + reaction + "CrossSection") 78 { 78 { 79 reactionName = reaction; 79 reactionName = reaction; 80 BuildPhysicsTable(*pd); 80 BuildPhysicsTable(*pd); 81 } 81 } 82 82 83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHED 83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHEData() = default; 84 84 85 void G4ParticleHPJENDLHEData::BuildPhysicsTabl 85 void G4ParticleHPJENDLHEData::BuildPhysicsTable(const G4ParticleDefinition& aP) 86 { 86 { 87 particleName = aP.GetParticleName(); 87 particleName = aP.GetParticleName(); 88 88 89 const G4String& baseName = G4FindDataDir("G4 << 89 G4String baseName = G4FindDataDir("G4NEUTRONHPDATA"); 90 const G4String& dirName = baseName + "/JENDL << 90 G4String dirName = baseName + "/JENDL_HE/" + particleName + "/" + reactionName; 91 const G4String& aFSType = "/CrossSection/"; << 91 G4String aFSType = "/CrossSection/"; 92 G4ParticleHPNames theNames; 92 G4ParticleHPNames theNames; 93 93 94 G4String filename; 94 G4String filename; 95 95 96 // Create JENDL_HE data 96 // Create JENDL_HE data 97 // Create map element or isotope 97 // Create map element or isotope 98 98 99 std::size_t numberOfElements = G4Element::Ge 99 std::size_t numberOfElements = G4Element::GetNumberOfElements(); 100 100 101 // make a PhysicsVector for each element 101 // make a PhysicsVector for each element 102 102 103 auto theElementTable = G4Element::GetElement 103 auto theElementTable = G4Element::GetElementTable(); 104 vElement.clear(); 104 vElement.clear(); 105 vElement.resize(numberOfElements); 105 vElement.resize(numberOfElements); 106 for (std::size_t i = 0; i < numberOfElements 106 for (std::size_t i = 0; i < numberOfElements; ++i) { 107 G4Element* theElement = (*theElementTable) 107 G4Element* theElement = (*theElementTable)[i]; 108 vElement[i] = false; 108 vElement[i] = false; 109 109 110 // isotope 110 // isotope 111 auto nIso = (G4int)(*theElementTable)[i]-> 111 auto nIso = (G4int)(*theElementTable)[i]->GetNumberOfIsotopes(); 112 auto Z = (G4int)(*theElementTable)[i]->Get 112 auto Z = (G4int)(*theElementTable)[i]->GetZ(); 113 for (G4int i1 = 0; i1 < nIso; ++i1) { 113 for (G4int i1 = 0; i1 < nIso; ++i1) { 114 G4int A = theElement->GetIsotope(i1)->Ge 114 G4int A = theElement->GetIsotope(i1)->GetN(); 115 115 116 if (isThisNewIsotope(Z, A)) { 116 if (isThisNewIsotope(Z, A)) { 117 std::stringstream ss; 117 std::stringstream ss; 118 ss << dirName << aFSType << Z << "_" << A << 118 ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName(Z - 1); 119 filename = ss.str(); 119 filename = ss.str(); 120 std::fstream file; 120 std::fstream file; 121 file.open(filename, std::fstream::in); 121 file.open(filename, std::fstream::in); 122 G4int dummy; 122 G4int dummy; 123 file >> dummy; 123 file >> dummy; 124 if (file.good()) { 124 if (file.good()) { 125 vElement[i] = true; 125 vElement[i] = true; 126 126 127 // read the file 127 // read the file 128 G4PhysicsVector* aPhysVec = readAFile(&fil 128 G4PhysicsVector* aPhysVec = readAFile(&file); 129 registAPhysicsVector(Z, A, aPhysVec); 129 registAPhysicsVector(Z, A, aPhysVec); 130 } 130 } 131 file.close(); 131 file.close(); 132 } 132 } 133 } 133 } 134 } 134 } 135 } 135 } 136 136 137 void G4ParticleHPJENDLHEData::DumpPhysicsTable 137 void G4ParticleHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition&) 138 {} 138 {} 139 139 140 G4double G4ParticleHPJENDLHEData::GetCrossSect 140 G4double G4ParticleHPJENDLHEData::GetCrossSection(const G4DynamicParticle* aP, 141 141 const G4Element* anE, G4double) 142 { 142 { 143 // Primary energy >20MeV 143 // Primary energy >20MeV 144 // Thus not taking into account of Doppler b 144 // Thus not taking into account of Doppler broadening 145 // also not taking into account of Target th 145 // also not taking into account of Target thermal motions 146 146 147 G4double result = 0; 147 G4double result = 0; 148 148 149 G4double ek = aP->GetKineticEnergy(); 149 G4double ek = aP->GetKineticEnergy(); 150 150 151 auto nIso = (G4int)anE->GetNumberOfIsotopes( 151 auto nIso = (G4int)anE->GetNumberOfIsotopes(); 152 auto Z = (G4int)anE->GetZ(); 152 auto Z = (G4int)anE->GetZ(); 153 for (G4int i1 = 0; i1 < nIso; ++i1) { 153 for (G4int i1 = 0; i1 < nIso; ++i1) { 154 G4int A = anE->GetIsotope(i1)->GetN(); 154 G4int A = anE->GetIsotope(i1)->GetN(); 155 G4double frac = anE->GetRelativeAbundanceV 155 G4double frac = anE->GetRelativeAbundanceVector()[i1]; 156 // This case does NOT request "*perCent". 156 // This case does NOT request "*perCent". 157 result += frac * getXSfromThisIsotope(Z, A 157 result += frac * getXSfromThisIsotope(Z, A, ek); 158 } 158 } 159 return result; 159 return result; 160 } 160 } 161 161 162 G4PhysicsVector* G4ParticleHPJENDLHEData::read 162 G4PhysicsVector* G4ParticleHPJENDLHEData::readAFile(std::fstream* file) 163 { 163 { 164 G4int dummy; 164 G4int dummy; 165 G4int len; 165 G4int len; 166 *file >> dummy; 166 *file >> dummy; 167 *file >> len; 167 *file >> len; 168 168 169 std::vector<G4double> v_e; 169 std::vector<G4double> v_e; 170 std::vector<G4double> v_xs; 170 std::vector<G4double> v_xs; 171 171 172 for (G4int i = 0; i < len; ++i) { 172 for (G4int i = 0; i < len; ++i) { 173 G4double e; 173 G4double e; 174 G4double xs; 174 G4double xs; 175 175 176 *file >> e; 176 *file >> e; 177 *file >> xs; 177 *file >> xs; 178 // data are written in eV and barn. 178 // data are written in eV and barn. 179 v_e.push_back(e * eV); 179 v_e.push_back(e * eV); 180 v_xs.push_back(xs * barn); 180 v_xs.push_back(xs * barn); 181 } 181 } 182 182 183 auto aPhysVec = new G4PhysicsFreeVector(stat 183 auto aPhysVec = new G4PhysicsFreeVector(static_cast<std::size_t>(len), v_e.front(), v_e.back()); 184 184 185 for (G4int i = 0; i < len; ++i) { 185 for (G4int i = 0; i < len; ++i) { 186 aPhysVec->PutValues(static_cast<std::size_ 186 aPhysVec->PutValues(static_cast<std::size_t>(i), v_e[i], v_xs[i]); 187 } 187 } 188 188 189 return aPhysVec; 189 return aPhysVec; 190 } 190 } 191 191 192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4 192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4int z, G4int a) 193 { 193 { 194 if (mIsotope.find(z) == mIsotope.end()) retu 194 if (mIsotope.find(z) == mIsotope.end()) return false; 195 if (mIsotope.find(z)->second->find(a) == mIs 195 if (mIsotope.find(z)->second->find(a) == mIsotope.find(z)->second->end()) return false; 196 return true; 196 return true; 197 } 197 } 198 198 199 void G4ParticleHPJENDLHEData::registAPhysicsVe 199 void G4ParticleHPJENDLHEData::registAPhysicsVector(G4int Z, G4int A, G4PhysicsVector* aPhysVec) 200 { 200 { 201 std::pair<G4int, G4PhysicsVector*> aPair = s 201 std::pair<G4int, G4PhysicsVector*> aPair = std::pair<G4int, G4PhysicsVector*>(A, aPhysVec); 202 auto itm = mIsotope.find(Z); 202 auto itm = mIsotope.find(Z); 203 if (itm != mIsotope.cend()) { 203 if (itm != mIsotope.cend()) { 204 itm->second->insert(aPair); 204 itm->second->insert(aPair); 205 } 205 } 206 else { 206 else { 207 auto aMap = new std::map<G4int, G4PhysicsV 207 auto aMap = new std::map<G4int, G4PhysicsVector*>; 208 aMap->insert(aPair); 208 aMap->insert(aPair); 209 mIsotope.insert(std::pair<G4int, std::map< 209 mIsotope.insert(std::pair<G4int, std::map<G4int, G4PhysicsVector*>*>(Z, aMap)); 210 } 210 } 211 } 211 } 212 212 213 G4double G4ParticleHPJENDLHEData::getXSfromThi 213 G4double G4ParticleHPJENDLHEData::getXSfromThisIsotope(G4int Z, G4int A, G4double ek) 214 { 214 { 215 G4double aXSection = 0.0; 215 G4double aXSection = 0.0; 216 216 217 G4PhysicsVector* aPhysVec; 217 G4PhysicsVector* aPhysVec; 218 if (mIsotope.find(Z)->second->find(A) != mIs 218 if (mIsotope.find(Z)->second->find(A) != mIsotope.find(Z)->second->end()) { 219 aPhysVec = mIsotope.find(Z)->second->find( 219 aPhysVec = mIsotope.find(Z)->second->find(A)->second; 220 aXSection = aPhysVec->Value(ek); 220 aXSection = aPhysVec->Value(ek); 221 } 221 } 222 else { 222 else { 223 // Select closest one in the same Z 223 // Select closest one in the same Z 224 G4int delta0 = 99; // no mean for 99 224 G4int delta0 = 99; // no mean for 99 225 for (auto it = mIsotope.find(Z)->second->c 225 for (auto it = mIsotope.find(Z)->second->cbegin(); it != mIsotope.find(Z)->second->cend(); ++it) 226 { 226 { 227 G4int delta = std::abs(A - it->first); 227 G4int delta = std::abs(A - it->first); 228 if (delta < delta0) delta0 = delta; 228 if (delta < delta0) delta0 = delta; 229 } 229 } 230 230 231 // Randomize of selection larger or smalle 231 // Randomize of selection larger or smaller than A 232 if (G4UniformRand() < 0.5) delta0 *= -1; 232 if (G4UniformRand() < 0.5) delta0 *= -1; 233 G4int A1 = A + delta0; 233 G4int A1 = A + delta0; 234 if (mIsotope.find(Z)->second->find(A1) != 234 if (mIsotope.find(Z)->second->find(A1) != mIsotope.find(Z)->second->cend()) { 235 aPhysVec = mIsotope.find(Z)->second->fin 235 aPhysVec = mIsotope.find(Z)->second->find(A1)->second; 236 } 236 } 237 else { 237 else { 238 A1 = A - delta0; 238 A1 = A - delta0; 239 aPhysVec = mIsotope.find(Z)->second->fin 239 aPhysVec = mIsotope.find(Z)->second->find(A1)->second; 240 } 240 } 241 241 242 aXSection = aPhysVec->Value(ek); 242 aXSection = aPhysVec->Value(ek); 243 // X^(2/3) factor 243 // X^(2/3) factor 244 aXSection *= G4Pow::GetInstance()->A23(1.0 244 aXSection *= G4Pow::GetInstance()->A23(1.0 * A / A1); 245 } 245 } 246 246 247 return aXSection; 247 return aXSection; 248 } 248 } 249 249