| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
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
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........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
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
164 G4double
165 G4CrossSectionDataStore::GetCrossSection(const
166 G4int
167 const G4Isotope* iso,
168 const
169 const G4Material* mat)
170 {
171 for (G4int i = nDataSetList-1; i >= 0; --i)
172 if (dataSetList[i]->IsIsoApplicable(dp, Z,
173 return dataSetList[i]->GetIsoCrossSectio
174 } else if(dataSetList[i]->IsElementApplica
175 return dataSetList[i]->GetElementCrossSe
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->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 }
269
270 //....oooOO0OOooo........oooOO0OOooo........oo
271
272 void
273 G4CrossSectionDataStore::BuildPhysicsTable(con
274 {
275 if (nDataSetList == 0) {
276 G4ExceptionDescription ed;
277 ed << "No cross section is registered for
278 << part.GetParticleName() << G4endl;
279 G4Exception("G4CrossSectionDataStore::Buil
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 = 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 }
301
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
344 for (G4int i = nDataSetList-1; i > 0; i--) {
345 elo = dataSetList[i]->GetMinKinEnergy()/Ge
346 ehi = dataSetList[i]->GetMaxKinEnergy()/Ge
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 {
404 if(p->ForAllAtomsAndEnergies()) {
405 dataSetList.clear();
406 nDataSetList = 0;
407 }
408 dataSetList.push_back(p);
409 ++nDataSetList;
410 }
411
412 //....oooOO0OOooo........oooOO0OOooo........oo
413
414 void G4CrossSectionDataStore::AddDataSet(G4VCr
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* >::it
425 dataSetList.insert(it , p);
426 ++nDataSetList;
427 }
428 }
429
430 //....oooOO0OOooo........oooOO0OOooo........oo
431