| 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 // particle_hp -- source file 26 // particle_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 // 070523 add neglecting doppler broadening on 30 // 070523 add neglecting doppler broadening on the fly. T. Koi
31 // 070613 fix memory leaking by T. Koi 31 // 070613 fix memory leaking by T. Koi
32 // 071002 enable cross section dump by T. Koi 32 // 071002 enable cross section dump by T. Koi
33 // 080428 change checking point of "neglecting << 33 // 080428 change checking point of "neglecting doppler broadening" flag
34 // from GetCrossSection to BuildPhysics 34 // from GetCrossSection to BuildPhysicsTable by T. Koi
35 // 081024 G4NucleiPropertiesTable:: to G4Nucle 35 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
36 // 36 //
37 // P. Arce, June-2014 Conversion neutron_hp to 37 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
38 // 38 //
39 #include "G4ParticleHPInelasticData.hh" 39 #include "G4ParticleHPInelasticData.hh"
40 << 40 #include "G4ParticleHPManager.hh"
41 #include "G4ElementTable.hh" <<
42 #include "G4HadronicParameters.hh" <<
43 #include "G4Neutron.hh" 41 #include "G4Neutron.hh"
44 #include "G4NucleiProperties.hh" << 42 #include "G4ElementTable.hh"
45 #include "G4ParticleHPData.hh" 43 #include "G4ParticleHPData.hh"
46 #include "G4ParticleHPManager.hh" << 44 #include "G4HadronicParameters.hh"
47 #include "G4Pow.hh" 45 #include "G4Pow.hh"
48 46
49 G4ParticleHPInelasticData::G4ParticleHPInelast 47 G4ParticleHPInelasticData::G4ParticleHPInelasticData(G4ParticleDefinition* projectile)
50 : G4VCrossSectionDataSet("") 48 : G4VCrossSectionDataSet("")
51 { 49 {
52 const char* dataDirVariable; 50 const char* dataDirVariable;
53 G4String particleName; 51 G4String particleName;
54 if (projectile == G4Neutron::Neutron()) { << 52 if( projectile == G4Neutron::Neutron() ) {
55 dataDirVariable = "G4NEUTRONHPDATA"; << 53 dataDirVariable = "G4NEUTRONHPDATA";
56 } << 54 }else if( projectile == G4Proton::Proton() ) {
57 else if (projectile == G4Proton::Proton()) { <<
58 dataDirVariable = "G4PROTONHPDATA"; 55 dataDirVariable = "G4PROTONHPDATA";
59 particleName = "Proton"; 56 particleName = "Proton";
60 } << 57 }else if( projectile == G4Deuteron::Deuteron() ) {
61 else if (projectile == G4Deuteron::Deuteron( <<
62 dataDirVariable = "G4DEUTERONHPDATA"; 58 dataDirVariable = "G4DEUTERONHPDATA";
63 particleName = "Deuteron"; 59 particleName = "Deuteron";
64 } << 60 }else if( projectile == G4Triton::Triton() ) {
65 else if (projectile == G4Triton::Triton()) { <<
66 dataDirVariable = "G4TRITONHPDATA"; 61 dataDirVariable = "G4TRITONHPDATA";
67 particleName = "Triton"; 62 particleName = "Triton";
68 } << 63 }else if( projectile == G4He3::He3() ) {
69 else if (projectile == G4He3::He3()) { <<
70 dataDirVariable = "G4HE3HPDATA"; 64 dataDirVariable = "G4HE3HPDATA";
71 particleName = "He3"; 65 particleName = "He3";
72 } << 66 }else if( projectile == G4Alpha::Alpha() ) {
73 else if (projectile == G4Alpha::Alpha()) { <<
74 dataDirVariable = "G4ALPHAHPDATA"; 67 dataDirVariable = "G4ALPHAHPDATA";
75 particleName = "Alpha"; 68 particleName = "Alpha";
76 } << 69 } else {
77 else { << 70 G4String message("G4ParticleHPInelasticData may only be called for neutron, proton, deuteron, triton, He3 or alpha, while it is called for " + projectile->GetParticleName());
78 G4String message( << 71 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
79 "G4ParticleHPInelasticData may only be c << 72 }
80 "alpha, while it is called for " << 73 // G4cout << this << " G4ParticleHPInelasticData::G4ParticleHPInelasticData " << projectile->GetParticleName() << " DATADIR " << dataDirVariable << G4endl;//GDEB
81 + projectile->GetParticleName()); << 74 G4String dataName = projectile->GetParticleName()+"HPInelasticXS";
82 throw G4HadronicException(__FILE__, __LINE << 75 dataName.at(0) = toupper(dataName.at(0)) ;
83 } << 76 SetName( dataName );
84 G4String dataName = projectile->GetParticleN << 77
85 dataName.at(0) = (char)std::toupper(dataName << 78 if ( !std::getenv(dataDirVariable) && !std::getenv( "G4PARTICLEHPDATA" ) ){
86 SetName(dataName); << 79 G4String message("Please setenv G4PARTICLEHPDATA (recommended) or, at least setenv " +
87 << 80 G4String(dataDirVariable) + " to point to the " + projectile->GetParticleName() + " cross-section files.");
88 if ((G4FindDataDir(dataDirVariable) == nullp << 81 throw G4HadronicException(__FILE__, __LINE__,message.c_str());
89 { <<
90 G4String message("Please setenv G4PARTICLE <<
91 + G4String(dataDirVariabl <<
92 + projectile->GetParticle <<
93 throw G4HadronicException(__FILE__, __LINE <<
94 } 82 }
95 83
96 G4String dirName; 84 G4String dirName;
97 if (G4FindDataDir(dataDirVariable) != nullpt << 85 if ( std::getenv(dataDirVariable) ) {
98 dirName = G4FindDataDir(dataDirVariable); << 86 dirName = std::getenv(dataDirVariable);
99 } << 87 } else {
100 else { << 88 G4String baseName = std::getenv( "G4PARTICLEHPDATA" );
101 G4String baseName = G4FindDataDir("G4PARTI << 89 dirName = baseName + "/" + particleName;
102 dirName = baseName + "/" + particleName; << 90 }
103 } << 91 #ifdef G4VERBOSE
104 #ifdef G4VERBOSE << 92 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) {
105 if (G4HadronicParameters::Instance()->GetVer << 93 G4cout << "@@@ G4ParticleHPInelasticData instantiated for particle " << projectile->GetParticleName() << " data directory variable is " << dataDirVariable << " pointing to " << dirName << G4endl;
106 G4cout << "@@@ G4ParticleHPInelasticData i << 94 }
107 << projectile->GetParticleName() << << 95 #endif
108 << " pointing to " << dirName << G4 << 96
109 } << 97 SetMinKinEnergy( 0*CLHEP::MeV );
110 #endif << 98 SetMaxKinEnergy( 20*CLHEP::MeV );
111 << 99
112 SetMinKinEnergy(0 * CLHEP::MeV); << 100 onFlightDB = true;
113 SetMaxKinEnergy(20 * CLHEP::MeV); << 101 theCrossSections = 0;
>> 102 theProjectile=projectile;
>> 103
>> 104 theHPData = NULL;
>> 105 instanceOfWorker = false;
>> 106 if ( G4Threading::IsMasterThread() ) {
>> 107 theHPData = new G4ParticleHPData( theProjectile );
>> 108 } else {
>> 109 instanceOfWorker = true;
>> 110 }
>> 111 element_cache = NULL;
>> 112 material_cache = NULL;
>> 113 ke_cache = 0.0;
>> 114 xs_cache = 0.0;
>> 115 }
>> 116
>> 117 G4ParticleHPInelasticData::~G4ParticleHPInelasticData()
>> 118 {
>> 119 if ( theCrossSections != NULL && instanceOfWorker != true ) {
>> 120 theCrossSections->clearAndDestroy();
>> 121 delete theCrossSections;
>> 122 theCrossSections = NULL;
>> 123 }
>> 124 if ( theHPData != NULL && instanceOfWorker != true ) {
>> 125 delete theHPData;
>> 126 theHPData = NULL;
>> 127 }
>> 128 }
114 129
115 theCrossSections = nullptr; << 130 G4bool G4ParticleHPInelasticData::IsIsoApplicable( const G4DynamicParticle* dp ,
116 theProjectile = projectile; << 131 G4int /*Z*/ , G4int /*A*/ ,
>> 132 const G4Element* /*elm*/ ,
>> 133 const G4Material* /*mat*/ )
>> 134 {
>> 135 G4double eKin = dp->GetKineticEnergy();
>> 136 if ( eKin > GetMaxKinEnergy()
>> 137 || eKin < GetMinKinEnergy()
>> 138 || dp->GetDefinition() != theProjectile ) return false;
117 139
118 theHPData = nullptr; << 140 return true;
119 instanceOfWorker = false; <<
120 if (G4Threading::IsMasterThread()) { <<
121 theHPData = new G4ParticleHPData(theProjec <<
122 } <<
123 else { <<
124 instanceOfWorker = true; <<
125 } <<
126 element_cache = nullptr; <<
127 material_cache = nullptr; <<
128 ke_cache = 0.0; <<
129 xs_cache = 0.0; <<
130 } 141 }
131 142
132 G4ParticleHPInelasticData::~G4ParticleHPInelas << 143 G4double G4ParticleHPInelasticData::GetIsoCrossSection( const G4DynamicParticle* dp ,
>> 144 G4int /*Z*/ , G4int /*A*/ ,
>> 145 const G4Isotope* /*iso*/ ,
>> 146 const G4Element* element ,
>> 147 const G4Material* material )
133 { 148 {
134 if (theCrossSections != nullptr && !instance << 149 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
135 theCrossSections->clearAndDestroy(); <<
136 delete theCrossSections; <<
137 theCrossSections = nullptr; <<
138 } <<
139 if (theHPData != nullptr && !instanceOfWorke <<
140 delete theHPData; <<
141 theHPData = nullptr; <<
142 } <<
143 } <<
144 150
145 G4bool G4ParticleHPInelasticData::IsIsoApplica << 151 ke_cache = dp->GetKineticEnergy();
146 << 152 element_cache = element;
147 << 153 material_cache = material;
148 { << 154 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
149 G4double eKin = dp->GetKineticEnergy(); << 155 xs_cache = xs;
150 return eKin <= GetMaxKinEnergy() && eKin >= << 156 return xs;
151 && dp->GetDefinition() == theProjecti <<
152 } 157 }
153 158
154 G4double G4ParticleHPInelasticData::GetIsoCros << 159 /*
155 << 160 G4bool G4ParticleHPInelasticData::IsApplicable(const G4DynamicParticle*aP, const G4Element*)
156 << 161 {
157 << 162 G4bool result = true;
158 { << 163 G4double eKin = aP->GetKineticEnergy();
159 if (dp->GetKineticEnergy() == ke_cache && el << 164 if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
160 return xs_cache; << 165 return result;
161 <<
162 ke_cache = dp->GetKineticEnergy(); <<
163 element_cache = element; <<
164 material_cache = material; <<
165 G4double xs = GetCrossSection(dp, element, m <<
166 xs_cache = xs; <<
167 return xs; <<
168 } 166 }
>> 167 */
169 168
170 void G4ParticleHPInelasticData::BuildPhysicsTa << 169 //void G4ParticleHPInelasticData::BuildPhysicsTableHP(G4ParticleDefinition* projectile,const char* /* dataDirVariable */)
>> 170 void G4ParticleHPInelasticData::BuildPhysicsTable( const G4ParticleDefinition& projectile )
171 { 171 {
172 if (G4Threading::IsWorkerThread()) { << 172 // if(&projectile!=G4Neutron::Neutron())
173 theCrossSections = G4ParticleHPManager::Ge << 173 // throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
174 return; << 174
175 } << 175 //080428
176 if (theHPData == nullptr) << 176 if ( G4ParticleHPManager::GetInstance()->GetNeglectDoppler() )
177 theHPData = G4ParticleHPData::Instance(con << 177 {
>> 178 onFlightDB = false;
>> 179 #ifdef G4VERBOSE
>> 180 if ( G4HadronicParameters::Instance()->GetVerboseLevel() > 0 ) {
>> 181 G4cout << "Find a flag of \"G4PHP_NEGLECT_DOPPLER\"." << G4endl;
>> 182 G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of inelastic scattering of neutrons (<20MeV)." << G4endl;
>> 183 }
>> 184 #endif
>> 185 }
>> 186
>> 187 if ( G4Threading::IsWorkerThread() ) {
>> 188 theCrossSections = G4ParticleHPManager::GetInstance()->GetInelasticCrossSections( &projectile );
>> 189 return;
>> 190 } else {
>> 191 if ( theHPData == NULL ) theHPData = G4ParticleHPData::Instance( const_cast<G4ParticleDefinition*> ( &projectile ) );
>> 192 }
>> 193
178 194
179 std::size_t numberOfElements = G4Element::Ge << 195
180 if (theCrossSections == nullptr) << 196 size_t numberOfElements = G4Element::GetNumberOfElements();
181 theCrossSections = new G4PhysicsTable(numb << 197 // theCrossSections = new G4PhysicsTable( numberOfElements );
>> 198 // TKDB
>> 199 //if ( theCrossSections == 0 )
>> 200 //{ theCrossSections = new G4PhysicsTable( numberOfElements ); }
>> 201 if ( theCrossSections == NULL )
>> 202 theCrossSections = new G4PhysicsTable( numberOfElements );
182 else 203 else
183 theCrossSections->clearAndDestroy(); 204 theCrossSections->clearAndDestroy();
184 205
185 // make a PhysicsVector for each element 206 // make a PhysicsVector for each element
186 207
187 auto theElementTable = G4Element::GetElement << 208 //G4ParticleHPData* hpData = new G4ParticleHPData(projectile); //NEW
188 for (std::size_t i = 0; i < numberOfElements << 209 static G4ThreadLocal G4ElementTable *theElementTable = 0 ;
>> 210 if (!theElementTable) theElementTable= G4Element::GetElementTable();
>> 211 for( size_t i=0; i<numberOfElements; ++i )
>> 212 {
>> 213 //NEW G4PhysicsVector* physVec = G4ParticleHPData::
>> 214 //NEW Instance(projectile, dataDirVariable)->MakePhysicsVector((*theElementTable)[i], this);
>> 215 //G4PhysicsVector* physVec = hpData->MakePhysicsVector((*theElementTable)[i], this);
189 G4PhysicsVector* physVec = theHPData->Make 216 G4PhysicsVector* physVec = theHPData->MakePhysicsVector((*theElementTable)[i], this);
190 theCrossSections->push_back(physVec); 217 theCrossSections->push_back(physVec);
191 } 218 }
192 G4ParticleHPManager::GetInstance()->Register << 219
>> 220 G4ParticleHPManager::GetInstance()->RegisterInelasticCrossSections( &projectile , theCrossSections );
193 } 221 }
194 222
195 void G4ParticleHPInelasticData::DumpPhysicsTab << 223 void G4ParticleHPInelasticData::DumpPhysicsTable(const G4ParticleDefinition& projectile)
196 { 224 {
197 #ifdef G4VERBOSE << 225 if(&projectile!=theProjectile)
198 if (G4HadronicParameters::Instance()->GetVer << 226 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use ParticleHP data for a wrong projectile!!!");
>> 227
>> 228 #ifdef G4VERBOSE
>> 229 if ( G4HadronicParameters::Instance()->GetVerboseLevel() == 0 ) return;
>> 230 //
>> 231 // Dump element based cross section
>> 232 // range 10e-5 eV to 20 MeV
>> 233 // 10 point per decade
>> 234 // in barn
>> 235 //
>> 236
>> 237 G4cout << G4endl;
>> 238 G4cout << G4endl;
>> 239 G4cout << "Inelastic Cross Section of Neutron HP"<< G4endl;
>> 240 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
>> 241 G4cout << G4endl;
>> 242 G4cout << "Name of Element" << G4endl;
>> 243 G4cout << "Energy[eV] XS[barn]" << G4endl;
>> 244 G4cout << G4endl;
>> 245
>> 246 size_t numberOfElements = G4Element::GetNumberOfElements();
>> 247 static G4ThreadLocal G4ElementTable *theElementTable = 0 ;
>> 248 if (!theElementTable) theElementTable= G4Element::GetElementTable();
>> 249
>> 250 for ( size_t i = 0 ; i < numberOfElements ; ++i )
>> 251 {
>> 252
>> 253 G4cout << (*theElementTable)[i]->GetName() << G4endl;
>> 254
>> 255 G4int ie = 0;
>> 256
>> 257 for ( ie = 0 ; ie < 130 ; ie++ )
>> 258 {
>> 259 G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA ( 10.0 , ie/10.0 ) *CLHEP::eV;
>> 260 G4bool outOfRange = false;
>> 261
>> 262 if ( eKinetic < 20*CLHEP::MeV )
>> 263 {
>> 264 G4cout << eKinetic/CLHEP::eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/CLHEP::barn << G4endl;
>> 265 }
199 266
200 // Dump element based cross section <<
201 // range 10e-5 eV to 20 MeV <<
202 // 10 point per decade <<
203 // in barn <<
204 <<
205 G4cout << G4endl; <<
206 G4cout << G4endl; <<
207 G4cout << "Inelastic Cross Section of Neutro <<
208 G4cout << "(Pointwise cross-section at 0 Kel <<
209 G4cout << G4endl; <<
210 G4cout << "Name of Element" << G4endl; <<
211 G4cout << "Energy[eV] XS[barn]" << G4endl; <<
212 G4cout << G4endl; <<
213 <<
214 std::size_t numberOfElements = G4Element::Ge <<
215 auto theElementTable = G4Element::GetElement <<
216 <<
217 for (std::size_t i = 0; i < numberOfElements <<
218 G4cout << (*theElementTable)[i]->GetName() <<
219 <<
220 G4int ie = 0; <<
221 <<
222 for (ie = 0; ie < 130; ie++) { <<
223 G4double eKinetic = 1.0e-5 * G4Pow::GetI <<
224 G4bool outOfRange = false; <<
225 <<
226 if (eKinetic < 20 * CLHEP::MeV) { <<
227 G4cout << eKinetic / CLHEP::eV << " " <<
228 << (*((*theCrossSections)(i))). <<
229 << G4endl; <<
230 } 267 }
231 } << 268
232 G4cout << G4endl; << 269 G4cout << G4endl;
233 } << 270 }
234 #endif << 271
>> 272 //G4cout << "G4ParticleHPInelasticData::DumpPhysicsTable still to be implemented"<<G4endl;
>> 273 #endif
235 } 274 }
236 275
237 G4double G4ParticleHPInelasticData::GetCrossSe << 276 #include "G4NucleiProperties.hh"
238 << 277
>> 278 G4double G4ParticleHPInelasticData::
>> 279 GetCrossSection(const G4DynamicParticle* projectile, const G4Element*anE, G4double aT)
239 { 280 {
240 G4double result = 0; 281 G4double result = 0;
241 G4bool outOfRange; 282 G4bool outOfRange;
242 std::size_t index = anE->GetIndex(); << 283 G4int index = anE->GetIndex();
243 284
244 // prepare neutron 285 // prepare neutron
245 G4double eKinetic = projectile->GetKineticEn 286 G4double eKinetic = projectile->GetKineticEnergy();
246 287
247 if (G4ParticleHPManager::GetInstance()->GetN << 288 if ( !onFlightDB )
248 // NEGLECT_DOPPLER << 289 {
249 G4double factor = 1.0; << 290 //NEGLECT_DOPPLER
250 if (eKinetic < aT * CLHEP::k_Boltzmann) { << 291 G4double factor = 1.0;
251 // below 0.1 eV neutrons << 292 if ( eKinetic < aT * CLHEP::k_Boltzmann )
252 // Have to do some, but now just igonre. << 293 {
253 // Will take care after performance chec << 294 // below 0.1 eV neutrons
254 // factor = factor * targetV; << 295 // Have to do some, but now just igonre.
255 } << 296 // Will take care after performance check.
256 return ((*((*theCrossSections)(index))).Ge << 297 // factor = factor * targetV;
257 } << 298 }
258 << 299 return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
259 G4ReactionProduct theNeutron(projectile->Get << 300
260 theNeutron.SetMomentum(projectile->GetMoment << 301 }
261 theNeutron.SetKineticEnergy(eKinetic); << 302
>> 303 G4ReactionProduct theNeutron( projectile->GetDefinition() );
>> 304 theNeutron.SetMomentum( projectile->GetMomentum() );
>> 305 theNeutron.SetKineticEnergy( eKinetic );
262 306
263 // prepare thermal nucleus 307 // prepare thermal nucleus
264 G4Nucleus aNuc; 308 G4Nucleus aNuc;
265 G4double eps = 0.0001; 309 G4double eps = 0.0001;
266 G4double theA = anE->GetN(); 310 G4double theA = anE->GetN();
267 G4double theZ = anE->GetZ(); 311 G4double theZ = anE->GetZ();
268 G4double eleMass; << 312 G4double eleMass;
269 eleMass = G4NucleiProperties::GetNuclearMass << 313 eleMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps) );
270 << 314
271 <<
272 G4ReactionProduct boosted; 315 G4ReactionProduct boosted;
273 G4double aXsection; 316 G4double aXsection;
274 << 317
275 // MC integration loop 318 // MC integration loop
276 G4int counter = 0; 319 G4int counter = 0;
277 G4int failCount = 0; 320 G4int failCount = 0;
278 G4double buffer = 0; << 321 G4double buffer = 0; G4int size = G4int(std::max(10., aT/60*CLHEP::kelvin));
279 G4int size = G4int(std::max(10., aT / 60 * C << 322 G4ThreeVector neutronVelocity = 1./theProjectile->GetPDGMass()*theNeutron.GetMomentum();
280 G4ThreeVector neutronVelocity = 1. / theProj <<
281 G4double neutronVMag = neutronVelocity.mag() 323 G4double neutronVMag = neutronVelocity.mag();
282 324
>> 325 // G4cout << " G4ParticleHPInelasticData 2 " << size << G4endl;//GDEB
283 #ifndef G4PHP_DOPPLER_LOOP_ONCE 326 #ifndef G4PHP_DOPPLER_LOOP_ONCE
284 while (counter == 0 << 327 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer) // Loop checking, 11.05.2015, T. Koi
285 || std::abs(buffer - result / std::ma <<
286 > 0.01 * buffer) // Loop checki <<
287 { 328 {
288 if (counter != 0) buffer = result / counte << 329 if(counter) buffer = result/counter;
289 while (counter < size) // Loop checking, << 330 while (counter<size) // Loop checking, 11.05.2015, T. Koi
290 { 331 {
291 ++counter; << 332 counter ++;
292 #endif 333 #endif
293 G4ReactionProduct aThermalNuc = << 334 //G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass/theProjectile->GetPDGMass(), aT );
294 aNuc.GetThermalNucleus(eleMass / G4Neu << 335 //G4Nucleus::GetThermalNucleus requests normalization of mass in neutron mass
>> 336 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus( eleMass/G4Neutron::Neutron()->GetPDGMass(), aT );
295 boosted.Lorentz(theNeutron, aThermalNuc) 337 boosted.Lorentz(theNeutron, aThermalNuc);
296 G4double theEkin = boosted.GetKineticEne 338 G4double theEkin = boosted.GetKineticEnergy();
297 aXsection = (*((*theCrossSections)(index 339 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
298 if (aXsection < 0) { << 340 // G4cout << " G4ParticleHPInelasticData aXsection " << aXsection << " index " << index << " theEkin " << theEkin << " outOfRange " << outOfRange <<G4endl;//GDEB
299 if (failCount < 1000) { << 341 if(aXsection <0)
300 ++failCount; << 342 {
>> 343 if(failCount<1000)
>> 344 {
>> 345 failCount++;
301 #ifndef G4PHP_DOPPLER_LOOP_ONCE 346 #ifndef G4PHP_DOPPLER_LOOP_ONCE
302 --counter; << 347 counter--;
303 continue; << 348 continue;
304 #endif 349 #endif
305 } << 350 }
306 else { << 351 else
307 aXsection = 0; << 352 {
308 } << 353 aXsection = 0;
>> 354 }
309 } 355 }
310 // velocity correction. 356 // velocity correction.
311 G4ThreeVector targetVelocity = 1. / aThe << 357 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
312 aXsection *= (targetVelocity - neutronVe << 358 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
313 result += aXsection; 359 result += aXsection;
314 } 360 }
315 #ifndef G4PHP_DOPPLER_LOOP_ONCE 361 #ifndef G4PHP_DOPPLER_LOOP_ONCE
316 size += size; 362 size += size;
317 } 363 }
318 result /= counter; 364 result /= counter;
319 #endif 365 #endif
320 366
>> 367 /*
>> 368 // Checking impact of G4PHP_NEGLECT_DOPPLER
>> 369 G4cout << " result " << result << " "
>> 370 << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " "
>> 371 << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl;
>> 372 */
>> 373 // G4cout << this << " G4ParticleHPInelasticData result " << result << G4endl; //GDEB
>> 374
321 return result; 375 return result;
322 } 376 }
323 377
324 G4int G4ParticleHPInelasticData::GetVerboseLev << 378 G4int G4ParticleHPInelasticData::GetVerboseLevel() const
325 { 379 {
326 return G4ParticleHPManager::GetInstance()->G << 380 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
327 } 381 }
328 << 382 void G4ParticleHPInelasticData::SetVerboseLevel( G4int newValue )
329 void G4ParticleHPInelasticData::SetVerboseLeve <<
330 { 383 {
331 G4ParticleHPManager::GetInstance()->SetVerbo << 384 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
332 } 385 }
333 <<
334 void G4ParticleHPInelasticData::CrossSectionDe 386 void G4ParticleHPInelasticData::CrossSectionDescription(std::ostream& outFile) const
335 { 387 {
336 outFile << "Extension of High Precision cros << 388 outFile << "Extension of High Precision cross section for inelastic reaction of proton, deuteron, triton, He3 and alpha below 20MeV\n";
337 "deuteron, triton, He3 and alpha <<
338 } 389 }
339 390