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