| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4CrossSectionDataStore
33 //
34 // Modifications:
35 // 23.01.2009 V.Ivanchenko add destruction of data sets
36 // 29.04.2010 G.Folger modifictaions for integer A & Z
37 // 14.03.2011 V.Ivanchenko fixed DumpPhysicsTable
38 // 15.08.2011 G.Folger, V.Ivanchenko, T.Koi, D.Wright redesign the class
39 // 07.03.2013 M.Maire cosmetic in DumpPhysicsTable
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
43
44 #include "G4CrossSectionDataStore.hh"
45 #include "G4SystemOfUnits.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........oooOO0OOooo........oooOO0OOooo.....
61
62 G4CrossSectionDataStore::G4CrossSectionDataStore()
63 : nist(G4NistManager::Instance())
64 {}
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
67
68 G4double
69 G4CrossSectionDataStore::ComputeCrossSection(const G4DynamicParticle* dp,
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->GetNumberOfElements();
78 const G4double* nAtomsPerVolume = mat->GetVecNbOfAtomsPerVolume();
79
80 if(xsecelm.size() < nElements) { xsecelm.resize(nElements); }
81
82 for(G4int i=0; i<(G4int)nElements; ++i) {
83 G4double xs =
84 nAtomsPerVolume[i]*GetCrossSection(dp, mat->GetElement(i), mat);
85 matCrossSection += std::max(xs, 0.0);
86 xsecelm[i] = matCrossSection;
87 }
88 return matCrossSection;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
92
93 G4double G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* dp,
94 const G4Element* elm,
95 const G4Material* mat)
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, mat))
103 {
104 // element wise cross section
105 return dataSetList[i]->GetElementCrossSection(dp, Z, mat);
106 }
107
108 // isotope wise cross section
109 G4int nIso = (G4int)elm->GetNumberOfIsotopes();
110
111 // user-defined isotope abundances
112 const G4double* abundVector = elm->GetRelativeAbundanceVector();
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(), iso, elm, mat, i);
122 }
123 return sigma;
124 }
125
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
127
128 G4double
129 G4CrossSectionDataStore::GetIsoCrossSection(const G4DynamicParticle* dp,
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 that dataSetList[idx]
137 // depend on isotopes, so first isotopes are checked
138 if(dataSetList[idx]->IsIsoApplicable(dp, Z, A, elm, mat) ) {
139 return dataSetList[idx]->GetIsoCrossSection(dp, Z, A, iso, elm, mat);
140 }
141
142 // no isotope wise cross section - check other datasets
143 for (G4int j = nDataSetList-1; j >= 0; --j) {
144 if(dataSetList[j]->IsElementApplicable(dp, Z, mat)) {
145 return dataSetList[j]->GetElementCrossSection(dp, Z, mat);
146 } else if (dataSetList[j]->IsIsoApplicable(dp, Z, A, elm, mat)) {
147 return dataSetList[j]->GetIsoCrossSection(dp, Z, A, iso, elm, mat);
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->GetName();
156 ed << " E(MeV)=" << dp->GetKineticEnergy()/MeV << G4endl;
157 G4Exception("G4CrossSectionDataStore::GetIsoCrossSection", "had001",
158 FatalException, ed);
159 return 0.0;
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
163
164 G4double
165 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* dp,
166 G4int Z, G4int A,
167 const G4Isotope* iso,
168 const G4Element* elm,
169 const G4Material* mat)
170 {
171 for (G4int i = nDataSetList-1; i >= 0; --i) {
172 if (dataSetList[i]->IsIsoApplicable(dp, Z, A, elm, mat) ) {
173 return dataSetList[i]->GetIsoCrossSection(dp, Z, A, iso, elm, mat);
174 } else if(dataSetList[i]->IsElementApplicable(dp, Z, mat)) {
175 return dataSetList[i]->GetElementCrossSection(dp, Z, mat);
176 }
177 }
178 G4ExceptionDescription ed;
179 ed << "No isotope cross section found for "
180 << dp->GetDefinition()->GetParticleName()
181 << " off target Element " << elm->GetName()
182 << " Z= " << Z << " A= " << A;
183 if(nullptr != mat) ed << " from " << mat->GetName();
184 ed << " E(MeV)=" << dp->GetKineticEnergy()/MeV << G4endl;
185 G4Exception("G4CrossSectionDataStore::GetCrossSection", "had001",
186 FatalException, ed);
187 return 0.0;
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
191
192 const G4Element*
193 G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* dp,
194 const G4Material* mat,
195 G4Nucleus& target)
196 {
197 if(nullptr != forcedElement) { return forcedElement; }
198 std::size_t nElements = mat->GetNumberOfElements();
199 const G4Element* anElement = mat->GetElement(0);
200
201 // select element from a compound
202 if(1 < nElements) {
203 G4double cross = matCrossSection*G4UniformRand();
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, Z, mat)) {
217
218 //----------------------------------------------------------------
219 // element-wise cross section
220 // isotope cross section is not computed
221 //----------------------------------------------------------------
222 std::size_t nIso = anElement->GetNumberOfIsotopes();
223 iso = anElement->GetIsotope(0);
224
225 // more than 1 isotope
226 if(1 < nIso) {
227 iso = dataSetList[i]->SelectIsotope(anElement,
228 dp->GetKineticEnergy(),
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->GetNumberOfIsotopes();
238 iso = anElement->GetIsotope(0);
239
240 // more than 1 isotope
241 if(1 < nIso) {
242 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
243 if(xseciso.size() < nIso) { xseciso.resize(nIso); }
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(), iso, anElement, mat, 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 }
269
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
271
272 void
273 G4CrossSectionDataStore::BuildPhysicsTable(const G4ParticleDefinition& part)
274 {
275 if (nDataSetList == 0) {
276 G4ExceptionDescription ed;
277 ed << "No cross section is registered for "
278 << part.GetParticleName() << G4endl;
279 G4Exception("G4CrossSectionDataStore::BuildPhysicsTable", "had001",
280 FatalException, ed);
281 return;
282 }
283 matParticle = ∂
284 for (G4int i=0; i<nDataSetList; ++i) {
285 dataSetList[i]->BuildPhysicsTable(part);
286 }
287 const G4MaterialTable* theMatTable = G4Material::GetMaterialTable();
288 std::size_t nelm = 0;
289 std::size_t niso = 0;
290 for(auto mat : *theMatTable) {
291 std::size_t nElements = mat->GetNumberOfElements();
292 nelm = std::max(nelm, nElements);
293 for(G4int j=0; j<(G4int)nElements; ++j) {
294 niso = std::max(niso, mat->GetElement(j)->GetNumberOfIsotopes());
295 }
296 }
297 // define vectors for a run
298 xsecelm.resize(nelm, 0.0);
299 xseciso.resize(niso, 0.0);
300 }
301
302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
303
304 void
305 G4CrossSectionDataStore::DumpPhysicsTable(const G4ParticleDefinition& part)
306 {
307 // Print out all cross section data sets used and the energies at
308 // which they apply
309
310 if (nDataSetList == 0) {
311 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: "
312 << " no data sets registered" << G4endl;
313 return;
314 }
315
316 for (G4int i = nDataSetList-1; i >= 0; --i) {
317 G4double e1 = dataSetList[i]->GetMinKinEnergy();
318 G4double e2 = dataSetList[i]->GetMaxKinEnergy();
319 G4cout
320 << " Cr_sctns: " << std::setw(25) << dataSetList[i]->GetName() << ": "
321 << G4BestUnit(e1, "Energy") << " ---> "
322 << G4BestUnit(e2, "Energy") << "\n";
323 if (dataSetList[i]->GetName() == "G4CrossSectionPairGG") {
324 dataSetList[i]->DumpPhysicsTable(part);
325 G4cout << G4endl;
326 }
327 }
328 }
329
330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
331
332 void G4CrossSectionDataStore::DumpHtml(const G4ParticleDefinition& /* pD */,
333 std::ofstream& outFile) const
334 {
335 // Write cross section data set info to html physics list
336 // documentation page
337
338 G4double ehi = 0;
339 G4double elo = 0;
340 auto param = G4HadronicParameters::Instance();
341 G4String physListName = param->GetPhysListName();
342 G4String dirName = param->GetPhysListDocDir();
343
344 for (G4int i = nDataSetList-1; i > 0; i--) {
345 elo = dataSetList[i]->GetMinKinEnergy()/GeV;
346 ehi = dataSetList[i]->GetMaxKinEnergy()/GeV;
347 outFile << " <li><b><a href=\"" << physListName << "_"
348 << dataSetList[i]->GetName() << ".html\"> "
349 << dataSetList[i]->GetName() << "</a> from "
350 << elo << " GeV to " << ehi << " GeV </b></li>\n";
351 PrintCrossSectionHtml(dataSetList[i], physListName, dirName);
352 }
353
354 G4double defaultHi = dataSetList[0]->GetMaxKinEnergy()/GeV;
355 if (ehi < defaultHi) {
356 outFile << " <li><b><a href=\"" << dataSetList[0]->GetName()
357 << ".html\"> "
358 << dataSetList[0]->GetName() << "</a> from "
359 << ehi << " GeV to " << defaultHi << " GeV </b></li>\n";
360 PrintCrossSectionHtml(dataSetList[0], physListName, dirName);
361 }
362 }
363
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
365
366 void G4CrossSectionDataStore::PrintCrossSectionHtml(const G4VCrossSectionDataSet *cs,
367 const G4String& physListName,
368 const G4String& dirName) const
369 {
370
371 G4String pathName = dirName + "/" + physListName + "_" + HtmlFileName(cs->GetName());
372 std::ofstream outCS;
373 outCS.open(pathName);
374 outCS << "<html>\n";
375 outCS << "<head>\n";
376 outCS << "<title>Description of " << cs->GetName()
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........oooOO0OOooo........oooOO0OOooo.....
388
389 G4String G4CrossSectionDataStore::HtmlFileName(const G4String & in) const
390 {
391 G4String str(in);
392 // replace blanks by _ C++11 version:
393 std::transform(str.begin(), str.end(), str.begin(), [](char ch) {
394 return ch == ' ' ? '_' : ch;
395 });
396 str=str + ".html";
397 return str;
398 }
399
400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
401
402 void G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* p)
403 {
404 if(p->ForAllAtomsAndEnergies()) {
405 dataSetList.clear();
406 nDataSetList = 0;
407 }
408 dataSetList.push_back(p);
409 ++nDataSetList;
410 }
411
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
413
414 void G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* p, std::size_t i)
415 {
416 if(p->ForAllAtomsAndEnergies()) {
417 dataSetList.clear();
418 dataSetList.push_back(p);
419 nDataSetList = 1;
420 } else if ( i >= dataSetList.size() ) {
421 dataSetList.push_back(p);
422 ++nDataSetList;
423 } else {
424 std::vector< G4VCrossSectionDataSet* >::iterator it = dataSetList.end() - i;
425 dataSetList.insert(it , p);
426 ++nDataSetList;
427 }
428 }
429
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
431