Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPInelasticData.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/particle_hp/src/G4ParticleHPInelasticData.cc (Version 11.3.0) and /processes/hadronic/models/particle_hp/src/G4ParticleHPInelasticData.cc (Version 10.5.p1)


  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