| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // Geant4 source file
30 //
31 // File name: G4ParticleHPIsoProbabilityTable.cc
32 //
33 // Authors: Marek Zmeskal (CTU, Czech Technical University in Prague, Czech Republic)
34 // Loic Thulliez (CEA France)
35 //
36 // Creation date: 4 June 2024
37 //
38 // Description: Class for the probability table of the given isotope
39 // and for the given temperature.
40 // It reads the files with probability tables and
41 // finds the correct cross-section.
42 //
43 // Modifications:
44 //
45 // -------------------------------------------------------------------
46 //
47 //
48
49 #include "G4ParticleHPIsoProbabilityTable.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "Randomize.hh"
52 #include "G4Exception.hh"
53 #include "G4NucleiProperties.hh"
54 #include "G4ReactionProduct.hh"
55 #include "G4ParticleHPManager.hh"
56 #include "G4ParticleHPVector.hh"
57 #include "G4DynamicParticle.hh"
58 #include "G4Element.hh"
59 #include "G4Nucleus.hh"
60 #include "G4Neutron.hh"
61 #include "G4ParticleHPChannel.hh"
62 #include "G4ParticleHPChannelList.hh"
63
64 #include <string>
65 #include <fstream>
66
67 ///--------------------------------------------------------------------------------------
68 G4ParticleHPIsoProbabilityTable::~G4ParticleHPIsoProbabilityTable()
69 {
70 for ( auto it = theProbabilities->cbegin(); it != theProbabilities->cend(); ++it ) {
71 delete* it;
72 }
73 for ( auto it = theElasticData->cbegin(); it != theElasticData->cend(); ++it ) {
74 delete* it;
75 }
76 for ( auto it = theCaptureData->cbegin(); it != theCaptureData->cend(); ++it ) {
77 delete* it;
78 }
79 for ( auto it = theFissionData->cbegin(); it != theFissionData->cend(); ++it ) {
80 delete* it;
81 }
82 delete theEnergies;
83 delete theProbabilities;
84 delete theElasticData;
85 delete theCaptureData;
86 delete theFissionData;
87 }
88
89 ///--------------------------------------------------------------------------------------
90 void G4ParticleHPIsoProbabilityTable::Init( G4int, G4int, G4int, G4double, const G4String& ) {}
91
92 ///--------------------------------------------------------------------------------------
93 G4double G4ParticleHPIsoProbabilityTable::GetCorrelatedIsoCrossSectionPT( const G4DynamicParticle*, G4int,
94 const G4Element*, G4double&, G4double&, std::thread::id& )
95 {
96 G4Exception( "G4ParticleHPIsoProbabilityTable::GetCorrelatedIsoCrossSectionPT", "hadhp01",
97 FatalException, "The base method G4ParticleHPIsoProbabilityTable::GetCorrelatedIsoCrossSectionPT"
98 "is trying to be accessed, which is not allowed,"
99 "the derived class for NJOY or CALENDF was not properly declared." );
100 return 0.0;
101 }
102
103 ///--------------------------------------------------------------------------------------
104 G4double G4ParticleHPIsoProbabilityTable::GetIsoCrossSectionPT( const G4DynamicParticle*, G4int,
105 const G4Element*, G4double&, std::map< std::thread::id, G4double >&, std::thread::id& )
106 {
107 G4Exception( "G4ParticleHPIsoProbabilityTable::GetIsoCrossSectionPT", "hadhp01",
108 FatalException, "The base method G4ParticleHPIsoProbabilityTable::GetIsoCrossSectionPT"
109 "is trying to be accessed, which is not allowed,"
110 "the derived class for NJOY or CALENDF was not properly declared." );
111 return 0.0;
112 }
113
114 ///--------------------------------------------------------------------------------------
115 G4double G4ParticleHPIsoProbabilityTable::GetDopplerBroadenedElasticXS( const G4DynamicParticle* dP, G4int indexEl, G4int isotopeJ ) {
116 G4double result = 0.0;
117 G4ReactionProduct theNeutron( dP->GetDefinition() );
118 theNeutron.SetMomentum( dP->GetMomentum() );
119 theNeutron.SetKineticEnergy( dP->GetKineticEnergy() );
120 // prepare thermal nucleus
121 G4Nucleus aNuc;
122 G4double eleMass = G4NucleiProperties::GetNuclearMass(A, Z) / G4Neutron::Neutron()->GetPDGMass();
123 G4ReactionProduct boosted;
124 G4double aXsection;
125 G4int counter = 0;
126 G4double buffer = 0.0;
127 G4int size = G4int( std::max( 10.0, T / 60.0 * kelvin ) );
128 G4ThreeVector neutronVelocity = 1.0 / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
129 G4double neutronVMag = neutronVelocity.mag();
130 while ( counter == 0 || std::abs( buffer - result / std::max(1, counter) ) > 0.03 * buffer ) {
131 if ( counter ) buffer = result / counter;
132 while ( counter < size ) {
133 ++counter;
134 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass, T );
135 boosted.Lorentz( theNeutron, aThermalNuc );
136 G4double theEkin = boosted.GetKineticEnergy();
137 aXsection = (*G4ParticleHPManager::GetInstance()->GetElasticFinalStates())[indexEl]->GetWeightedXsec( theEkin, isotopeJ );
138 // velocity correction
139 G4ThreeVector targetVelocity = 1.0 / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
140 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
141 result += aXsection;
142 }
143 size += size;
144 }
145 result /= counter;
146 return result;
147 }
148
149 ///--------------------------------------------------------------------------------------
150 G4double G4ParticleHPIsoProbabilityTable::GetDopplerBroadenedCaptureXS( const G4DynamicParticle* dP, G4int indexEl, G4int isotopeJ ) {
151 G4double result = 0.0;
152 G4ReactionProduct theNeutron( dP->GetDefinition() );
153 theNeutron.SetMomentum( dP->GetMomentum() );
154 theNeutron.SetKineticEnergy( dP->GetKineticEnergy() );
155 // prepare thermal nucleus
156 G4Nucleus aNuc;
157 G4double eleMass = G4NucleiProperties::GetNuclearMass(A, Z) / G4Neutron::Neutron()->GetPDGMass();
158 G4ReactionProduct boosted;
159 G4double aXsection;
160 G4int counter = 0;
161 G4double buffer = 0.0;
162 G4int size = G4int( std::max( 10.0, T / 60.0 * kelvin));
163 G4ThreeVector neutronVelocity = 1.0 / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
164 G4double neutronVMag = neutronVelocity.mag();
165 while ( counter == 0 || std::abs( buffer - result / std::max( 1, counter ) ) > 0.03 * buffer ) {
166 if ( counter ) buffer = result / counter;
167 while ( counter < size ) {
168 ++counter;
169 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass, T );
170 boosted.Lorentz( theNeutron, aThermalNuc );
171 G4double theEkin = boosted.GetKineticEnergy();
172 aXsection = (*G4ParticleHPManager::GetInstance()->GetCaptureFinalStates())[indexEl]->GetWeightedXsec( theEkin, isotopeJ );
173 // velocity correction
174 G4ThreeVector targetVelocity = 1.0 / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
175 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
176 result += aXsection;
177 }
178 size += size;
179 }
180 result /= counter;
181 return result;
182 }
183
184 ///--------------------------------------------------------------------------------------
185 G4double G4ParticleHPIsoProbabilityTable::GetDopplerBroadenedFissionXS( const G4DynamicParticle* dP, G4int indexEl, G4int isotopeJ ) {
186 G4double result = 0.0;
187 G4ReactionProduct theNeutron( dP->GetDefinition() );
188 theNeutron.SetMomentum( dP->GetMomentum() );
189 theNeutron.SetKineticEnergy( dP->GetKineticEnergy() );
190 // prepare thermal nucleus
191 G4Nucleus aNuc;
192 G4double eleMass;
193 eleMass = G4NucleiProperties::GetNuclearMass(A, Z) / G4Neutron::Neutron()->GetPDGMass();
194 G4ReactionProduct boosted;
195 G4double aXsection;
196 G4int counter = 0;
197 G4double buffer = 0.0;
198 G4int size = G4int( std::max( 10.0, T / 60.0 * kelvin ) );
199 G4ThreeVector neutronVelocity = 1.0 / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
200 G4double neutronVMag = neutronVelocity.mag();
201 while ( counter == 0 || std::abs( buffer - result / std::max( 1, counter ) ) > 0.01 * buffer ) {
202 if ( counter ) buffer = result / counter;
203 while ( counter < size ) {
204 ++counter;
205 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass, T );
206 boosted.Lorentz( theNeutron, aThermalNuc );
207 G4double theEkin = boosted.GetKineticEnergy();
208 aXsection = (*G4ParticleHPManager::GetInstance()->GetFissionFinalStates())[indexEl]->GetWeightedXsec( theEkin, isotopeJ );
209 // velocity correction
210 G4ThreeVector targetVelocity = 1.0 / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
211 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
212 result += aXsection;
213 }
214 size += size;
215 }
216 result /= counter;
217 return result;
218 }
219
220 ///--------------------------------------------------------------------------------------
221 G4double G4ParticleHPIsoProbabilityTable::GetDopplerBroadenedInelasticXS( const G4DynamicParticle* dP, G4int indexEl, G4int isotopeJ ) {
222 G4double result = 0.0;
223 G4ReactionProduct theNeutron( dP->GetDefinition() );
224 theNeutron.SetMomentum( dP->GetMomentum() );
225 theNeutron.SetKineticEnergy( dP->GetKineticEnergy() );
226 // prepare thermal nucleus
227 G4Nucleus aNuc;
228 G4double eleMass = G4NucleiProperties::GetNuclearMass(A, Z) / G4Neutron::Neutron()->GetPDGMass();
229 G4ReactionProduct boosted;
230 G4double aXsection;
231 G4int counter = 0;
232 G4int failCount = 0;
233 G4double buffer = 0.0;
234 G4int size = G4int( std::max( 10.0, T / 60.0 * kelvin ) );
235 G4ThreeVector neutronVelocity = 1.0 / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
236 G4double neutronVMag = neutronVelocity.mag();
237 #ifndef G4PHP_DOPPLER_LOOP_ONCE
238 while ( counter == 0 || std::abs( buffer - result / std::max( 1, counter ) ) > 0.01 * buffer ) {
239 if ( counter ) buffer = result / counter;
240 while ( counter < size ) {
241 ++counter;
242 #endif
243 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass / G4Neutron::Neutron()->GetPDGMass(), T );
244 boosted.Lorentz( theNeutron, aThermalNuc );
245 G4double theEkin = boosted.GetKineticEnergy();
246 aXsection = ((*G4ParticleHPManager::GetInstance()->GetInelasticFinalStates( dP->GetDefinition() ))[indexEl]
247 ->GetWeightedXsec( theEkin, isotopeJ ) ) * barn;
248 if ( aXsection < 0.0 ) {
249 if ( failCount < 1000 ) {
250 ++failCount;
251 #ifndef G4PHP_DOPPLER_LOOP_ONCE
252 --counter;
253 continue;
254 #endif
255 } else {
256 aXsection = 0.0;
257 }
258 }
259 // velocity correction
260 G4ThreeVector targetVelocity = 1.0 / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
261 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
262 result += aXsection;
263 }
264 #ifndef G4PHP_DOPPLER_LOOP_ONCE
265 size += size;
266 }
267 result /= counter;
268 #endif
269 return result;
270 }
271