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