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