| 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