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.0.p2)


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