Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4NeutronInelasticXS.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/cross_sections/src/G4NeutronInelasticXS.cc (Version 11.3.0) and /processes/hadronic/cross_sections/src/G4NeutronInelasticXS.cc (Version 10.3)


  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 // $Id: G4NeutronInelasticXS.cc 93682 2015-10-28 10:09:49Z gcosmo $
                                                   >>  27 //
 26 // -------------------------------------------     28 // -------------------------------------------------------------------
 27 //                                                 29 //
 28 // GEANT4 Class file                               30 // GEANT4 Class file
 29 //                                                 31 //
 30 //                                                 32 //
 31 // File name:    G4NeutronInelasticXS              33 // File name:    G4NeutronInelasticXS
 32 //                                                 34 //
 33 // Author  Ivantchenko, Geant4, 3-Aug-09           35 // Author  Ivantchenko, Geant4, 3-Aug-09
 34 //                                                 36 //
                                                   >>  37 // Modifications:
                                                   >>  38 //
 35                                                    39 
 36 #include "G4NeutronInelasticXS.hh"                 40 #include "G4NeutronInelasticXS.hh"
 37 #include "G4Neutron.hh"                            41 #include "G4Neutron.hh"
 38 #include "G4DynamicParticle.hh"                    42 #include "G4DynamicParticle.hh"
 39 #include "G4ElementTable.hh"                   << 
 40 #include "G4Material.hh"                       << 
 41 #include "G4Element.hh"                            43 #include "G4Element.hh"
                                                   >>  44 #include "G4ElementTable.hh"
 42 #include "G4PhysicsLogVector.hh"                   45 #include "G4PhysicsLogVector.hh"
 43 #include "G4CrossSectionDataSetRegistry.hh"    <<  46 #include "G4PhysicsVector.hh"
 44 #include "G4ComponentGGHadronNucleusXsc.hh"        47 #include "G4ComponentGGHadronNucleusXsc.hh"
 45 #include "G4HadronicParameters.hh"             <<  48 #include "G4HadronNucleonXsc.hh"
                                                   >>  49 #include "G4NistManager.hh"
                                                   >>  50 #include "G4Proton.hh"
 46 #include "Randomize.hh"                            51 #include "Randomize.hh"
 47 #include "G4Neutron.hh"                        << 
 48 #include "G4SystemOfUnits.hh"                  << 
 49 #include "G4IsotopeList.hh"                    << 
 50 #include "G4NuclearRadii.hh"                   << 
 51 #include "G4AutoLock.hh"                       << 
 52                                                    52 
                                                   >>  53 #include <iostream>
 53 #include <fstream>                                 54 #include <fstream>
 54 #include <sstream>                                 55 #include <sstream>
 55 #include <thread>                              << 
 56                                                    56 
 57 G4double G4NeutronInelasticXS::coeff[] = {1.0} <<  57 // factory
 58 G4double G4NeutronInelasticXS::lowcoeff[] = {1 <<  58 #include "G4CrossSectionFactory.hh"
 59 G4ElementData* G4NeutronInelasticXS::data = nu <<  59 //
 60 G4String G4NeutronInelasticXS::gDataDirectory  <<  60 G4_DECLARE_XS_FACTORY(G4NeutronInelasticXS);
 61                                                    61 
 62 static std::once_flag applyOnce;               <<  62 using namespace std;
 63                                                    63 
 64 namespace                                      <<  64 const G4int G4NeutronInelasticXS::amin[] = {
 65 {                                              <<  65   0,
 66   G4Mutex nInelasticXSMutex = G4MUTEX_INITIALI <<  66   0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
 67 }                                              <<  67   0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
                                                   >>  68   0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
                                                   >>  69   0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
                                                   >>  70   0, 0, 0, 0, 0, 0,107,106, 0,112, //41-50
                                                   >>  71   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
                                                   >>  72   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
                                                   >>  73   0, 0, 0,180, 0, 0, 0, 0, 0, 0, //71-80
                                                   >>  74   0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
                                                   >>  75   0,235};
                                                   >>  76 const G4int G4NeutronInelasticXS::amax[] = {
                                                   >>  77   0,
                                                   >>  78   0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
                                                   >>  79   0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
                                                   >>  80   0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
                                                   >>  81   0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
                                                   >>  82   0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
                                                   >>  83   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
                                                   >>  84   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
                                                   >>  85   0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
                                                   >>  86   0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
                                                   >>  87   0,238};
                                                   >>  88 
                                                   >>  89 G4double G4NeutronInelasticXS::coeff[] = {1.0};
                                                   >>  90 
                                                   >>  91 G4ElementData* G4NeutronInelasticXS::data = 0;
 68                                                    92 
 69 G4NeutronInelasticXS::G4NeutronInelasticXS()       93 G4NeutronInelasticXS::G4NeutronInelasticXS() 
 70   : G4VCrossSectionDataSet(Default_Name()),        94   : G4VCrossSectionDataSet(Default_Name()),
 71     neutron(G4Neutron::Neutron()),             <<  95     proton(G4Proton::Proton())
 72     elimit(20*CLHEP::MeV),                     << 
 73     lowElimit(1.0e-7*CLHEP::eV)                << 
 74 {                                                  96 {
 75   verboseLevel = 0;                            <<  97   //  verboseLevel = 0;
 76   if (verboseLevel > 0) {                      <<  98   if(verboseLevel > 0){
 77     G4cout << "G4NeutronInelasticXS::G4Neutron     99     G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < " 
 78             << MAXZINEL << G4endl;             << 100       << MAXZINEL << G4endl;
 79   }                                               101   }
 80   loglowElimit = G4Log(lowElimit);             << 102   ggXsection = new G4ComponentGGHadronNucleusXsc();
 81   if (nullptr == data) {                       << 103   fNucleon = new G4HadronNucleonXsc();
 82     data = new G4ElementData(MAXZINEL);        << 104   isMaster = false;
 83     data->SetName("nInelastic");               << 105 }
 84     FindDirectoryPath();                       << 
 85   }                                            << 
 86   ggXsection =                                 << 
 87     G4CrossSectionDataSetRegistry::Instance()- << 
 88   if(ggXsection == nullptr)                    << 
 89     ggXsection = new G4ComponentGGHadronNucleu << 
 90                                                   106 
 91   SetForAllAtomsAndEnergies(true);             << 107 G4NeutronInelasticXS::~G4NeutronInelasticXS()
                                                   >> 108 {
                                                   >> 109   //G4cout << "G4NeutronInelasticXS::~G4NeutronInelasticXS() " 
                                                   >> 110   // << " isMaster= " << isMaster << "  data: " << data << G4endl;
                                                   >> 111   delete fNucleon;
                                                   >> 112   delete ggXsection;
                                                   >> 113   if(isMaster) { delete data; data = 0; }
 92 }                                                 114 }
 93                                                   115 
 94 void G4NeutronInelasticXS::CrossSectionDescrip    116 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
 95 {                                                 117 {
 96   outFile << "G4NeutronInelasticXS calculates     118   outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
 97           << "cross section on nuclei using da    119           << "cross section on nuclei using data from the high precision\n"
 98           << "neutron database.  These data ar    120           << "neutron database.  These data are simplified and smoothed over\n"
 99           << "the resonance region in order to    121           << "the resonance region in order to reduce CPU time.\n"
100           << "For high energy Glauber-Gribov c << 122           << "G4NeutronInelasticXS is valid for energies up to 20 MeV, for\n"
                                                   >> 123           << "nuclei through U.\n";
101 }                                                 124 }
102                                                   125 
103 G4bool                                            126 G4bool 
104 G4NeutronInelasticXS::IsElementApplicable(cons    127 G4NeutronInelasticXS::IsElementApplicable(const G4DynamicParticle*, 
105                                           G4in << 128             G4int, const G4Material*)
106 {                                                 129 {
107   return true;                                    130   return true;
108 }                                                 131 }
109                                                   132 
110 G4bool                                            133 G4bool 
111 G4NeutronInelasticXS::IsIsoApplicable(const G4    134 G4NeutronInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
112                                       G4int, G << 135               G4int /*ZZ*/, G4int /*AA*/,
113                                       const G4 << 136               const G4Element*, const G4Material*)
114 {                                                 137 {
115   return true;                                    138   return true;
116 }                                                 139 }
117                                                   140 
118 G4double                                       << 141 G4double G4NeutronInelasticXS::GetElementCrossSection(
119 G4NeutronInelasticXS::GetElementCrossSection(c << 142          const G4DynamicParticle* aParticle,
120                                              G << 143    G4int Z, const G4Material*)
121 {                                                 144 {
122   return ElementCrossSection(aParticle->GetKin << 145   G4double xs = 0.0;
123                              aParticle->GetLog << 146   G4double ekin = aParticle->GetKineticEnergy();
124 }                                              << 
125                                                   147 
126 G4double                                       << 148   if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
127 G4NeutronInelasticXS::ComputeCrossSectionPerEl << 149   else if(Z < 1)    { Z = 1; }
128                                                << 
129                                                << 
130                                                << 
131 {                                              << 
132   return ElementCrossSection(ekin, loge, elm-> << 
133 }                                              << 
134                                                   150 
135 G4double                                       << 151   G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
136 G4NeutronInelasticXS::ElementCrossSection(G4do << 
137 {                                              << 
138   G4int Z = std::min(ZZ, MAXZINEL-1);          << 
139   G4double ekin = eKin;                        << 
140   G4double loge = logE;                        << 
141   G4double xs;                                 << 
142                                                   152 
143   // very low energy limit                     << 153   G4PhysicsVector* pv = data->GetElementData(Z);
144   if (ekin < lowElimit) {                      << 154   //  G4cout  << "G4NeutronInelasticXS::GetCrossSection e= " << ekin 
145     ekin = lowElimit;                          << 155   //  << " Z= " << Z << G4endl;
146     loge = loglowElimit;                       << 156 
                                                   >> 157   // element was not initialised
                                                   >> 158   if(!pv) {
                                                   >> 159     Initialise(Z);
                                                   >> 160     pv = data->GetElementData(Z);
                                                   >> 161     if(!pv) { return xs; }
147   }                                               162   }
148   // pv should exist                           << 
149   auto pv = GetPhysicsVector(Z);               << 
150                                                   163 
151   const G4double e0 = pv->Energy(0);           << 164   G4double e1 = pv->Energy(0);
152   if (ekin <= e0) {                            << 165   if(ekin <= e1) { return xs; }
153     xs = (*pv)[0];                             << 166 
154     if (xs > 0.0) { xs *= std::sqrt(e0/ekin);  << 167   G4double e2 = pv->GetMaxEnergy();
155   } else if (ekin <= pv->GetMaxEnergy()) {     << 168 
156     xs = pv->LogVectorValue(ekin, loge);       << 169   if(ekin <= e2) { 
157   } else {                                     << 170     xs = pv->Value(ekin); 
158     xs = coeff[Z]*ggXsection->GetInelasticElem << 171   } else if(1 == Z) { 
159                                                << 172     fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
                                                   >> 173     xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
                                                   >> 174   } else {          
                                                   >> 175     ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
                                                   >> 176     xs = coeff[Z]*ggXsection->GetInelasticGlauberGribovXsc();
160   }                                               177   }
161                                                   178 
162 #ifdef G4VERBOSE                               << 179   if(verboseLevel > 0) {
163   if(verboseLevel > 1) {                       << 180     G4cout  << "ekin= " << ekin << ",  XSinel= " << xs << G4endl;
164     G4cout  << "G4NeutronInelasticXS::ElementC << 
165             << " Ekin(MeV)= " << ekin/CLHEP::M << 
166             << ", ElmXSinel(b)= " << xs/CLHEP: << 
167             << G4endl;                         << 
168   }                                               181   }
169 #endif                                         << 
170   return xs;                                      182   return xs;
171 }                                                 183 }
172                                                   184 
173 G4double                                       << 185 G4double G4NeutronInelasticXS::GetIsoCrossSection(
174 G4NeutronInelasticXS::ComputeIsoCrossSection(G << 186          const G4DynamicParticle* aParticle, 
175                                              c << 187    G4int Z, G4int A,
176                                              G << 188    const G4Isotope*, const G4Element*,
177                                              c << 189    const G4Material*)
178                                              c << 190 {
179 {                                              << 191   if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
180   return IsoCrossSection(ekin, loge, Z, A);    << 192   else if(Z < 1)    { Z = 1; }
181 }                                              << 
182                                                << 
183 G4double                                       << 
184 G4NeutronInelasticXS::GetIsoCrossSection(const << 
185                                          G4int << 
186                                          const << 
187                                          const << 
188 {                                              << 
189   return IsoCrossSection(aParticle->GetKinetic << 
190                          aParticle->GetLogKine << 
191 }                                              << 
192                                                << 
193 G4double                                       << 
194 G4NeutronInelasticXS::IsoCrossSection(G4double << 
195                                       G4int ZZ << 
196 {                                              << 
197   G4double xs;                                 << 
198   G4int Z = std::min(ZZ, MAXZINEL-1);          << 
199   G4double ekin = eKin;                        << 
200   G4double loge = logE;                        << 
201                                                << 
202   /*                                           << 
203   G4cout << "G4NeutronInelasticXS::IsoCrossSec << 
204          << Z << "  A= " << A << G4endl;       << 
205   G4cout << "  Amin= " << amin[Z] << " Amax= " << 
206          << " E(MeV)= " << ekin << " Ncomp="   << 
207          << data->GetNumberOfComponents(Z) <<  << 
208   */                                           << 
209   GetPhysicsVector(Z);                         << 
210                                                   193 
211   // use isotope cross section if applicable   << 194   return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
212   if (ekin <= elimit && data->GetNumberOfCompo << 195 }
213     auto pviso = data->GetComponentDataByID(Z, << 
214     if (nullptr != pviso) {                    << 
215       const G4double e0 = pviso->Energy(0);    << 
216       if (ekin > e0) {                         << 
217         xs = pviso->LogVectorValue(ekin, loge) << 
218       } else {                                 << 
219         xs = (*pviso)[0];                      << 
220         if (xs > 0.0) { xs *= std::sqrt(e0/eki << 
221       }                                        << 
222 #ifdef G4VERBOSE                               << 
223       if(verboseLevel > 1) {                   << 
224         G4cout << "G4NeutronInelasticXS::IsoXS << 
225                << ekin/CLHEP::MeV              << 
226                << "  xs(b)= " << xs/CLHEP::bar << 
227                << "  Z= " << Z << "  A= " << A << 
228       }                                        << 
229 #endif                                         << 
230       return xs;                               << 
231     }                                          << 
232   }                                            << 
233                                                   196 
234   // use element x-section                     << 197 G4double 
235   xs = ElementCrossSection(ekin, loge, Z)*A/ae << 198 G4NeutronInelasticXS::IsoCrossSection(G4double ekin, G4int Z, G4int A)
                                                   >> 199 {
                                                   >> 200   G4double xs = 0.0;
                                                   >> 201   /*
                                                   >> 202   G4cout << "IsoCrossSection  Z= " << Z << "  A= " << A 
                                                   >> 203          << "  Amin= " << amin[Z] << " Amax= " << amax[Z]
                                                   >> 204          << " E(MeV)= " << ekin << G4endl;
                                                   >> 205   */
                                                   >> 206   // element was not initialised
                                                   >> 207   G4PhysicsVector* pv = data->GetElementData(Z);
                                                   >> 208   if(!pv) {
                                                   >> 209     Initialise(Z);
                                                   >> 210     pv = data->GetElementData(Z);
                                                   >> 211   }
236                                                   212 
237 #ifdef G4VERBOSE                               << 213   // isotope cross section exist
238   if(verboseLevel > 1) {                       << 214   if(pv && amin[Z] > 0  && A >= amin[Z] && A <= amax[Z]) {
239     G4cout  << "G4NeutronInelasticXS::IsoXS: Z << 215     pv = data->GetComponentDataByIndex(Z, A - amin[Z]);
240             << " Ekin(MeV)= " << ekin/CLHEP::M << 216     if(pv && ekin > pv->Energy(0)) { xs = pv->Value(ekin); }
241             << ", ElmXS(b)= " << xs/CLHEP::bar << 217   }
                                                   >> 218   if(verboseLevel > 0){
                                                   >> 219     G4cout  << "ekin= " << ekin << ",  xs= " << xs << G4endl;
242   }                                               220   }
243 #endif                                         << 
244   return xs;                                      221   return xs;
245 }                                                 222 }
246                                                   223 
247 const G4Isotope* G4NeutronInelasticXS::SelectI << 224 G4Isotope* G4NeutronInelasticXS::SelectIsotope(
248       const G4Element* anElement, G4double kin << 225            const G4Element* anElement, G4double kinEnergy)
249 {                                                 226 {
250   std::size_t nIso = anElement->GetNumberOfIso << 227   size_t nIso = anElement->GetNumberOfIsotopes();
251   const G4Isotope* iso = anElement->GetIsotope << 228   G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
252   if(1 == nIso) { return iso; }                << 229   G4Isotope* iso = (*isoVector)[0];
253                                                   230 
254   // more than 1 isotope                       << 231   //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
255   G4int Z = anElement->GetZasInt();            << 
256   if (nullptr == data->GetElementData(Z)) { In << 
257                                                   232 
258   const G4double* abundVector = anElement->Get << 233   // more than 1 isotope
259   G4double q = G4UniformRand();                << 234   if(1 < nIso) {
260   G4double sum = 0.0;                          << 235     G4int Z = G4lrint(anElement->GetZ());
261   std::size_t j;                               << 236     //G4cout << "SelectIsotope Z= " << Z << G4endl;
262                                                << 237 
263   // isotope wise cross section not available  << 238     G4double* abundVector = anElement->GetRelativeAbundanceVector();
264   if (Z >= MAXZINEL || 0 == data->GetNumberOfC << 239     G4double q = G4UniformRand();
265     for (j=0; j<nIso; ++j) {                   << 240     G4double sum = 0.0;
266       sum += abundVector[j];                   << 241 
267       if(q <= sum) {                           << 242     // is there isotope wise cross section?
268         iso = anElement->GetIsotope((G4int)j); << 243     size_t j;
269         break;                                 << 244     if(0 == amin[Z] || Z >= MAXZINEL) {
                                                   >> 245       for (j = 0; j<nIso; ++j) {
                                                   >> 246   sum += abundVector[j];
                                                   >> 247   if(q <= sum) {
                                                   >> 248     iso = (*isoVector)[j];
                                                   >> 249     break;
                                                   >> 250   }
270       }                                           251       }
271     }                                          << 252     } else {
272     return iso;                                << 
273   }                                            << 
274                                                   253 
275   // use isotope cross sections                << 254       // element may be not initialised in unit test
276   auto nn = temp.size();                       << 255       if(!data->GetElementData(Z)) { Initialise(Z); }
277   if(nn < nIso) { temp.resize(nIso, 0.); }     << 256       size_t nn = temp.size();
278                                                << 257       if(nn < nIso) { temp.resize(nIso, 0.); }
279   for (j=0; j<nIso; ++j) {                     << 258 
280     // G4cout << j << "-th isotope " << anElem << 259       for (j=0; j<nIso; ++j) {
281     //       <<  " abund= " << abundVector[j]  << 260   //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN() 
282     sum += abundVector[j]*IsoCrossSection(kinE << 261   //       <<  " abund= " << abundVector[j] << G4endl;
283                                           anEl << 262         sum += abundVector[j]*IsoCrossSection(kinEnergy, Z, 
284     temp[j] = sum;                             << 263                 (*isoVector)[j]->GetN());
285   }                                            << 264         temp[j] = sum;
286   sum *= q;                                    << 265       }
287   for (j = 0; j<nIso; ++j) {                   << 266       sum *= q;
288     if (temp[j] >= sum) {                      << 267       for (j = 0; j<nIso; ++j) {
289       iso = anElement->GetIsotope((G4int)j);   << 268         if(temp[j] >= sum) {
290       break;                                   << 269           iso = (*isoVector)[j];
                                                   >> 270           break;
                                                   >> 271   }
                                                   >> 272       }
291     }                                             273     }
292   }                                               274   }
293   return iso;                                     275   return iso;
294 }                                                 276 }
295                                                   277 
296 void                                              278 void 
297 G4NeutronInelasticXS::BuildPhysicsTable(const     279 G4NeutronInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
298 {                                                 280 {
299   if (verboseLevel > 0) {                      << 281   if(verboseLevel > 0){
300     G4cout << "G4NeutronInelasticXS::BuildPhys    282     G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for " 
301            << p.GetParticleName() << G4endl;   << 283      << p.GetParticleName() << G4endl;
302   }                                               284   }
303   if (p.GetParticleName() != "neutron") {      << 285   if(p.GetParticleName() != "neutron") { 
304     G4ExceptionDescription ed;                    286     G4ExceptionDescription ed;
305     ed << p.GetParticleName() << " is a wrong     287     ed << p.GetParticleName() << " is a wrong particle type -"
306        << " only neutron is allowed";             288        << " only neutron is allowed";
307     G4Exception("G4NeutronInelasticXS::BuildPh    289     G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
308                 FatalException, ed, "");       << 290     FatalException, ed, "");
309     return;                                       291     return; 
310   }                                               292   }
311   // it is possible re-initialisation for the  << 
312   const G4ElementTable* table = G4Element::Get << 
313                                                   293 
314   // initialise static tables only once        << 294   if(!data) { 
315   std::call_once(applyOnce, [this]() { isIniti << 295     isMaster = true;
                                                   >> 296     data = new G4ElementData(); 
                                                   >> 297     data->SetName("NeutronInelastic");
                                                   >> 298     temp.resize(13,0.0);
                                                   >> 299   }
316                                                   300 
317   if (isInitializer) {                         << 301   // it is possible re-initialisation for the new run
318     G4AutoLock l(&nInelasticXSMutex);          << 302   if(isMaster) {
319                                                   303 
320     // Upload data for elements used in geomet << 304     // check environment variable 
321     for ( auto const & elm : *table ) {        << 305     // Build the complete string identifying the file with the data set
322       G4int Z = std::max( 1, std::min( elm->Ge << 306     char* path = getenv("G4NEUTRONXSDATA");
323       if ( nullptr == data->GetElementData(Z)  << 307 
                                                   >> 308     G4DynamicParticle* dynParticle = 
                                                   >> 309       new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
                                                   >> 310 
                                                   >> 311     // Access to elements
                                                   >> 312     const G4ElementTable* theElmTable = G4Element::GetElementTable();
                                                   >> 313     size_t numOfElm = G4Element::GetNumberOfElements();
                                                   >> 314     if(numOfElm > 0) {
                                                   >> 315       for(size_t i=0; i<numOfElm; ++i) {
                                                   >> 316   G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
                                                   >> 317   if(Z < 1)              { Z = 1; }
                                                   >> 318   else if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
                                                   >> 319   //G4cout << "Z= " << Z << G4endl;
                                                   >> 320   // Initialisation 
                                                   >> 321   if(!(data->GetElementData(Z))) { 
                                                   >> 322     Initialise(Z, dynParticle, path); 
                                                   >> 323   }
                                                   >> 324       }
324     }                                             325     }
325     l.unlock();                                << 326     delete dynParticle;
326   }                                               327   }
327   // prepare isotope selection                 << 
328   std::size_t nIso = temp.size();              << 
329   for ( auto const & elm : *table ) {          << 
330     std::size_t n = elm->GetNumberOfIsotopes() << 
331     if (n > nIso) { nIso = n; }                << 
332   }                                            << 
333   temp.resize(nIso, 0.0);                      << 
334 }                                                 328 }
335                                                   329 
336 const G4String& G4NeutronInelasticXS::FindDire << 330 void 
                                                   >> 331 G4NeutronInelasticXS::Initialise(G4int Z, G4DynamicParticle* dp, 
                                                   >> 332          const char* p)
337 {                                                 333 {
338   // build the complete string identifying the << 334   if(data->GetElementData(Z) || Z < 1 || Z >= MAXZINEL) { return; }
339   if (gDataDirectory.empty()) {                << 335   const char* path = p;
340     std::ostringstream ost;                    << 336   if(!p) {
341     ost << G4HadronicParameters::Instance()->G << 337     // check environment variable 
342     gDataDirectory = ost.str();                << 338     // Build the complete string identifying the file with the data set
                                                   >> 339     path = getenv("G4NEUTRONXSDATA");
                                                   >> 340     if (!path) {
                                                   >> 341       G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
                                                   >> 342       FatalException,
                                                   >> 343                   "Environment variable G4NEUTRONXSDATA is not defined");
                                                   >> 344       return;
                                                   >> 345     }
                                                   >> 346   }
                                                   >> 347   G4DynamicParticle* dynParticle = dp;
                                                   >> 348   if(!dp) {
                                                   >> 349     dynParticle = 
                                                   >> 350       new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
343   }                                               351   }
344   return gDataDirectory;                       << 
345 }                                              << 
346                                                << 
347 void G4NeutronInelasticXS::InitialiseOnFly(G4i << 
348 {                                              << 
349   G4AutoLock l(&nInelasticXSMutex);            << 
350   Initialise(Z);                               << 
351   l.unlock();                                  << 
352 }                                              << 
353                                                   352 
354 void G4NeutronInelasticXS::Initialise(G4int Z) << 353   G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
355 {                                              << 
356   if (nullptr != data->GetElementData(Z)) { re << 
357                                                   354 
358   // upload element data                          355   // upload element data 
359   std::ostringstream ost;                         356   std::ostringstream ost;
360   ost << FindDirectoryPath() << Z;             << 357   ost << path << "/inelast" << Z ;
361   G4PhysicsVector* v = RetrieveVector(ost, tru    358   G4PhysicsVector* v = RetrieveVector(ost, true);
362   data->InitialiseForElement(Z, v);               359   data->InitialiseForElement(Z, v);
363   if (verboseLevel > 1) {                      << 360   /*
364     G4cout  << "G4NeutronInelasticXS::Initiali << 361   G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z 
365             << " A= " << aeff[Z] << "  Amin= " << 362    << " A= " << Amean << "  Amin= " << amin[Z] 
366             << "  Amax= " << amax[Z] << G4endl << 363    << "  Amax= " << amax[Z] << G4endl;
367   }                                            << 364   */
368   // upload isotope data                          365   // upload isotope data
369   G4bool noComp = true;                        << 366   if(amin[Z] > 0) {
370   if (amin[Z] < amax[Z]) {                     << 367     size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
                                                   >> 368     data->InitialiseForComponent(Z, nmax);
371                                                   369 
372     for (G4int A=amin[Z]; A<=amax[Z]; ++A) {   << 370     for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
373       std::ostringstream ost1;                    371       std::ostringstream ost1;
374       ost1 << gDataDirectory << Z << "_" << A; << 372       ost1 << path << "/inelast" << Z << "_" << A;
375       G4PhysicsVector* v1 = RetrieveVector(ost    373       G4PhysicsVector* v1 = RetrieveVector(ost1, false);
376       if (nullptr != v1) {                     << 374       data->AddComponent(Z, A, v1); 
377         if (noComp) {                          << 
378           G4int nmax = amax[Z] - A + 1;        << 
379           data->InitialiseForComponent(Z, nmax << 
380           noComp = false;                      << 
381         }                                      << 
382         data->AddComponent(Z, A, v1);          << 
383       }                                        << 
384     }                                             375     }
385   }                                               376   }
386   // no components case                        << 
387   if (noComp) { data->InitialiseForComponent(Z << 
388                                                   377 
389   // smooth transition                            378   // smooth transition 
390   G4double sig1 = (*v)[v->GetVectorLength()-1] << 379   G4double emax = v->GetMaxEnergy();
391   G4double ehigh= v->GetMaxEnergy();           << 380   G4double sig1 = (*v)[v->GetVectorLength() - 1];
392   G4double sig2 = ggXsection->GetInelasticElem << 381   dynParticle->SetKineticEnergy(emax);
393                               ehigh, Z, aeff[Z << 382   G4double sig2 = 0.0;
394   coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0;    << 383   if(1 == Z) {
                                                   >> 384     fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
                                                   >> 385     sig2 = fNucleon->GetInelasticHadronNucleonXsc();
                                                   >> 386   } else {
                                                   >> 387     ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
                                                   >> 388     sig2 = ggXsection->GetInelasticGlauberGribovXsc();
                                                   >> 389   }
                                                   >> 390   if(sig2 > 0.) { coeff[Z] = sig1/sig2; } 
                                                   >> 391   if(!dp) { delete dynParticle; }
395 }                                                 392 }
396                                                   393 
397 G4PhysicsVector*                                  394 G4PhysicsVector* 
398 G4NeutronInelasticXS::RetrieveVector(std::ostr    395 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
399 {                                                 396 {
400   G4PhysicsLogVector* v = nullptr;             << 397   G4PhysicsLogVector* v = 0;
401   std::ifstream filein(ost.str().c_str());        398   std::ifstream filein(ost.str().c_str());
402   if (!filein.is_open()) {                     << 399   if (!(filein)) {
403     if(warn) {                                    400     if(warn) { 
404       G4ExceptionDescription ed;                  401       G4ExceptionDescription ed;
405       ed << "Data file <" << ost.str().c_str()    402       ed << "Data file <" << ost.str().c_str()
406          << "> is not opened!";                << 403    << "> is not opened!";
407       G4Exception("G4NeutronInelasticXS::Retri    404       G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
408                   FatalException, ed, "Check G << 405       FatalException, ed, "Check G4NEUTRONXSDATA");
409     }                                             406     }
410   } else {                                        407   } else {
411     if(verboseLevel > 1) {                        408     if(verboseLevel > 1) {
412       G4cout << "File " << ost.str()              409       G4cout << "File " << ost.str() 
413              << " is opened by G4NeutronInelas << 410        << " is opened by G4NeutronInelasticXS" << G4endl;
414     }                                             411     }
415     // retrieve data from DB                      412     // retrieve data from DB
416     v = new G4PhysicsLogVector();                 413     v = new G4PhysicsLogVector();
417     if(!v->Retrieve(filein, true)) {              414     if(!v->Retrieve(filein, true)) {
418       G4ExceptionDescription ed;                  415       G4ExceptionDescription ed;
419       ed << "Data file <" << ost.str().c_str()    416       ed << "Data file <" << ost.str().c_str()
420          << "> is not retrieved!";             << 417    << "> is not retrieved!";
421       G4Exception("G4NeutronInelasticXS::Retri    418       G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
422                   FatalException, ed, "Check G << 419       FatalException, ed, "Check G4NEUTRONXSDATA");
423     }                                             420     }
424   }                                               421   }
425   return v;                                       422   return v;
426 }                                                 423 }
427                                                   424