| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 070618 fix memory leaking by T. Koi 30 // 070618 fix memory leaking by T. Koi
31 // 071002 enable cross section dump by T. Koi 31 // 071002 enable cross section dump by T. Koi
32 // 081024 G4NucleiPropertiesTable:: to G4Nucle 32 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33 // 081124 Protect invalid read which caused ru 33 // 081124 Protect invalid read which caused run time errors by T. Koi
34 // 100729 Add safty for 0 lenght cross section 34 // 100729 Add safty for 0 lenght cross sections by T. Ko
35 // P. Arce, June-2014 Conversion neutron_hp to 35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 // 36 //
37 #include "G4ParticleHPFissionData.hh" 37 #include "G4ParticleHPFissionData.hh"
38 << 38 #include "G4ParticleHPManager.hh"
39 #include "G4ElementTable.hh" << 39 #include "G4PhysicalConstants.hh"
40 #include "G4HadronicParameters.hh" << 40 #include "G4SystemOfUnits.hh"
41 #include "G4Neutron.hh" 41 #include "G4Neutron.hh"
42 #include "G4NucleiProperties.hh" << 42 #include "G4ElementTable.hh"
43 #include "G4ParticleHPData.hh" 43 #include "G4ParticleHPData.hh"
44 #include "G4ParticleHPManager.hh" 44 #include "G4ParticleHPManager.hh"
45 #include "G4PhysicalConstants.hh" << 45 #include "G4HadronicParameters.hh"
46 #include "G4Pow.hh" 46 #include "G4Pow.hh"
47 #include "G4SystemOfUnits.hh" <<
48 47
49 G4ParticleHPFissionData::G4ParticleHPFissionDa << 48 G4ParticleHPFissionData::G4ParticleHPFissionData()
>> 49 :G4VCrossSectionDataSet("NeutronHPFissionXS")
50 { 50 {
51 SetMinKinEnergy(0 * MeV); << 51 SetMinKinEnergy( 0*MeV );
52 SetMaxKinEnergy(20 * MeV); << 52 SetMaxKinEnergy( 20*MeV );
53 53
54 theCrossSections = nullptr; << 54 theCrossSections = 0;
55 instanceOfWorker = false; << 55 onFlightDB = true;
56 if (G4Threading::IsWorkerThread()) { << 56 instanceOfWorker = false;
57 instanceOfWorker = true; << 57 if ( G4Threading::IsWorkerThread() ) {
58 } << 58 instanceOfWorker = true;
59 element_cache = nullptr; << 59 }
60 material_cache = nullptr; << 60 element_cache = NULL;
61 ke_cache = 0.0; << 61 material_cache = NULL;
62 xs_cache = 0.0; << 62 ke_cache = 0.0;
>> 63 xs_cache = 0.0;
>> 64 //BuildPhysicsTable(*G4Neutron::Neutron());
63 } 65 }
64 << 66
65 G4ParticleHPFissionData::~G4ParticleHPFissionD 67 G4ParticleHPFissionData::~G4ParticleHPFissionData()
66 { 68 {
67 if (theCrossSections != nullptr && !instance << 69 if ( theCrossSections != NULL && instanceOfWorker != true ) {
68 theCrossSections->clearAndDestroy(); << 70 theCrossSections->clearAndDestroy();
69 delete theCrossSections; << 71 delete theCrossSections;
70 theCrossSections = nullptr; << 72 theCrossSections = NULL;
71 } << 73 }
72 } 74 }
73 75
74 G4bool G4ParticleHPFissionData::IsIsoApplicabl << 76 G4bool G4ParticleHPFissionData::IsIsoApplicable( const G4DynamicParticle* dp ,
75 << 77 G4int /*Z*/ , G4int /*A*/ ,
76 << 78 const G4Element* /*elm*/ ,
>> 79 const G4Material* /*mat*/ )
77 { 80 {
78 G4double eKin = dp->GetKineticEnergy(); << 81 G4double eKin = dp->GetKineticEnergy();
79 return eKin <= GetMaxKinEnergy() && eKin >= << 82 if ( eKin > GetMaxKinEnergy()
80 && dp->GetDefinition() == G4Neutron:: << 83 || eKin < GetMinKinEnergy()
>> 84 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
>> 85
>> 86 return true;
81 } 87 }
82 88
83 G4double G4ParticleHPFissionData::GetIsoCrossS << 89 G4double G4ParticleHPFissionData::GetIsoCrossSection( const G4DynamicParticle* dp ,
84 << 90 G4int /*Z*/ , G4int /*A*/ ,
85 << 91 const G4Isotope* /*iso*/ ,
86 << 92 const G4Element* element ,
>> 93 const G4Material* material )
87 { 94 {
88 if (dp->GetKineticEnergy() == ke_cache && el << 95 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
89 return xs_cache; <<
90 96
91 ke_cache = dp->GetKineticEnergy(); << 97 ke_cache = dp->GetKineticEnergy();
92 element_cache = element; << 98 element_cache = element;
93 material_cache = material; << 99 material_cache = material;
94 G4double xs = GetCrossSection(dp, element, m << 100 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
95 xs_cache = xs; << 101 xs_cache = xs;
96 return xs; << 102 return xs;
97 } 103 }
98 104
99 void G4ParticleHPFissionData::BuildPhysicsTabl << 105 /*
>> 106 G4bool G4ParticleHPFissionData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
100 { 107 {
101 if (G4Threading::IsWorkerThread()) { << 108 G4bool result = true;
102 theCrossSections = G4ParticleHPManager::Ge << 109 G4double eKin = aP->GetKineticEnergy();
103 return; << 110 if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
104 } << 111 return result;
>> 112 }
>> 113 */
>> 114
>> 115 void G4ParticleHPFissionData::BuildPhysicsTable(const G4ParticleDefinition& aP)
>> 116 {
>> 117
>> 118 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() ) {
>> 119 onFlightDB = false;
>> 120 #ifdef G4VERBOSE
>> 121 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) {
>> 122 G4cout << "Find a flag of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
>> 123 G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of fission reaction of neutrons (<20MeV)." << G4endl;
>> 124 }
>> 125 #endif
>> 126 }
105 127
106 std::size_t numberOfElements = G4Element::Ge << 128 if(&aP!=G4Neutron::Neutron())
107 if (theCrossSections == nullptr) << 129 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
108 theCrossSections = new G4PhysicsTable(numb << 130
109 else << 131 if ( G4Threading::IsWorkerThread() ) {
110 theCrossSections->clearAndDestroy(); << 132 theCrossSections = G4ParticleHPManager::GetInstance()->GetFissionCrossSections();
>> 133 return;
>> 134 }
>> 135
>> 136 size_t numberOfElements = G4Element::GetNumberOfElements();
>> 137 //theCrossSections = new G4PhysicsTable( numberOfElements );
>> 138 // TKDB
>> 139 //if ( theCrossSections == NULL ) theCrossSections = new G4PhysicsTable( numberOfElements );
>> 140 if ( theCrossSections == NULL )
>> 141 theCrossSections = new G4PhysicsTable( numberOfElements );
>> 142 else
>> 143 theCrossSections->clearAndDestroy();
111 144
112 // make a PhysicsVector for each element 145 // make a PhysicsVector for each element
113 146
114 auto theElementTable = G4Element::GetElement << 147 static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
115 for (std::size_t i = 0; i < numberOfElements << 148 for( size_t i=0; i<numberOfElements; ++i )
116 G4PhysicsVector* physVec = G4ParticleHPDat << 149 {
117 ->MakePhysics << 150 G4PhysicsVector* physVec = G4ParticleHPData::
>> 151 Instance(G4Neutron::Neutron())->MakePhysicsVector((*theElementTable)[i], this);
118 theCrossSections->push_back(physVec); 152 theCrossSections->push_back(physVec);
119 } 153 }
120 154
121 G4ParticleHPManager::GetInstance()->Register << 155 G4ParticleHPManager::GetInstance()->RegisterFissionCrossSections( theCrossSections );
122 } 156 }
123 157
124 void G4ParticleHPFissionData::DumpPhysicsTable << 158 void G4ParticleHPFissionData::DumpPhysicsTable(const G4ParticleDefinition& aP)
125 { 159 {
126 #ifdef G4VERBOSE << 160 if(&aP!=G4Neutron::Neutron())
127 if (G4HadronicParameters::Instance()->GetVer << 161 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
>> 162
>> 163 #ifdef G4VERBOSE
>> 164 if ( G4HadronicParameters::Instance()->GetVerboseLevel() == 0 ) return;
>> 165
>> 166 //
>> 167 // Dump element based cross section
>> 168 // range 10e-5 eV to 20 MeV
>> 169 // 10 point per decade
>> 170 // in barn
>> 171 //
>> 172
>> 173 G4cout << G4endl;
>> 174 G4cout << G4endl;
>> 175 G4cout << "Fission Cross Section of Neutron HP"<< G4endl;
>> 176 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
>> 177 G4cout << G4endl;
>> 178 G4cout << "Name of Element" << G4endl;
>> 179 G4cout << "Energy[eV] XS[barn]" << G4endl;
>> 180 G4cout << G4endl;
>> 181
>> 182 size_t numberOfElements = G4Element::GetNumberOfElements();
>> 183 static G4ThreadLocal G4ElementTable *theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
>> 184
>> 185 for ( size_t i = 0 ; i < numberOfElements ; ++i )
>> 186 {
>> 187
>> 188 G4cout << (*theElementTable)[i]->GetName() << G4endl;
>> 189
>> 190 if ( (*((*theCrossSections)(i))).GetVectorLength() == 0 )
>> 191 {
>> 192 G4cout << "The cross-section data of the fission of this element is not available." << G4endl;
>> 193 G4cout << G4endl;
>> 194 continue;
>> 195 }
128 196
129 // << 197 G4int ie = 0;
130 // Dump element based cross section <<
131 // range 10e-5 eV to 20 MeV <<
132 // 10 point per decade <<
133 // in barn <<
134 // <<
135 G4cout << G4endl; <<
136 G4cout << G4endl; <<
137 G4cout << "Fission Cross Section of Neutron <<
138 G4cout << "(Pointwise cross-section at 0 Kel <<
139 G4cout << G4endl; <<
140 G4cout << "Name of Element" << G4endl; <<
141 G4cout << "Energy[eV] XS[barn]" << G4endl; <<
142 G4cout << G4endl; <<
143 198
144 std::size_t numberOfElements = G4Element::Ge << 199 for ( ie = 0 ; ie < 130 ; ie++ )
145 auto theElementTable = G4Element::GetElement << 200 {
>> 201 G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *eV;
>> 202 G4bool outOfRange = false;
>> 203
>> 204 if ( eKinetic < 20*MeV )
>> 205 {
>> 206 G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
>> 207 }
146 208
147 for (std::size_t i = 0; i < numberOfElements << 209 }
148 G4cout << (*theElementTable)[i]->GetName() <<
149 210
150 if ((*((*theCrossSections)(i))).GetVectorL <<
151 G4cout << "The cross-section data of the <<
152 G4cout << G4endl; 211 G4cout << G4endl;
153 continue; << 212 }
154 } <<
155 213
156 for (G4int ie = 0; ie < 130; ++ie) { << 214 //G4cout << "G4ParticleHPFissionData::DumpPhysicsTable still to be implemented"<<G4endl;
157 G4double eKinetic = 1.0e-5 * G4Pow::GetI << 215 #endif
158 G4bool outOfRange = false; <<
159 <<
160 if (eKinetic < 20 * MeV) { <<
161 G4cout << eKinetic / eV << " " <<
162 << (*((*theCrossSections)(i))). <<
163 } <<
164 } <<
165 <<
166 G4cout << G4endl; <<
167 } <<
168 #endif <<
169 } 216 }
170 217
171 G4double G4ParticleHPFissionData::GetCrossSect << 218 #include "G4NucleiProperties.hh"
172 << 219
>> 220 G4double G4ParticleHPFissionData::
>> 221 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
173 { 222 {
174 G4double result = 0; 223 G4double result = 0;
175 if (anE->GetZ() < 88) return result; << 224 if(anE->GetZ()<88) return result;
176 G4bool outOfRange; 225 G4bool outOfRange;
177 auto index = (G4int)anE->GetIndex(); << 226 G4int index = anE->GetIndex();
178 227
179 if (((*theCrossSections)(index))->GetVectorL << 228 // 100729 TK add safety
>> 229 if ( ( ( *theCrossSections )( index ) )->GetVectorLength() == 0 ) return result;
180 230
181 // prepare neutron 231 // prepare neutron
182 G4double eKinetic = aP->GetKineticEnergy(); 232 G4double eKinetic = aP->GetKineticEnergy();
183 G4ReactionProduct theNeutronRP(aP->GetDefini << 233 G4ReactionProduct theNeutronRP( aP->GetDefinition() );
184 theNeutronRP.SetMomentum(aP->GetMomentum()); << 234 theNeutronRP.SetMomentum( aP->GetMomentum() );
185 theNeutronRP.SetKineticEnergy(eKinetic); << 235 theNeutronRP.SetKineticEnergy( eKinetic );
186 << 236
187 if (G4ParticleHPManager::GetInstance()->GetN << 237 if ( !onFlightDB ) {
188 // NEGLECT_DOPPLER << 238 //NEGLECT_DOPPLER
189 G4double factor = 1.0; << 239 G4double factor = 1.0;
190 if (eKinetic < aT * k_Boltzmann) { << 240 if ( eKinetic < aT * k_Boltzmann ) {
191 // below 0.1 eV neutrons << 241 // below 0.1 eV neutrons
192 // Have to do some, but now just igonre. << 242 // Have to do some, but now just igonre.
193 // Will take care after performance chec << 243 // Will take care after performance check.
194 // factor = factor * targetV; << 244 // factor = factor * targetV;
195 } << 245 }
196 return ((*((*theCrossSections)(index))).Ge << 246 return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
197 } 247 }
198 248
199 // prepare thermal nucleus 249 // prepare thermal nucleus
200 G4Nucleus aNuc; 250 G4Nucleus aNuc;
201 G4double eps = 0.0001; 251 G4double eps = 0.0001;
202 G4double theA = anE->GetN(); 252 G4double theA = anE->GetN();
203 G4double theZ = anE->GetZ(); 253 G4double theZ = anE->GetZ();
204 G4double eleMass; << 254 G4double eleMass;
205 eleMass = (G4NucleiProperties::GetNuclearMas << 255 eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
206 << 256 ) / G4Neutron::Neutron()->GetPDGMass();
207 / G4Neutron::Neutron()->GetPDGMass << 257
208 <<
209 G4ReactionProduct boosted; 258 G4ReactionProduct boosted;
210 G4double aXsection; 259 G4double aXsection;
211 << 260
212 // MC integration loop 261 // MC integration loop
213 G4int counter = 0; 262 G4int counter = 0;
214 G4double buffer = 0; 263 G4double buffer = 0;
215 G4int size = G4int(std::max(10., aT / 60 * k << 264 G4int size = G4int(std::max(10., aT/60*kelvin));
216 G4ThreeVector neutronVelocity = << 265 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutronRP.GetMomentum();
217 1. / G4Neutron::Neutron()->GetPDGMass() * <<
218 G4double neutronVMag = neutronVelocity.mag() 266 G4double neutronVMag = neutronVelocity.mag();
219 267
220 while (counter == 0 << 268 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer) // Loop checking, 11.05.2015, T. Koi
221 || std::abs(buffer - result / std::ma <<
222 > 0.01 * buffer) // Loop checki <<
223 { 269 {
224 if (counter != 0) buffer = result / counte << 270 if(counter) buffer = result/counter;
225 while (counter < size) // Loop checking, << 271 while (counter<size) // Loop checking, 11.05.2015, T. Koi
226 { 272 {
227 counter++; << 273 counter ++;
228 G4ReactionProduct aThermalNuc = aNuc.Get 274 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
229 boosted.Lorentz(theNeutronRP, aThermalNu 275 boosted.Lorentz(theNeutronRP, aThermalNuc);
230 G4double theEkin = boosted.GetKineticEne 276 G4double theEkin = boosted.GetKineticEnergy();
231 aXsection = (*((*theCrossSections)(index 277 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
232 // velocity correction. 278 // velocity correction.
233 G4ThreeVector targetVelocity = 1. / aThe << 279 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
234 aXsection *= (targetVelocity - neutronVe << 280 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
235 result += aXsection; 281 result += aXsection;
236 } 282 }
237 size += size; 283 size += size;
238 } 284 }
239 result /= counter; 285 result /= counter;
240 return result; 286 return result;
241 } 287 }
242 288
243 G4int G4ParticleHPFissionData::GetVerboseLevel << 289 G4int G4ParticleHPFissionData::GetVerboseLevel() const
244 { 290 {
245 return G4ParticleHPManager::GetInstance()->G << 291 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
246 } 292 }
247 << 293 void G4ParticleHPFissionData::SetVerboseLevel( G4int newValue )
248 void G4ParticleHPFissionData::SetVerboseLevel( <<
249 { 294 {
250 G4ParticleHPManager::GetInstance()->SetVerbo << 295 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
251 } 296 }
252 <<
253 void G4ParticleHPFissionData::CrossSectionDesc 297 void G4ParticleHPFissionData::CrossSectionDescription(std::ostream& outFile) const
254 { 298 {
255 outFile << "High Precision cross data based << 299 outFile << "High Precision cross data based on Evaluated Nuclear Data Files (ENDF) for induced fission reaction of neutrons below 20MeV\n" ;
256 << "for induced fission reaction of <<
257 } 300 }
258 301