Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // 23 // 27 // ------------------------------------------- << 28 // 24 // 29 // GEANT4 Class file << 25 // GEANT4 physics class: G4CrossSectionDataStore >> 26 // F.W. Jones, TRIUMF, 19-NOV-97 30 // 27 // 31 // << 32 // File name: G4CrossSectionDataStore << 33 // << 34 // Modifications: << 35 // 23.01.2009 V.Ivanchenko add destruction of << 36 // 29.04.2010 G.Folger modifictaions for i << 37 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTa << 38 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D << 39 // 07.03.2013 M.Maire cosmetic in DumpPhysicsT << 40 // << 41 //....oooOO0OOooo........oooOO0OOooo........oo << 42 //....oooOO0OOooo........oooOO0OOooo........oo << 43 28 44 #include "G4CrossSectionDataStore.hh" 29 #include "G4CrossSectionDataStore.hh" 45 #include "G4SystemOfUnits.hh" << 30 #include "G4HadronicException.hh" 46 #include "G4UnitsTable.hh" << 47 #include "Randomize.hh" << 48 #include "G4Nucleus.hh" << 49 << 50 #include "G4DynamicParticle.hh" << 51 #include "G4Isotope.hh" << 52 #include "G4Element.hh" << 53 #include "G4Material.hh" << 54 #include "G4MaterialTable.hh" << 55 #include "G4NistManager.hh" << 56 #include "G4HadronicParameters.hh" << 57 #include <algorithm> << 58 #include <typeinfo> << 59 << 60 //....oooOO0OOooo........oooOO0OOooo........oo << 61 << 62 G4CrossSectionDataStore::G4CrossSectionDataSto << 63 : nist(G4NistManager::Instance()) << 64 {} << 65 << 66 //....oooOO0OOooo........oooOO0OOooo........oo << 67 << 68 G4double << 69 G4CrossSectionDataStore::ComputeCrossSection(c << 70 const G4Material* mat) << 71 { << 72 currentMaterial = mat; << 73 matParticle = dp->GetDefinition(); << 74 matKinEnergy = dp->GetKineticEnergy(); << 75 matCrossSection = 0.0; << 76 << 77 std::size_t nElements = mat->GetNumberOfElem << 78 const G4double* nAtomsPerVolume = mat->GetVe << 79 << 80 if(xsecelm.size() < nElements) { xsecelm.res << 81 << 82 for(G4int i=0; i<(G4int)nElements; ++i) { << 83 G4double xs = << 84 nAtomsPerVolume[i]*GetCrossSection(dp, m << 85 matCrossSection += std::max(xs, 0.0); << 86 xsecelm[i] = matCrossSection; << 87 } << 88 return matCrossSection; << 89 } << 90 << 91 //....oooOO0OOooo........oooOO0OOooo........oo << 92 << 93 G4double G4CrossSectionDataStore::GetCrossSect << 94 << 95 << 96 { << 97 // first check the most last cross section << 98 G4int i = nDataSetList-1; << 99 G4int Z = elm->GetZasInt(); << 100 << 101 if(elm->GetNaturalAbundanceFlag() && << 102 dataSetList[i]->IsElementApplicable(dp, Z << 103 { << 104 // element wise cross section << 105 return dataSetList[i]->GetElementCrossSect << 106 } << 107 31 108 // isotope wise cross section << 109 G4int nIso = (G4int)elm->GetNumberOfIsotopes << 110 << 111 // user-defined isotope abundances << 112 const G4double* abundVector = elm->GetRelati << 113 << 114 G4double sigma = 0.0; << 115 << 116 // isotope and element wise cross sections << 117 for(G4int j = 0; j < nIso; ++j) << 118 { << 119 const G4Isotope* iso = elm->GetIsotope(j); << 120 sigma += abundVector[j] * << 121 GetIsoCrossSection(dp, Z, iso->GetN(), i << 122 } << 123 return sigma; << 124 } << 125 << 126 //....oooOO0OOooo........oooOO0OOooo........oo << 127 << 128 G4double << 129 G4CrossSectionDataStore::GetIsoCrossSection(co << 130 G4int Z, G4int A, << 131 const G4Isotope* iso, << 132 const G4Element* elm, << 133 const G4Material* mat, << 134 G4int idx) << 135 { << 136 // this methods is called after the check th << 137 // depend on isotopes, so first isotopes are << 138 if(dataSetList[idx]->IsIsoApplicable(dp, Z, << 139 return dataSetList[idx]->GetIsoCrossSectio << 140 } << 141 << 142 // no isotope wise cross section - check oth << 143 for (G4int j = nDataSetList-1; j >= 0; --j) << 144 if(dataSetList[j]->IsElementApplicable(dp, << 145 return dataSetList[j]->GetElementCrossSe << 146 } else if (dataSetList[j]->IsIsoApplicable << 147 return dataSetList[j]->GetIsoCrossSectio << 148 } << 149 } << 150 G4ExceptionDescription ed; << 151 ed << "No isotope cross section found for " << 152 << dp->GetDefinition()->GetParticleName() << 153 << " off target Element " << elm->GetName << 154 << " Z= " << Z << " A= " << A; << 155 if(nullptr != mat) ed << " from " << mat->Ge << 156 ed << " E(MeV)=" << dp->GetKineticEnergy()/M << 157 G4Exception("G4CrossSectionDataStore::GetIso << 158 FatalException, ed); << 159 return 0.0; << 160 } << 161 << 162 //....oooOO0OOooo........oooOO0OOooo........oo << 163 32 164 G4double 33 G4double 165 G4CrossSectionDataStore::GetCrossSection(const << 34 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle, 166 G4int << 35 const G4Element* anElement, 167 const G4Isotope* iso, << 36 G4double aTemperature) 168 const << 37 { 169 const G4Material* mat) << 38 if (NDataSetList == 0) 170 { << 39 { 171 for (G4int i = nDataSetList-1; i >= 0; --i) << 40 throw G4HadronicException(__FILE__, __LINE__, 172 if (dataSetList[i]->IsIsoApplicable(dp, Z, << 41 "G4CrossSectionDataStore: no data sets registered"); 173 return dataSetList[i]->GetIsoCrossSectio << 42 return DBL_MIN; 174 } else if(dataSetList[i]->IsElementApplica << 43 } 175 return dataSetList[i]->GetElementCrossSe << 44 for (G4int i = NDataSetList-1; i >= 0; i--) { 176 } << 45 if (DataSetList[i]->IsApplicable(aParticle, anElement)) 177 } << 46 return DataSetList[i]->GetCrossSection(aParticle, anElement, aTemperature); 178 G4ExceptionDescription ed; << 47 } 179 ed << "No isotope cross section found for " << 48 throw G4HadronicException(__FILE__, __LINE__, 180 << dp->GetDefinition()->GetParticleName() << 49 "G4CrossSectionDataStore: no applicable data set found " 181 << " off target Element " << elm->GetName << 50 "for particle/element"); 182 << " Z= " << Z << " A= " << A; << 51 return DBL_MIN; 183 if(nullptr != mat) ed << " from " << mat->Ge << 184 ed << " E(MeV)=" << dp->GetKineticEnergy()/M << 185 G4Exception("G4CrossSectionDataStore::GetCro << 186 FatalException, ed); << 187 return 0.0; << 188 } << 189 << 190 //....oooOO0OOooo........oooOO0OOooo........oo << 191 << 192 const G4Element* << 193 G4CrossSectionDataStore::SampleZandA(const G4D << 194 const G4M << 195 G4Nucleus& target) << 196 { << 197 if(nullptr != forcedElement) { return forced << 198 std::size_t nElements = mat->GetNumberOfElem << 199 const G4Element* anElement = mat->GetElement << 200 << 201 // select element from a compound << 202 if(1 < nElements) { << 203 G4double cross = matCrossSection*G4Uniform << 204 for(G4int i=0; i<(G4int)nElements; ++i) { << 205 if(cross <= xsecelm[i]) { << 206 anElement = mat->GetElement(i); << 207 break; << 208 } << 209 } << 210 } << 211 << 212 G4int Z = anElement->GetZasInt(); << 213 const G4Isotope* iso = nullptr; << 214 << 215 G4int i = nDataSetList-1; << 216 if (dataSetList[i]->IsElementApplicable(dp, << 217 << 218 //---------------------------------------- << 219 // element-wise cross section << 220 // isotope cross section is not computed << 221 //---------------------------------------- << 222 std::size_t nIso = anElement->GetNumberOfI << 223 iso = anElement->GetIsotope(0); << 224 << 225 // more than 1 isotope << 226 if(1 < nIso) { << 227 iso = dataSetList[i]->SelectIsotope(anEl << 228 dp-> << 229 dp->GetLogKineticEnergy()); << 230 } << 231 } else { << 232 << 233 //---------------------------------------- << 234 // isotope-wise cross section << 235 // isotope cross section is computed << 236 //---------------------------------------- << 237 std::size_t nIso = anElement->GetNumberOfI << 238 iso = anElement->GetIsotope(0); << 239 << 240 // more than 1 isotope << 241 if(1 < nIso) { << 242 const G4double* abundVector = anElement- << 243 if(xseciso.size() < nIso) { xseciso.resi << 244 << 245 G4double cross = 0.0; << 246 G4int j; << 247 for (j = 0; j<(G4int)nIso; ++j) { << 248 G4double xsec = 0.0; << 249 if(abundVector[j] > 0.0) { << 250 iso = anElement->GetIsotope(j); << 251 xsec = abundVector[j]* << 252 GetIsoCrossSection(dp, Z, iso->GetN(), i << 253 } << 254 cross += xsec; << 255 xseciso[j] = cross; << 256 } << 257 cross *= G4UniformRand(); << 258 for (j = 0; j<(G4int)nIso; ++j) { << 259 if(cross <= xseciso[j]) { << 260 iso = anElement->GetIsotope(j); << 261 break; << 262 } << 263 } << 264 } << 265 } << 266 target.SetIsotope(iso); << 267 return anElement; << 268 } 52 } 269 53 270 //....oooOO0OOooo........oooOO0OOooo........oo << 271 54 272 void 55 void 273 G4CrossSectionDataStore::BuildPhysicsTable(con << 56 G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet) 274 { 57 { 275 if (nDataSetList == 0) { << 58 if (NDataSetList == NDataSetMax) { 276 G4ExceptionDescription ed; << 59 G4cout << "WARNING: G4CrossSectionDataStore::AddDataSet: "<<G4endl; 277 ed << "No cross section is registered for << 60 G4cout << " reached maximum number of data sets"; 278 << part.GetParticleName() << G4endl; << 61 G4cout << " data set not added !!!!!!!!!!!!!!!!"; 279 G4Exception("G4CrossSectionDataStore::Buil << 62 return; 280 FatalException, ed); << 63 } 281 return; << 64 DataSetList[NDataSetList] = aDataSet; 282 } << 65 NDataSetList++; 283 matParticle = ∂ << 284 for (G4int i=0; i<nDataSetList; ++i) { << 285 dataSetList[i]->BuildPhysicsTable(part); << 286 } << 287 const G4MaterialTable* theMatTable = G4Mater << 288 std::size_t nelm = 0; << 289 std::size_t niso = 0; << 290 for(auto mat : *theMatTable) { << 291 std::size_t nElements = mat->GetNumberOfEl << 292 nelm = std::max(nelm, nElements); << 293 for(G4int j=0; j<(G4int)nElements; ++j) { << 294 niso = std::max(niso, mat->GetElement(j) << 295 } << 296 } << 297 // define vectors for a run << 298 xsecelm.resize(nelm, 0.0); << 299 xseciso.resize(niso, 0.0); << 300 } 66 } 301 67 302 //....oooOO0OOooo........oooOO0OOooo........oo << 303 << 304 void << 305 G4CrossSectionDataStore::DumpPhysicsTable(cons << 306 { << 307 // Print out all cross section data sets use << 308 // which they apply << 309 << 310 if (nDataSetList == 0) { << 311 G4cout << "WARNING - G4CrossSectionDataSto << 312 << " no data sets registered" << G4endl; << 313 return; << 314 } << 315 << 316 for (G4int i = nDataSetList-1; i >= 0; --i) << 317 G4double e1 = dataSetList[i]->GetMinKinEne << 318 G4double e2 = dataSetList[i]->GetMaxKinEne << 319 G4cout << 320 << " Cr_sctns: " << std::setw(25) << << 321 << G4BestUnit(e1, "Energy") << " ---> " << 322 << G4BestUnit(e2, "Energy") << "\n"; << 323 if (dataSetList[i]->GetName() == "G4CrossS << 324 dataSetList[i]->DumpPhysicsTable(part); << 325 G4cout << G4endl; << 326 } << 327 } << 328 } << 329 << 330 //....oooOO0OOooo........oooOO0OOooo........oo << 331 << 332 void G4CrossSectionDataStore::DumpHtml(const G << 333 std::of << 334 { << 335 // Write cross section data set info to html << 336 // documentation page << 337 << 338 G4double ehi = 0; << 339 G4double elo = 0; << 340 auto param = G4HadronicParameters::Instance( << 341 G4String physListName = param->GetPhysListNa << 342 G4String dirName = param->GetPhysListDocDir( << 343 68 344 for (G4int i = nDataSetList-1; i > 0; i--) { << 69 void 345 elo = dataSetList[i]->GetMinKinEnergy()/Ge << 70 G4CrossSectionDataStore:: 346 ehi = dataSetList[i]->GetMaxKinEnergy()/Ge << 71 BuildPhysicsTable(const G4ParticleDefinition& aParticleType) 347 outFile << " <li><b><a href=\"" << ph << 348 << dataSetList[i]->GetName() << ".html\" << 349 << dataSetList[i]->GetName() << "< << 350 << elo << " GeV to " << ehi << " G << 351 PrintCrossSectionHtml(dataSetList[i], phys << 352 } << 353 << 354 G4double defaultHi = dataSetList[0]->GetMaxK << 355 if (ehi < defaultHi) { << 356 outFile << " <li><b><a href=\"" << da << 357 << ".html\"> " << 358 << dataSetList[0]->GetName() << "< << 359 << ehi << " GeV to " << defaultHi << 360 PrintCrossSectionHtml(dataSetList[0], phys << 361 } << 362 } << 363 << 364 //....oooOO0OOooo........oooOO0OOooo........oo << 365 << 366 void G4CrossSectionDataStore::PrintCrossSectio << 367 << 368 << 369 { << 370 << 371 G4String pathName = dirName + "/" + physList << 372 std::ofstream outCS; << 373 outCS.open(pathName); << 374 outCS << "<html>\n"; << 375 outCS << "<head>\n"; << 376 outCS << "<title>Description of " << cs->Get << 377 << "</title>\n"; << 378 outCS << "</head>\n"; << 379 outCS << "<body>\n"; << 380 << 381 cs->CrossSectionDescription(outCS); << 382 << 383 outCS << "</body>\n"; << 384 outCS << "</html>\n"; << 385 } << 386 << 387 //....oooOO0OOooo........oooOO0OOooo........oo << 388 << 389 G4String G4CrossSectionDataStore::HtmlFileName << 390 { << 391 G4String str(in); << 392 // replace blanks by _ C++11 version: << 393 std::transform(str.begin(), str.end(), str. << 394 return ch == ' ' ? '_' : ch; << 395 }); << 396 str=str + ".html"; << 397 return str; << 398 } << 399 << 400 //....oooOO0OOooo........oooOO0OOooo........oo << 401 << 402 void G4CrossSectionDataStore::AddDataSet(G4VCr << 403 { 72 { 404 if(p->ForAllAtomsAndEnergies()) { << 73 if (NDataSetList == 0) 405 dataSetList.clear(); << 74 { 406 nDataSetList = 0; << 75 G4Exception("G4CrossSectionDataStore", "007", FatalException, 407 } << 76 "BuildPhysicsTable: no data sets registered"); 408 dataSetList.push_back(p); << 77 return; 409 ++nDataSetList; << 78 } >> 79 for (G4int i = NDataSetList-1; i >= 0; i--) { >> 80 DataSetList[i]->BuildPhysicsTable(aParticleType); >> 81 } 410 } 82 } 411 83 412 //....oooOO0OOooo........oooOO0OOooo........oo << 413 84 414 void G4CrossSectionDataStore::AddDataSet(G4VCr << 85 void >> 86 G4CrossSectionDataStore:: >> 87 DumpPhysicsTable(const G4ParticleDefinition& aParticleType) 415 { 88 { 416 if(p->ForAllAtomsAndEnergies()) { << 89 if (NDataSetList == 0) { 417 dataSetList.clear(); << 90 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: no data sets registered"<<G4endl; 418 dataSetList.push_back(p); << 91 return; 419 nDataSetList = 1; << 92 } 420 } else if ( i >= dataSetList.size() ) { << 93 for (G4int i = NDataSetList-1; i >= 0; i--) { 421 dataSetList.push_back(p); << 94 DataSetList[i]->DumpPhysicsTable(aParticleType); 422 ++nDataSetList; << 95 } 423 } else { << 424 std::vector< G4VCrossSectionDataSet* >::it << 425 dataSetList.insert(it , p); << 426 ++nDataSetList; << 427 } << 428 } 96 } 429 << 430 //....oooOO0OOooo........oooOO0OOooo........oo << 431 97