| 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 source file
30 //
31 // File name: G4ParticleHPIsoProbabilityT
32 //
33 // Authors: Marek Zmeskal (CTU, Czech Tec
34 // Loic Thulliez (CEA France)
35 //
36 // Creation date: 4 June 2024
37 //
38 // Description: Class for the probability
39 // and for the given tempera
40 // It reads the files with p
41 // finds the correct cross-s
42 //
43 // Modifications:
44 //
45 // -------------------------------------------
46 //
47 //
48
49 #include "G4ParticleHPIsoProbabilityTable_NJOY
50 #include "G4SystemOfUnits.hh"
51 #include "G4DynamicParticle.hh"
52 #include "G4ParticleHPVector.hh"
53 #include "Randomize.hh"
54 #include "G4NucleiProperties.hh"
55 #include "G4ReactionProduct.hh"
56 #include "G4ParticleHPManager.hh"
57 #include "G4ParticleHPChannel.hh"
58 #include "G4ParticleHPChannelList.hh"
59 #include "G4Nucleus.hh"
60 #include "G4Element.hh"
61
62 #include <string>
63 #include <sstream>
64
65 ///-------------------------------------------
66 G4ParticleHPIsoProbabilityTable_NJOY::G4Partic
67 tableOrder = 0;
68 lssf_flag = -1;
69 }
70
71 ///-------------------------------------------
72 G4ParticleHPIsoProbabilityTable_NJOY::~G4Parti
73
74 ///-------------------------------------------
75 void G4ParticleHPIsoProbabilityTable_NJOY::Ini
76 Z = theZ;
77 A = theA;
78 m = them;
79 T = theT;
80 G4cout << "The NJOY probability tables are b
81 std::string strZ = std::to_string(Z);
82 std::string strA = std::to_string(A);
83 filename = strZ + "_" + strA;
84 if ( m != 0 ) {
85 std::string strm = std::to_string(m);
86 filename += "_m" + strm;
87 }
88 G4String fullPathFileName = dirName + filena
89 std::istringstream theData( std::ios::in );
90 G4ParticleHPManager::GetInstance()->GetDataS
91 if ( theData.good() ) {
92 G4double emin;
93 G4double emax;
94 G4double energy;
95 theData >> emin >> emax;
96 Emin = emin * eV;
97 Emax = emax * eV;
98 theData >> nEnergies >> tableOrder >> lssf
99 theEnergies = new G4ParticleHPVector(nEner
100 theProbabilities = new std::vector< std::v
101 theElasticData = new std::vector< std::vec
102 theCaptureData = new std::vector< std::vec
103 theFissionData = new std::vector< std::vec
104 G4double probability, total, elastic, capt
105 for ( G4int i = 0; i < nEnergies; i++ ) {
106 theData >> energy;
107 theEnergies->SetEnergy( i, energy * eV )
108 std::vector< G4double >* vecprob = new s
109 std::vector< G4double >* vecela = new s
110 std::vector< G4double >* veccap = new s
111 std::vector< G4double >* vecfiss = new s
112 for ( G4int j = 0; j < tableOrder; j++ )
113 theData >> probability >> total >> ela
114 vecprob->push_back( probability );
115 vecela->push_back( elastic );
116 veccap->push_back( capture );
117 vecfiss->push_back( fission );
118 }
119 theProbabilities->push_back( vecprob );
120 theElasticData->push_back( vecela );
121 theCaptureData->push_back( veccap );
122 theFissionData->push_back( vecfiss );
123 }
124 G4cout << "Probability tables found and su
125 } else {
126 G4cout << "No probability tables found for
127 }
128 }
129
130 ///-------------------------------------------
131 G4double G4ParticleHPIsoProbabilityTable_NJOY:
132 const G4Element* ele, G4double& kineticEnerg
133 if ( kineticEnergy == energy_cache[id] ) {
134 if ( MTnumber == 2 ) { // elasti
135 return xsela_cache[id];
136 } else if ( MTnumber == 102 ) { // radiat
137 return xscap_cache[id];
138 } else if ( MTnumber == 18 ) { // fissio
139 return xsfiss_cache[id];
140 }
141 }
142 energy_cache[id] = kineticEnergy;
143 if ( kineticEnergy < Emin || kineticEnergy
144 // if the kinetic energy is outside of the
145 G4int indexEl = (G4int)ele->GetIndex();
146 G4int isotopeJ = -1; // index of isotope
147 for ( G4int j = 0; j < (G4int)ele->GetNumb
148 if ( A == (G4int)ele->GetIsotope(j)->Get
149 isotopeJ = j;
150 break;
151 }
152 }
153 G4double frac = ele->GetRelativeAbundanceV
154 G4double weightedelasticXS;
155 G4double weightedcaptureXS;
156 if ( G4ParticleHPManager::GetInstance()->G
157 weightedelasticXS = (*G4ParticleHPManage
158 weightedcaptureXS = (*G4ParticleHPManage
159 } else {
160 weightedelasticXS = this->GetDopplerBroa
161 weightedcaptureXS = this->GetDopplerBroa
162 }
163 xsela_cache[id] = weightedelasticXS / frac
164 xscap_cache[id] = weightedcaptureXS / frac
165 if ( Z < 88 ) {
166 xsfiss_cache[id] = 0.0;
167 } else {
168 if ( G4ParticleHPManager::GetInstance()-
169 G4double weightedfissionXS = (*G4Parti
170 xsfiss_cache[id] = weightedfissionXS /
171 } else {
172 G4double weightedfissionXS = this->Get
173 xsfiss_cache[id] = weightedfissionXS /
174 }
175 }
176 } else if ( lssf_flag == 0 ) {
177 G4int indexE = theEnergies->GetEnergyIndex
178 std::vector< G4double >* theProbability1 =
179 std::vector< G4double >* theProbability2 =
180 G4double ene1 = theEnergies->GetEnergy(ind
181 G4double ene2 = theEnergies->GetEnergy(ind
182 G4double rand = random_number;
183 G4int indexP1;
184 G4int indexP2;
185 for ( indexP1 = 0; indexP1 < tableOrder; i
186 if ( rand <= theProbability1->at(indexP1
187 }
188 for ( indexP2 = 0; indexP2 < tableOrder; i
189 if ( rand <= theProbability2->at(indexP2
190 }
191 G4double ela1 = theElasticData->at(indexE
192 G4double ela2 = theElasticData->at(indexE)
193 G4double cap1 = theCaptureData->at(indexE
194 G4double cap2 = theCaptureData->at(indexE)
195 xsela_cache[id] = theInt.Lin( kineticEnerg
196 xscap_cache[id] = theInt.Lin( kineticEnerg
197 if ( Z < 88 ) {
198 xsfiss_cache[id] = 0.0;
199 } else {
200 G4double fiss1 = theFissionData->at(inde
201 G4double fiss2 = theFissionData->at(inde
202 xsfiss_cache[id] = theInt.Lin( kineticEn
203 }
204 } else if ( lssf_flag == 1 ) {
205 G4int indexE = theEnergies->GetEnergyIndex
206 std::vector< G4double >* theProbability1 =
207 std::vector< G4double >* theProbability2 =
208 G4double ene1 = theEnergies->GetEnergy(ind
209 G4double ene2 = theEnergies->GetEnergy(ind
210 G4double rand = random_number;
211 G4int indexP1;
212 G4int indexP2;
213 for ( indexP1 = 0; indexP1 < tableOrder; i
214 if ( rand <= theProbability1->at(indexP1
215 }
216 for ( indexP2 = 0; indexP2 < tableOrder; i
217 if ( rand <= theProbability2->at(indexP2
218 }
219 G4int indexEl = (G4int)ele->GetIndex();
220 G4int isotopeJ = -1; // index of isotope
221 for ( G4int j = 0; j < (G4int)ele->GetNumb
222 if ( A == (G4int)ele->GetIsotope(j)->Get
223 isotopeJ = j;
224 break;
225 }
226 }
227 G4double frac = ele->GetRelativeAbundanceV
228 G4double weightedelasticXS;
229 G4double weightedcaptureXS;
230 if ( G4ParticleHPManager::GetInstance()->G
231 weightedelasticXS = (*G4ParticleHPManage
232 weightedcaptureXS = (*G4ParticleHPManage
233 } else {
234 weightedelasticXS = this->GetDopplerBroa
235 weightedcaptureXS = this->GetDopplerBroa
236 }
237 G4double ela1 = theElasticData->at(indexE
238 G4double ela2 = theElasticData->at(indexE)
239 G4double cap1 = theCaptureData->at(indexE
240 G4double cap2 = theCaptureData->at(indexE)
241 G4double elasticXS = weightedelasticXS / f
242 G4double captureXS = weightedcaptureXS / f
243 xsela_cache[id] = theInt.Lin( kineticEnerg
244 xscap_cache[id] = theInt.Lin( kineticEnerg
245 if ( Z < 88 ) {
246 xsfiss_cache[id] = 0.0;
247 } else {
248 if ( G4ParticleHPManager::GetInstance()-
249 G4double weightedfissionXS = (*G4Parti
250 G4double fissionXS = weightedfissionXS
251 G4double fiss1 = theFissionData->at(in
252 G4double fiss2 = theFissionData->at(in
253 xsfiss_cache[id] = theInt.Lin( kinetic
254 } else {
255 G4double weightedfissionXS = this->Get
256 G4double fissionXS = weightedfissionXS
257 G4double fiss1 = theFissionData->at(in
258 G4double fiss2 = theFissionData->at(in
259 xsfiss_cache[id] = theInt.Lin( kinetic
260 }
261 }
262 }
263 if ( MTnumber == 2 ) { // elastic
264 return xsela_cache[id];
265 } else if ( MTnumber == 102 ) { // radiativ
266 return xscap_cache[id];
267 } else if ( MTnumber == 18 ) { // fission
268 return xsfiss_cache[id];
269 } else {
270 G4cout << "Reaction was not found, returns
271 return 0;
272 }
273 }
274
275 ///-------------------------------------------
276 G4double G4ParticleHPIsoProbabilityTable_NJOY:
277 const G4Element* ele, G4double& kineticEnerg
278 energy_cache[id] = kineticEnergy;
279 if ( kineticEnergy < Emin || kineticEnergy
280 // if the kinetic energy is outside of the
281 G4int indexEl = (G4int)ele->GetIndex();
282 G4int isotopeJ = -1; // index of isotope
283 for ( G4int j = 0; j < (G4int)ele->GetNumb
284 if ( A == (G4int)ele->GetIsotope(j)->Get
285 isotopeJ = j;
286 break;
287 }
288 }
289 G4double frac = ele->GetRelativeAbundanceV
290 G4double weightedelasticXS;
291 G4double weightedcaptureXS;
292 if ( G4ParticleHPManager::GetInstance()->G
293 weightedelasticXS = (*G4ParticleHPManage
294 weightedcaptureXS = (*G4ParticleHPManage
295 } else {
296 weightedelasticXS = this->GetDopplerBroa
297 weightedcaptureXS = this->GetDopplerBroa
298 }
299 xsela_cache[id] = weightedelasticXS / frac
300 xscap_cache[id] = weightedcaptureXS / frac
301 if ( Z < 88 ) {
302 xsfiss_cache[id] = 0.0;
303 } else {
304 if ( G4ParticleHPManager::GetInstance()-
305 G4double weightedfissionXS = (*G4Parti
306 xsfiss_cache[id] = weightedfissionXS /
307 } else {
308 G4double weightedfissionXS = this->Get
309 xsfiss_cache[id] = weightedfissionXS /
310 }
311 }
312 } else if ( lssf_flag == 0 ) {
313 G4int indexE = theEnergies->GetEnergyIndex
314 std::vector< G4double >* theProbability1 =
315 std::vector< G4double >* theProbability2 =
316 G4double ene1 = theEnergies->GetEnergy(ind
317 G4double ene2 = theEnergies->GetEnergy(ind
318 G4double rand = G4UniformRand();
319 random_number_cache[id] = rand;
320 G4int indexP1;
321 G4int indexP2;
322 for ( indexP1 = 0; indexP1 < tableOrder; i
323 if ( rand <= theProbability1->at(indexP1
324 }
325 for ( indexP2 = 0; indexP2 < tableOrder; i
326 if ( rand <= theProbability2->at(indexP2
327 }
328 G4double ela1 = theElasticData->at(indexE
329 G4double ela2 = theElasticData->at(indexE)
330 G4double cap1 = theCaptureData->at(indexE
331 G4double cap2 = theCaptureData->at(indexE)
332 xsela_cache[id] = theInt.Lin( kineticEnerg
333 xscap_cache[id] = theInt.Lin( kineticEnerg
334 if ( Z < 88 ) {
335 xsfiss_cache[id] = 0.0;
336 } else {
337 G4double fiss1 = theFissionData->at(inde
338 G4double fiss2 = theFissionData->at(inde
339 xsfiss_cache[id] = theInt.Lin( kineticEn
340 }
341 } else if ( lssf_flag == 1 ) {
342 G4int indexE = theEnergies->GetEnergyIndex
343 std::vector< G4double >* theProbability1 =
344 std::vector< G4double >* theProbability2 =
345 G4double ene1 = theEnergies->GetEnergy(ind
346 G4double ene2 = theEnergies->GetEnergy(ind
347 G4double rand = G4UniformRand();
348 random_number_cache[id] = rand;
349 G4int indexP1;
350 G4int indexP2;
351 for ( indexP1 = 0; indexP1 < tableOrder; i
352 if ( rand <= theProbability1->at(indexP1
353 }
354 for ( indexP2 = 0; indexP2 < tableOrder; i
355 if ( rand <= theProbability2->at(indexP2
356 }
357 G4int indexEl = (G4int)ele->GetIndex();
358 G4int isotopeJ = -1; // index of isotope
359 for ( G4int j = 0; j < (G4int)ele->GetNumb
360 if ( A == (G4int)ele->GetIsotope(j)->Get
361 isotopeJ = j;
362 break;
363 }
364 }
365 G4double frac = ele->GetRelativeAbundanceV
366 G4double weightedelasticXS;
367 G4double weightedcaptureXS;
368 if ( G4ParticleHPManager::GetInstance()->G
369 weightedelasticXS = (*G4ParticleHPManage
370 weightedcaptureXS = (*G4ParticleHPManage
371 } else {
372 weightedelasticXS = this->GetDopplerBroa
373 weightedcaptureXS = this->GetDopplerBroa
374 }
375 G4double ela1 = theElasticData->at(indexE
376 G4double ela2 = theElasticData->at(indexE)
377 G4double cap1 = theCaptureData->at(indexE
378 G4double cap2 = theCaptureData->at(indexE)
379 G4double elasticXS = weightedelasticXS / f
380 G4double captureXS = weightedcaptureXS / f
381 xsela_cache[id] = theInt.Lin( kineticEnerg
382 xscap_cache[id] = theInt.Lin( kineticEnerg
383 if ( Z < 88 ) {
384 xsfiss_cache[id] = 0.0;
385 } else {
386 if ( G4ParticleHPManager::GetInstance()-
387 G4double weightedfissionXS = (*G4Par
388 G4double fissionXS = weightedfission
389 G4double fiss1 = theFissionData->at(
390 G4double fiss2 = theFissionData->at(
391 xsfiss_cache[id] = theInt.Lin(kineti
392 } else {
393 G4double weightedfissionXS = this->G
394 G4double fissionXS = weightedfission
395 G4double fiss1 = theFissionData->at(
396 G4double fiss2 = theFissionData->at(
397 xsfiss_cache[id] = theInt.Lin( kinet
398 }
399 }
400 }
401 if ( MTnumber == 2 ) { // elastic
402 return xsela_cache[id];
403 } else if ( MTnumber == 102 ) { // radiativ
404 return xscap_cache[id];
405 } else if ( MTnumber == 18 ) { // fission
406 return xsfiss_cache[id];
407 } else {
408 G4cout << "Reaction was not found, returns
409 return 0;
410 }
411 }
412