| 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 // Class Description
27 // Cross-section data set for a high precision
28 // libraries) description of elastic scatterin
29 // Class Description - End
30
31 // 15-Nov-06 First Implementation is done by T
32 // P. Arce, June-2014 Conversion neutron_hp to
33 //
34 #include "G4ParticleHPJENDLHEData.hh"
35
36 #include "G4ElementTable.hh"
37 #include "G4ParticleHPData.hh"
38 #include "G4PhysicsFreeVector.hh"
39 #include "G4Pow.hh"
40 #include "G4SystemOfUnits.hh"
41
42 G4bool G4ParticleHPJENDLHEData::IsApplicable(c
43 {
44 G4bool result = true;
45 G4double eKin = aP->GetKineticEnergy();
46 // if(eKin>20*MeV||aP->GetDefinition()!=G4Ne
47 if (eKin < 20 * MeV || 3 * GeV < eKin || aP-
48 result = false;
49 }
50 // Element Check
51 else if (!(vElement[anE->GetIndex()]))
52 result = false;
53
54 return result;
55 }
56
57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa
58 {
59 for (auto& itZ : mIsotope) {
60 std::map<G4int, G4PhysicsVector*>* pointer
61 if (pointer_map != nullptr) {
62 for (auto& itA : *pointer_map) {
63 G4PhysicsVector* pointerPhysicsVector
64 if (pointerPhysicsVector != nullptr) {
65 delete pointerPhysicsVector;
66 itA.second = NULL;
67 }
68 }
69 delete pointer_map;
70 itZ.second = NULL;
71 }
72 }
73 mIsotope.clear();
74 }
75
76 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa
77 : G4VCrossSectionDataSet("JENDLHE" + reactio
78 {
79 reactionName = reaction;
80 BuildPhysicsTable(*pd);
81 }
82
83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHED
84
85 void G4ParticleHPJENDLHEData::BuildPhysicsTabl
86 {
87 particleName = aP.GetParticleName();
88
89 const G4String& baseName = G4FindDataDir("G4
90 const G4String& dirName = baseName + "/JENDL
91 const G4String& aFSType = "/CrossSection/";
92 G4ParticleHPNames theNames;
93
94 G4String filename;
95
96 // Create JENDL_HE data
97 // Create map element or isotope
98
99 std::size_t numberOfElements = G4Element::Ge
100
101 // make a PhysicsVector for each element
102
103 auto theElementTable = G4Element::GetElement
104 vElement.clear();
105 vElement.resize(numberOfElements);
106 for (std::size_t i = 0; i < numberOfElements
107 G4Element* theElement = (*theElementTable)
108 vElement[i] = false;
109
110 // isotope
111 auto nIso = (G4int)(*theElementTable)[i]->
112 auto Z = (G4int)(*theElementTable)[i]->Get
113 for (G4int i1 = 0; i1 < nIso; ++i1) {
114 G4int A = theElement->GetIsotope(i1)->Ge
115
116 if (isThisNewIsotope(Z, A)) {
117 std::stringstream ss;
118 ss << dirName << aFSType << Z << "_" << A <<
119 filename = ss.str();
120 std::fstream file;
121 file.open(filename, std::fstream::in);
122 G4int dummy;
123 file >> dummy;
124 if (file.good()) {
125 vElement[i] = true;
126
127 // read the file
128 G4PhysicsVector* aPhysVec = readAFile(&fil
129 registAPhysicsVector(Z, A, aPhysVec);
130 }
131 file.close();
132 }
133 }
134 }
135 }
136
137 void G4ParticleHPJENDLHEData::DumpPhysicsTable
138 {}
139
140 G4double G4ParticleHPJENDLHEData::GetCrossSect
141
142 {
143 // Primary energy >20MeV
144 // Thus not taking into account of Doppler b
145 // also not taking into account of Target th
146
147 G4double result = 0;
148
149 G4double ek = aP->GetKineticEnergy();
150
151 auto nIso = (G4int)anE->GetNumberOfIsotopes(
152 auto Z = (G4int)anE->GetZ();
153 for (G4int i1 = 0; i1 < nIso; ++i1) {
154 G4int A = anE->GetIsotope(i1)->GetN();
155 G4double frac = anE->GetRelativeAbundanceV
156 // This case does NOT request "*perCent".
157 result += frac * getXSfromThisIsotope(Z, A
158 }
159 return result;
160 }
161
162 G4PhysicsVector* G4ParticleHPJENDLHEData::read
163 {
164 G4int dummy;
165 G4int len;
166 *file >> dummy;
167 *file >> len;
168
169 std::vector<G4double> v_e;
170 std::vector<G4double> v_xs;
171
172 for (G4int i = 0; i < len; ++i) {
173 G4double e;
174 G4double xs;
175
176 *file >> e;
177 *file >> xs;
178 // data are written in eV and barn.
179 v_e.push_back(e * eV);
180 v_xs.push_back(xs * barn);
181 }
182
183 auto aPhysVec = new G4PhysicsFreeVector(stat
184
185 for (G4int i = 0; i < len; ++i) {
186 aPhysVec->PutValues(static_cast<std::size_
187 }
188
189 return aPhysVec;
190 }
191
192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4
193 {
194 if (mIsotope.find(z) == mIsotope.end()) retu
195 if (mIsotope.find(z)->second->find(a) == mIs
196 return true;
197 }
198
199 void G4ParticleHPJENDLHEData::registAPhysicsVe
200 {
201 std::pair<G4int, G4PhysicsVector*> aPair = s
202 auto itm = mIsotope.find(Z);
203 if (itm != mIsotope.cend()) {
204 itm->second->insert(aPair);
205 }
206 else {
207 auto aMap = new std::map<G4int, G4PhysicsV
208 aMap->insert(aPair);
209 mIsotope.insert(std::pair<G4int, std::map<
210 }
211 }
212
213 G4double G4ParticleHPJENDLHEData::getXSfromThi
214 {
215 G4double aXSection = 0.0;
216
217 G4PhysicsVector* aPhysVec;
218 if (mIsotope.find(Z)->second->find(A) != mIs
219 aPhysVec = mIsotope.find(Z)->second->find(
220 aXSection = aPhysVec->Value(ek);
221 }
222 else {
223 // Select closest one in the same Z
224 G4int delta0 = 99; // no mean for 99
225 for (auto it = mIsotope.find(Z)->second->c
226 {
227 G4int delta = std::abs(A - it->first);
228 if (delta < delta0) delta0 = delta;
229 }
230
231 // Randomize of selection larger or smalle
232 if (G4UniformRand() < 0.5) delta0 *= -1;
233 G4int A1 = A + delta0;
234 if (mIsotope.find(Z)->second->find(A1) !=
235 aPhysVec = mIsotope.find(Z)->second->fin
236 }
237 else {
238 A1 = A - delta0;
239 aPhysVec = mIsotope.find(Z)->second->fin
240 }
241
242 aXSection = aPhysVec->Value(ek);
243 // X^(2/3) factor
244 aXSection *= G4Pow::GetInstance()->A23(1.0
245 }
246
247 return aXSection;
248 }
249