| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
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 // 7 //
>> 8 // $Id: G4CrossSectionDataStore.cc,v 1.2 1999/12/15 14:52:06 gunter Exp $
>> 9 // GEANT4 tag $Name: geant4-03-01 $
26 // 10 //
27 // ------------------------------------------- <<
28 // 11 //
29 // GEANT4 Class file << 12 // GEANT4 physics class: G4CrossSectionDataStore
>> 13 // F.W. Jones, TRIUMF, 19-NOV-97
30 // 14 //
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 15
44 #include "G4CrossSectionDataStore.hh" 16 #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 17
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 18
164 G4double 19 G4double
165 G4CrossSectionDataStore::GetCrossSection(const << 20 G4CrossSectionDataStore::GetCrossSection(const G4DynamicParticle* aParticle,
166 G4int << 21 const G4Element* anElement)
167 const G4Isotope* iso, <<
168 const <<
169 const G4Material* mat) <<
170 { 22 {
171 for (G4int i = nDataSetList-1; i >= 0; --i) << 23 if (NDataSetList == 0) {
172 if (dataSetList[i]->IsIsoApplicable(dp, Z, << 24 G4Exception("G4CrossSectionDataStore: no data sets registered");
173 return dataSetList[i]->GetIsoCrossSectio << 25 return DBL_MIN;
174 } else if(dataSetList[i]->IsElementApplica << 26 }
175 return dataSetList[i]->GetElementCrossSe << 27 for (G4int i = NDataSetList-1; i >= 0; i--) {
176 } << 28 if (DataSetList[i]->IsApplicable(aParticle, anElement))
177 } << 29 return DataSetList[i]->GetCrossSection(aParticle, anElement);
178 G4ExceptionDescription ed; << 30 }
179 ed << "No isotope cross section found for " << 31 G4Exception("G4CrossSectionDataStore: no applicable data set found "
180 << dp->GetDefinition()->GetParticleName() << 32 "for particle/element");
181 << " off target Element " << elm->GetName << 33 return DBL_MIN;
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 } 34 }
189 35
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 36
272 void 37 void
273 G4CrossSectionDataStore::BuildPhysicsTable(con << 38 G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet)
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 { 39 {
370 << 40 if (NDataSetList == NDataSetMax) {
371 G4String pathName = dirName + "/" + physList << 41 G4Exception("G4CrossSectionDataStore::AddDataSet: "
372 std::ofstream outCS; << 42 "reached maximum number of data sets");
373 outCS.open(pathName); << 43 return;
374 outCS << "<html>\n"; << 44 }
375 outCS << "<head>\n"; << 45 DataSetList[NDataSetList] = aDataSet;
376 outCS << "<title>Description of " << cs->Get << 46 NDataSetList++;
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 } 47 }
399 48
400 //....oooOO0OOooo........oooOO0OOooo........oo <<
401 49
402 void G4CrossSectionDataStore::AddDataSet(G4VCr << 50 void
>> 51 G4CrossSectionDataStore::
>> 52 BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
403 { 53 {
404 if(p->ForAllAtomsAndEnergies()) { << 54 if (NDataSetList == 0) {
405 dataSetList.clear(); << 55 G4Exception("G4CrossSectionDataStore: no data sets registered");
406 nDataSetList = 0; << 56 return;
407 } << 57 }
408 dataSetList.push_back(p); << 58 for (G4int i = NDataSetList-1; i >= 0; i--) {
409 ++nDataSetList; << 59 DataSetList[i]->BuildPhysicsTable(aParticleType);
>> 60 }
410 } 61 }
411 62
412 //....oooOO0OOooo........oooOO0OOooo........oo <<
413 63
414 void G4CrossSectionDataStore::AddDataSet(G4VCr << 64 void
>> 65 G4CrossSectionDataStore::
>> 66 DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
415 { 67 {
416 if(p->ForAllAtomsAndEnergies()) { << 68 if (NDataSetList == 0) {
417 dataSetList.clear(); << 69 G4Exception("G4CrossSectionDataStore: no data sets registered");
418 dataSetList.push_back(p); << 70 return;
419 nDataSetList = 1; << 71 }
420 } else if ( i >= dataSetList.size() ) { << 72 for (G4int i = NDataSetList-1; i >= 0; i--) {
421 dataSetList.push_back(p); << 73 DataSetList[i]->DumpPhysicsTable(aParticleType);
422 ++nDataSetList; << 74 }
423 } else { <<
424 std::vector< G4VCrossSectionDataSet* >::it <<
425 dataSetList.insert(it , p); <<
426 ++nDataSetList; <<
427 } <<
428 } 75 }
429 <<
430 //....oooOO0OOooo........oooOO0OOooo........oo <<
431 76