| 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 "G4NucleiProperties.hh"
46 #include "G4Pow.hh" 47 #include "G4Pow.hh"
47 #include "G4SystemOfUnits.hh" <<
48 48
49 G4ParticleHPFissionData::G4ParticleHPFissionDa << 49 G4ParticleHPFissionData::G4ParticleHPFissionData()
>> 50 :G4VCrossSectionDataSet("NeutronHPFissionXS")
50 { 51 {
51 SetMinKinEnergy(0 * MeV); << 52 SetMinKinEnergy( 0*MeV );
52 SetMaxKinEnergy(20 * MeV); << 53 SetMaxKinEnergy( 20*MeV );
53 54
54 theCrossSections = nullptr; << 55 theCrossSections = 0;
55 instanceOfWorker = false; << 56 instanceOfWorker = false;
56 if (G4Threading::IsWorkerThread()) { << 57 if ( G4Threading::IsWorkerThread() ) {
57 instanceOfWorker = true; << 58 instanceOfWorker = true;
58 } << 59 }
59 element_cache = nullptr; << 60 element_cache = nullptr;
60 material_cache = nullptr; << 61 material_cache = nullptr;
61 ke_cache = 0.0; << 62 ke_cache = 0.0;
62 xs_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 != nullptr && instanceOfWorker != true ) {
68 theCrossSections->clearAndDestroy(); << 70 theCrossSections->clearAndDestroy();
69 delete theCrossSections; << 71 delete theCrossSections;
70 theCrossSections = nullptr; << 72 theCrossSections = nullptr;
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 void G4ParticleHPFissionData::BuildPhysicsTable(const G4ParticleDefinition& aP)
100 { 106 {
101 if (G4Threading::IsWorkerThread()) { << 107 if(&aP!=G4Neutron::Neutron())
102 theCrossSections = G4ParticleHPManager::Ge << 108 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
103 return; << 109
104 } << 110 if ( G4Threading::IsWorkerThread() ) {
>> 111 theCrossSections = G4ParticleHPManager::GetInstance()->GetFissionCrossSections();
>> 112 return;
>> 113 }
105 114
106 std::size_t numberOfElements = G4Element::Ge 115 std::size_t numberOfElements = G4Element::GetNumberOfElements();
107 if (theCrossSections == nullptr) << 116 if ( theCrossSections == nullptr )
108 theCrossSections = new G4PhysicsTable(numb << 117 theCrossSections = new G4PhysicsTable( numberOfElements );
109 else 118 else
110 theCrossSections->clearAndDestroy(); << 119 theCrossSections->clearAndDestroy();
111 120
112 // make a PhysicsVector for each element 121 // make a PhysicsVector for each element
113 122
114 auto theElementTable = G4Element::GetElement << 123 static G4ThreadLocal G4ElementTable *theElementTable = nullptr ;
115 for (std::size_t i = 0; i < numberOfElements << 124 if (!theElementTable) theElementTable= G4Element::GetElementTable();
116 G4PhysicsVector* physVec = G4ParticleHPDat << 125 for( std::size_t i=0; i<numberOfElements; ++i )
117 ->MakePhysics << 126 {
>> 127 G4PhysicsVector* physVec = G4ParticleHPData::
>> 128 Instance(G4Neutron::Neutron())->MakePhysicsVector((*theElementTable)[i], this);
118 theCrossSections->push_back(physVec); 129 theCrossSections->push_back(physVec);
119 } 130 }
120 131
121 G4ParticleHPManager::GetInstance()->Register << 132 G4ParticleHPManager::GetInstance()->RegisterFissionCrossSections( theCrossSections );
122 } 133 }
123 134
124 void G4ParticleHPFissionData::DumpPhysicsTable << 135 void G4ParticleHPFissionData::DumpPhysicsTable(const G4ParticleDefinition& aP)
125 { 136 {
>> 137 if(&aP!=G4Neutron::Neutron())
>> 138 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
>> 139
126 #ifdef G4VERBOSE 140 #ifdef G4VERBOSE
127 if (G4HadronicParameters::Instance()->GetVer << 141 if ( G4HadronicParameters::Instance()->GetVerboseLevel() == 0 ) return;
128 << 142
129 // 143 //
130 // Dump element based cross section 144 // Dump element based cross section
131 // range 10e-5 eV to 20 MeV 145 // range 10e-5 eV to 20 MeV
132 // 10 point per decade 146 // 10 point per decade
133 // in barn 147 // in barn
134 // 148 //
135 G4cout << G4endl; << 149 G4cout << G4endl;
136 G4cout << G4endl; << 150 G4cout << G4endl;
137 G4cout << "Fission Cross Section of Neutron << 151 G4cout << "Fission Cross Section of Neutron HP"<< G4endl;
138 G4cout << "(Pointwise cross-section at 0 Kel << 152 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
139 G4cout << G4endl; << 153 G4cout << G4endl;
140 G4cout << "Name of Element" << G4endl; << 154 G4cout << "Name of Element" << G4endl;
141 G4cout << "Energy[eV] XS[barn]" << G4endl; << 155 G4cout << "Energy[eV] XS[barn]" << G4endl;
142 G4cout << G4endl; << 156 G4cout << G4endl;
143 << 157
144 std::size_t numberOfElements = G4Element::Ge << 158 std::size_t numberOfElements = G4Element::GetNumberOfElements();
145 auto theElementTable = G4Element::GetElement << 159 static G4ThreadLocal G4ElementTable *theElementTable = nullptr ;
146 << 160 if (!theElementTable) theElementTable= G4Element::GetElementTable();
147 for (std::size_t i = 0; i < numberOfElements << 161
148 G4cout << (*theElementTable)[i]->GetName() << 162 for ( std::size_t i = 0 ; i < numberOfElements ; ++i )
149 << 163 {
150 if ((*((*theCrossSections)(i))).GetVectorL << 164 G4cout << (*theElementTable)[i]->GetName() << G4endl;
151 G4cout << "The cross-section data of the << 165
152 G4cout << G4endl; << 166 if ( (*((*theCrossSections)(i))).GetVectorLength() == 0 )
153 continue; << 167 {
154 } << 168 G4cout << "The cross-section data of the fission of this element is not available." << G4endl;
>> 169 G4cout << G4endl;
>> 170 continue;
>> 171 }
155 172
156 for (G4int ie = 0; ie < 130; ++ie) { << 173 for ( G4int ie = 0 ; ie < 130 ; ++ie )
157 G4double eKinetic = 1.0e-5 * G4Pow::GetI << 174 {
158 G4bool outOfRange = false; << 175 G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *eV;
159 << 176 G4bool outOfRange = false;
160 if (eKinetic < 20 * MeV) { << 177
161 G4cout << eKinetic / eV << " " << 178 if ( eKinetic < 20*MeV )
162 << (*((*theCrossSections)(i))). << 179 {
>> 180 G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
>> 181 }
163 } 182 }
164 } <<
165 183
166 G4cout << G4endl; << 184 G4cout << G4endl;
167 } << 185 }
168 #endif 186 #endif
169 } 187 }
170 188
171 G4double G4ParticleHPFissionData::GetCrossSect << 189 G4double G4ParticleHPFissionData::
172 << 190 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
173 { 191 {
174 G4double result = 0; 192 G4double result = 0;
175 if (anE->GetZ() < 88) return result; << 193 if(anE->GetZ()<88) return result;
176 G4bool outOfRange; 194 G4bool outOfRange;
177 auto index = (G4int)anE->GetIndex(); << 195 G4int index = (G4int)anE->GetIndex();
178 196
179 if (((*theCrossSections)(index))->GetVectorL << 197 if ( ( ( *theCrossSections )( index ) )->GetVectorLength() == 0 )
>> 198 return result;
180 199
181 // prepare neutron 200 // prepare neutron
182 G4double eKinetic = aP->GetKineticEnergy(); 201 G4double eKinetic = aP->GetKineticEnergy();
183 G4ReactionProduct theNeutronRP(aP->GetDefini << 202 G4ReactionProduct theNeutronRP( aP->GetDefinition() );
184 theNeutronRP.SetMomentum(aP->GetMomentum()); << 203 theNeutronRP.SetMomentum( aP->GetMomentum() );
185 theNeutronRP.SetKineticEnergy(eKinetic); << 204 theNeutronRP.SetKineticEnergy( eKinetic );
186 << 205
187 if (G4ParticleHPManager::GetInstance()->GetN << 206 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() )
188 // NEGLECT_DOPPLER << 207 {
189 G4double factor = 1.0; << 208 //NEGLECT_DOPPLER
190 if (eKinetic < aT * k_Boltzmann) { << 209 G4double factor = 1.0;
191 // below 0.1 eV neutrons << 210 if ( eKinetic < aT * k_Boltzmann ) {
192 // Have to do some, but now just igonre. << 211 // below 0.1 eV neutrons
193 // Will take care after performance chec << 212 // Have to do some, but now just igonre.
194 // factor = factor * targetV; << 213 // Will take care after performance check.
195 } << 214 // factor = factor * targetV;
196 return ((*((*theCrossSections)(index))).Ge << 215 }
>> 216 return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
197 } 217 }
198 218
199 // prepare thermal nucleus 219 // prepare thermal nucleus
200 G4Nucleus aNuc; 220 G4Nucleus aNuc;
201 G4double eps = 0.0001; 221 G4double eps = 0.0001;
202 G4double theA = anE->GetN(); 222 G4double theA = anE->GetN();
203 G4double theZ = anE->GetZ(); 223 G4double theZ = anE->GetZ();
204 G4double eleMass; << 224 G4double eleMass;
205 eleMass = (G4NucleiProperties::GetNuclearMas << 225 eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
206 << 226 ) / G4Neutron::Neutron()->GetPDGMass();
207 / G4Neutron::Neutron()->GetPDGMass << 227
208 <<
209 G4ReactionProduct boosted; 228 G4ReactionProduct boosted;
210 G4double aXsection; 229 G4double aXsection;
211 << 230
212 // MC integration loop 231 // MC integration loop
213 G4int counter = 0; 232 G4int counter = 0;
214 G4double buffer = 0; 233 G4double buffer = 0;
215 G4int size = G4int(std::max(10., aT / 60 * k << 234 G4int size = G4int(std::max(10., aT/60*kelvin));
216 G4ThreeVector neutronVelocity = << 235 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutronRP.GetMomentum();
217 1. / G4Neutron::Neutron()->GetPDGMass() * <<
218 G4double neutronVMag = neutronVelocity.mag() 236 G4double neutronVMag = neutronVelocity.mag();
219 237
220 while (counter == 0 << 238 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 { 239 {
224 if (counter != 0) buffer = result / counte << 240 if(counter) buffer = result/counter;
225 while (counter < size) // Loop checking, << 241 while (counter<size) // Loop checking, 11.05.2015, T. Koi
226 { 242 {
227 counter++; << 243 counter ++;
228 G4ReactionProduct aThermalNuc = aNuc.Get 244 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
229 boosted.Lorentz(theNeutronRP, aThermalNu 245 boosted.Lorentz(theNeutronRP, aThermalNuc);
230 G4double theEkin = boosted.GetKineticEne 246 G4double theEkin = boosted.GetKineticEnergy();
231 aXsection = (*((*theCrossSections)(index 247 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
232 // velocity correction. 248 // velocity correction.
233 G4ThreeVector targetVelocity = 1. / aThe << 249 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
234 aXsection *= (targetVelocity - neutronVe << 250 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
235 result += aXsection; 251 result += aXsection;
236 } 252 }
237 size += size; 253 size += size;
238 } 254 }
239 result /= counter; 255 result /= counter;
240 return result; 256 return result;
241 } 257 }
242 258
243 G4int G4ParticleHPFissionData::GetVerboseLevel << 259 G4int G4ParticleHPFissionData::GetVerboseLevel() const
244 { 260 {
245 return G4ParticleHPManager::GetInstance()->G << 261 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
246 } 262 }
247 263
248 void G4ParticleHPFissionData::SetVerboseLevel( << 264 void G4ParticleHPFissionData::SetVerboseLevel( G4int newValue )
249 { 265 {
250 G4ParticleHPManager::GetInstance()->SetVerbo << 266 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
251 } 267 }
252 268
253 void G4ParticleHPFissionData::CrossSectionDesc 269 void G4ParticleHPFissionData::CrossSectionDescription(std::ostream& outFile) const
254 { 270 {
255 outFile << "High Precision cross data based << 271 outFile << "High Precision cross data based on Evaluated Nuclear Data Files (ENDF)\n"
256 << "for induced fission reaction of << 272 << "for induced fission reaction of neutrons below 20MeV\n" ;
257 } 273 }
258 274