Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4NeutronCaptureXS.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/G4NeutronCaptureXS.cc (Version 11.3.0) and /processes/hadronic/cross_sections/src/G4NeutronCaptureXS.cc (Version 9.6.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$
                                                   >>  27 //
 26 // -------------------------------------------     28 // -------------------------------------------------------------------
 27 //                                                 29 //
 28 // GEANT4 Class file                               30 // GEANT4 Class file
 29 //                                                 31 //
 30 //                                                 32 //
 31 // File name:    G4NeutronCaptureXS                33 // File name:    G4NeutronCaptureXS
 32 //                                                 34 //
 33 // Author  Ivantchenko, Geant4, 3-Aug-09           35 // Author  Ivantchenko, Geant4, 3-Aug-09
 34 //                                                 36 //
 35 // Modifications:                                  37 // Modifications:
 36 //                                                 38 //
 37                                                    39 
 38 #include <fstream>                                 40 #include <fstream>
 39 #include <sstream>                                 41 #include <sstream>
 40 #include <thread>                              << 
 41                                                    42 
                                                   >>  43 #include "G4HadronicException.hh"
 42 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 43 #include "G4NeutronCaptureXS.hh"                   45 #include "G4NeutronCaptureXS.hh"
 44 #include "G4Material.hh"                       << 
 45 #include "G4Element.hh"                            46 #include "G4Element.hh"
                                                   >>  47 #include "G4ElementTable.hh"
 46 #include "G4PhysicsLogVector.hh"                   48 #include "G4PhysicsLogVector.hh"
                                                   >>  49 #include "G4PhysicsVector.hh"
 47 #include "G4DynamicParticle.hh"                    50 #include "G4DynamicParticle.hh"
 48 #include "G4ElementTable.hh"                   << 
 49 #include "G4IsotopeList.hh"                    << 
 50 #include "G4HadronicParameters.hh"             << 
 51 #include "Randomize.hh"                            51 #include "Randomize.hh"
 52 #include "G4Log.hh"                            << 
 53 #include "G4AutoLock.hh"                       << 
 54                                                    52 
 55 G4ElementData* G4NeutronCaptureXS::data = null <<  53 using namespace std;
 56 G4String G4NeutronCaptureXS::gDataDirectory =  << 
 57                                                    54 
 58 static std::once_flag applyOnce;               <<  55 const G4int G4NeutronCaptureXS::amin[] = {0,
 59                                                <<  56             0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
 60 namespace                                      <<  57             0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
 61 {                                              <<  58             0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
 62   G4Mutex neutronCaptureXSMutex = G4MUTEX_INIT <<  59             0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
 63   const G4int MAXZCAPTURE = 92;                <<  60             0, 0, 0, 0, 0, 0,107,106, 0,112, //41-50
 64 }                                              <<  61             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
                                                   >>  62             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
                                                   >>  63             0, 0, 0,180, 0, 0, 0, 0, 0, 0, //71-80
                                                   >>  64             0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
                                                   >>  65             0,235};
                                                   >>  66 const G4int G4NeutronCaptureXS::amax[] = {0,
                                                   >>  67             0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
                                                   >>  68             0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
                                                   >>  69             0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
                                                   >>  70             0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
                                                   >>  71             0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
                                                   >>  72             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
                                                   >>  73             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
                                                   >>  74             0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
                                                   >>  75             0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
                                                   >>  76             0,238};
 65                                                    77 
 66 G4NeutronCaptureXS::G4NeutronCaptureXS()           78 G4NeutronCaptureXS::G4NeutronCaptureXS() 
 67  : G4VCrossSectionDataSet(Default_Name()),     <<  79  : G4VCrossSectionDataSet("G4NeutronCaptureXS"),
 68    emax(20*CLHEP::MeV), elimit(1.0e-5*CLHEP::e <<  80    emax(20*MeV),maxZ(92)
 69 {                                                  81 {
 70   verboseLevel = 0;                            <<  82   //  verboseLevel = 0;
 71   if (verboseLevel > 0) {                      <<  83   if(verboseLevel > 0){
 72     G4cout  << "G4NeutronCaptureXS::G4NeutronC     84     G4cout  << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
 73       << MAXZCAPTURE << G4endl;                <<  85       << maxZ + 1 << G4endl;
 74   }                                                86   }
 75   logElimit = G4Log(elimit);                   <<  87   //data.resize(maxZ+1, 0);
 76   if (nullptr == data) {                       <<  88   data.SetName("NeutronCapture");
 77     data = new G4ElementData(MAXZCAPTURE+1);   <<  89   work.resize(13,0);
 78     data->SetName("nCapture");                 <<  90   temp.resize(13,0.0);
 79     FindDirectoryPath();                       <<  91   isInitialized = false;
                                                   >>  92 }
                                                   >>  93 
                                                   >>  94 G4NeutronCaptureXS::~G4NeutronCaptureXS()
                                                   >>  95 {
                                                   >>  96   /*
                                                   >>  97   for(G4int i=0; i<=maxZ; ++i) {
                                                   >>  98     delete data[i];
 80   }                                                99   }
                                                   >> 100   */
 81 }                                                 101 }
 82                                                   102 
 83 void G4NeutronCaptureXS::CrossSectionDescripti    103 void G4NeutronCaptureXS::CrossSectionDescription(std::ostream& outFile) const
 84 {                                                 104 {
 85   outFile << "G4NeutronCaptureXS calculates th    105   outFile << "G4NeutronCaptureXS calculates the neutron capture cross sections\n"
 86           << "on nuclei using data from the hi    106           << "on nuclei using data from the high precision neutron database.\n"
 87           << "These data are simplified and sm    107           << "These data are simplified and smoothed over the resonance region\n"
 88           << "in order to reduce CPU time. G4N << 108           << "in order to reduce CPU time.  G4NeutronCaptureXS is valid up to\n"
 89           << "above 20 MeV for all targets. Fo << 109           << "20 MeV for all targets through U.\n";
 90     << "Uranium is used.\n";                   << 
 91 }                                                 110 }
 92                                                   111  
 93 G4bool                                            112 G4bool 
 94 G4NeutronCaptureXS::IsElementApplicable(const     113 G4NeutronCaptureXS::IsElementApplicable(const G4DynamicParticle*, 
 95           G4int, const G4Material*)               114           G4int, const G4Material*)
 96 {                                                 115 {
 97   return true;                                    116   return true;
 98 }                                                 117 }
 99                                                   118 
100 G4bool                                            119 G4bool 
101 G4NeutronCaptureXS::IsIsoApplicable(const G4Dy    120 G4NeutronCaptureXS::IsIsoApplicable(const G4DynamicParticle*,
102             G4int, G4int,                      << 121             G4int /*ZZ*/, G4int /*AA*/,
103             const G4Element*, const G4Material    122             const G4Element*, const G4Material*)
104 {                                                 123 {
105   return true;                                 << 124   return false;
106 }                                                 125 }
107                                                   126 
108 G4double                                          127 G4double 
109 G4NeutronCaptureXS::GetElementCrossSection(con    128 G4NeutronCaptureXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
110              G4int Z, const G4Material*)          129              G4int Z, const G4Material*)
111 {                                                 130 {
112   G4double xs = 0.0;                              131   G4double xs = 0.0;
113   G4double ekin = aParticle->GetKineticEnergy(    132   G4double ekin = aParticle->GetKineticEnergy();
114   if (ekin < emax) {                           << 133   if(ekin > emax || Z < 1 || Z > maxZ) { return xs; }
115     xs = ElementCrossSection(ekin, aParticle-> << 134   const G4double elimit = 1.0e-10*eV;
116   }                                            << 135   if(ekin < elimit) { ekin = elimit; }
117   return xs;                                   << 136 
118 }                                              << 137   //  G4PhysicsVector* pv = data[Z];
                                                   >> 138   G4PhysicsVector* pv = data.GetElementData(Z);
                                                   >> 139 
                                                   >> 140   // element was not initialised
                                                   >> 141   if(!pv) {
                                                   >> 142     Initialise(Z);
                                                   >> 143     //    pv = data[Z];
                                                   >> 144     pv = data.GetElementData(Z);
                                                   >> 145     if(!pv) { return xs; }
                                                   >> 146   }
                                                   >> 147 
                                                   >> 148   G4double e1 = pv->Energy(0);
                                                   >> 149   if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); }
                                                   >> 150   else          { xs = pv->Value(ekin); }
119                                                   151 
120 G4double                                       << 152   if(verboseLevel > 0){
121 G4NeutronCaptureXS::ComputeCrossSectionPerElem << 153     G4cout  << "ekin= " << ekin << ",  xs= " << xs << G4endl;
122                           const G4ParticleDefi << 
123                           const G4Element* elm << 
124                           const G4Material*)   << 
125 {                                              << 
126   G4double xs = 0.0;                           << 
127   if (ekin < emax) {                           << 
128     xs = ElementCrossSection(ekin, loge, elm-> << 
129   }                                               154   }
130   return xs;                                      155   return xs;
131 }                                                 156 }
132                                                   157 
133 G4double                                       << 
134 G4NeutronCaptureXS::ElementCrossSection(G4doub << 
135 {                                              << 
136   G4int Z = std::min(ZZ, MAXZCAPTURE);         << 
137   G4double ekin = eKin;                        << 
138   G4double logEkin = logE;                     << 
139   if (ekin < elimit) {                         << 
140     ekin = elimit;                             << 
141     logEkin = logElimit;                       << 
142   }                                            << 
143                                                << 
144   auto pv = GetPhysicsVector(Z);               << 
145   const G4double e0 = pv->Energy(0);           << 
146   G4double xs = (ekin >= e0) ? pv->LogVectorVa << 
147     : (*pv)[0]*std::sqrt(e0/ekin);             << 
148                                                << 
149 #ifdef G4VERBOSE                               << 
150   if (verboseLevel > 1){                       << 
151     G4cout  << "Ekin= " << ekin/CLHEP::MeV     << 
152             << " ElmXScap(b)= " << xs/CLHEP::b << 
153   }                                            << 
154 #endif                                         << 
155   return xs;                                   << 
156 }                                              << 
157                                                << 
158 G4double                                       << 
159 G4NeutronCaptureXS::ComputeIsoCrossSection(G4d << 
160                    const G4ParticleDefinition* << 
161                    G4int Z, G4int A,           << 
162                    const G4Isotope*, const G4E << 
163                    const G4Material*)          << 
164 {                                              << 
165   return IsoCrossSection(ekin, loge, Z, A);    << 
166 }                                              << 
167                                                << 
168 G4double                                          158 G4double 
169 G4NeutronCaptureXS::GetIsoCrossSection(const G    159 G4NeutronCaptureXS::GetIsoCrossSection(const G4DynamicParticle* aParticle, 
170                G4int Z, G4int A,                  160                G4int Z, G4int A,
171                const G4Isotope*, const G4Eleme    161                const G4Isotope*, const G4Element*,
172                const G4Material*)                 162                const G4Material*)
173 {                                                 163 {
174   return IsoCrossSection(aParticle->GetKinetic << 
175                          aParticle->GetLogKine << 
176                          Z, A);                << 
177 }                                              << 
178                                                << 
179 G4double G4NeutronCaptureXS::IsoCrossSection(G << 
180                                              G << 
181 {                                              << 
182   G4double xs = 0.0;                              164   G4double xs = 0.0;
183   if (eKin > emax) { return xs; }              << 165   G4double ekin = aParticle->GetKineticEnergy();
                                                   >> 166   if(ekin > emax || Z < 1 || Z > maxZ) { return xs; }
                                                   >> 167   const G4double elimit = 1.0e-10*eV;
                                                   >> 168   if(ekin < elimit) { ekin = elimit; }
                                                   >> 169 
                                                   >> 170   //  G4PhysicsVector* pv = data[Z];
                                                   >> 171   G4PhysicsVector* pv = data.GetElementData(Z);
                                                   >> 172 
                                                   >> 173   // element was not initialised
                                                   >> 174   if(!pv) {
                                                   >> 175     Initialise(Z);
                                                   >> 176     //    pv = data[Z];
                                                   >> 177     pv = data.GetElementData(Z);
                                                   >> 178     if(!pv) { return xs; }
                                                   >> 179   }
                                                   >> 180   pv = data.GetComponentDataByID(Z, A);
                                                   >> 181   if(!pv) { return xs; }
                                                   >> 182 
                                                   >> 183   G4double e1 = pv->Energy(0);
                                                   >> 184   if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); }
                                                   >> 185   else          { xs = pv->Value(ekin); }
184                                                   186 
185   G4int Z = std::min(ZZ, MAXZCAPTURE);         << 187   if(verboseLevel > 0){
186   G4double ekin = eKin;                        << 188     G4cout  << "ekin= " << ekin << ",  xs= " << xs << G4endl;
187   G4double logEkin = logE;                     << 
188   if (ekin < elimit) {                         << 
189     ekin = elimit;                             << 
190     logEkin = logElimit;                       << 
191   }                                            << 
192                                                << 
193   auto pv = GetPhysicsVector(Z);               << 
194   if (pv == nullptr) { return xs; }            << 
195                                                << 
196   // use isotope x-section if possible         << 
197   if (data->GetNumberOfComponents(Z) > 0) {    << 
198     G4PhysicsVector* pviso = data->GetComponen << 
199     if(pviso != nullptr) {                     << 
200       const G4double e0 = pviso->Energy(0);    << 
201       xs = (ekin >= e0) ? pviso->LogVectorValu << 
202   : (*pviso)[0]*std::sqrt(e0/ekin);            << 
203 #ifdef G4VERBOSE                               << 
204       if(verboseLevel > 0) {                   << 
205   G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(M << 
206          << "  xs(b)= " << xs/barn             << 
207          << "  Z= " << Z << "  A= " << A << G4 << 
208       }                                        << 
209 #endif                                         << 
210       return xs;                               << 
211     }                                          << 
212   }                                               189   }
213   // isotope data are not available or applica << 
214   const G4double e0 = pv->Energy(0);           << 
215   xs = (ekin >= e0) ? pv->LogVectorValue(ekin, << 
216     : (*pv)[0]*std::sqrt(e0/ekin);             << 
217 #ifdef G4VERBOSE                               << 
218   if (verboseLevel > 0) {                      << 
219     G4cout << "G4NeutronCaptureXS::IsoXS: Ekin << 
220            << "  xs(b)= " << xs/barn           << 
221      << "  Z= " << Z << "  A= " << A << " no i << 
222   }                                            << 
223 #endif                                         << 
224   return xs;                                      190   return xs;
225 }                                                 191 }
226                                                   192 
227 const G4Isotope*                               << 193 G4Isotope* G4NeutronCaptureXS::SelectIsotope(const G4Element* anElement,
228 G4NeutronCaptureXS::SelectIsotope(const G4Elem << 194                G4double kinEnergy)
229           G4double kinEnergy, G4double logE)   << 
230 {                                                 195 {
231   G4int nIso = (G4int)anElement->GetNumberOfIs << 196   G4int nIso = anElement->GetNumberOfIsotopes();
232   const G4Isotope* iso = anElement->GetIsotope << 197   G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
233                                                << 198   G4Isotope* iso = (*isoVector)[0];
234   //G4cout << "SelectIsotope NIso= " << nIso < << 
235   if(1 == nIso) { return iso; }                << 
236                                                   199 
237   // more than 1 isotope                          200   // more than 1 isotope
238   G4int Z = anElement->GetZasInt();            << 201   if(1 < nIso) {
239   if (nullptr == data->GetElementData(Z)) { In << 202     G4int Z = G4lrint(anElement->GetZ());
240                                                << 203     if(Z > maxZ) { Z = maxZ; }
241   const G4double* abundVector = anElement->Get << 204     G4double* abundVector = anElement->GetRelativeAbundanceVector();
242   G4double q = G4UniformRand();                << 205     G4double q = G4UniformRand();
243   G4double sum = 0.0;                          << 206     G4double sum = 0.0;
244                                                << 207 
245   // is there isotope wise cross section?      << 208     // is there isotope wise cross section?
246   G4int j;                                     << 209     if(0 == amin[Z]) {
247   if (Z > MAXZCAPTURE || 0 == data->GetNumberO << 210       for (G4int j = 0; j<nIso; ++j) {
248     for (j = 0; j<nIso; ++j) {                 << 211   sum += abundVector[j];
249       sum += abundVector[j];                   << 212   if(q <= sum) {
250       if(q <= sum) {                           << 213     iso = (*isoVector)[j];
251   iso = anElement->GetIsotope(j);              << 214     break;
252   break;                                       << 215   }
                                                   >> 216       }
                                                   >> 217     } else {
                                                   >> 218       size_t nmax = data.GetNumberOfComponents(Z);
                                                   >> 219       if(temp.size() < nmax) { temp.resize(nmax,0.0); }
                                                   >> 220       for (size_t i=0; i<nmax; ++i) {
                                                   >> 221   G4int A = (*isoVector)[i]->GetN();
                                                   >> 222   G4PhysicsVector* v = data.GetComponentDataByID(Z, A);
                                                   >> 223         if(v) { sum += abundVector[i]*v->Value(kinEnergy); }
                                                   >> 224         temp[i] = sum;
                                                   >> 225       }
                                                   >> 226       sum *= q;
                                                   >> 227       for (size_t j = 0; j<nmax; ++j) {
                                                   >> 228         if(temp[j] >= sum) {
                                                   >> 229           iso = (*isoVector)[j];
                                                   >> 230           break;
                                                   >> 231   }
253       }                                           232       }
254     }                                          << 
255     return iso;                                << 
256   }                                            << 
257   G4int nn = (G4int)temp.size();               << 
258   if (nn < nIso) { temp.resize(nIso, 0.); }    << 
259                                                << 
260   for (j=0; j<nIso; ++j) {                     << 
261     sum += abundVector[j]*IsoCrossSection(kinE << 
262             anElement->GetIsotope(j)->GetN()); << 
263     temp[j] = sum;                             << 
264   }                                            << 
265   sum *= q;                                    << 
266   for (j = 0; j<nIso; ++j) {                   << 
267     if (temp[j] >= sum) {                      << 
268       iso = anElement->GetIsotope(j);          << 
269       break;                                   << 
270     }                                             233     }
271   }                                               234   }
272   return iso;                                     235   return iso;
273 }                                                 236 }
274                                                   237 
275 void                                              238 void 
276 G4NeutronCaptureXS::BuildPhysicsTable(const G4    239 G4NeutronCaptureXS::BuildPhysicsTable(const G4ParticleDefinition& p)
277 {                                                 240 {
278   if (verboseLevel > 0){                       << 241   if(isInitialized) { return; }
                                                   >> 242   if(verboseLevel > 0){
279     G4cout << "G4NeutronCaptureXS::BuildPhysic    243     G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for " 
280      << p.GetParticleName() << G4endl;            244      << p.GetParticleName() << G4endl;
281   }                                               245   }
282   if (p.GetParticleName() != "neutron") {      << 246   if(p.GetParticleName() != "neutron") { 
283     G4ExceptionDescription ed;                 << 247     throw G4HadronicException(__FILE__, __LINE__,"Wrong particle type");
284     ed << p.GetParticleName() << " is a wrong  << 
285        << " only neutron is allowed";          << 
286     G4Exception("G4NeutronCaptureXS::BuildPhys << 
287     FatalException, ed, "");                   << 
288     return;                                       248     return; 
289   }                                               249   }
                                                   >> 250   isInitialized = true;
290                                                   251 
291   // it is possible re-initialisation for the  << 252   // check environment variable 
292   const G4ElementTable* table = G4Element::Get << 253   // Build the complete string identifying the file with the data set
293                                                << 254   char* path = getenv("G4NEUTRONXSDATA");
294   // initialise static tables only once        << 255   if (!path){
295   std::call_once(applyOnce, [this]() { isIniti << 256     throw G4HadronicException(__FILE__, __LINE__, 
296                                                << 257             "G4NEUTRONXSDATA environment variable not defined");
297   if (isInitializer) {                         << 258     return;
298     G4AutoLock l(&neutronCaptureXSMutex);      << 259   }
299     // Access to elements                      << 260 
300     for ( auto const & elm : *table ) {        << 261   // Access to elements
301       G4int Z = std::max( 1, std::min( elm->Ge << 262   const G4ElementTable* theElmTable = G4Element::GetElementTable();
302       if ( nullptr == data->GetElementData(Z)  << 263   size_t numOfElm = G4Element::GetNumberOfElements();
                                                   >> 264   if(numOfElm > 0) {
                                                   >> 265     for(size_t i=0; i<numOfElm; ++i) {
                                                   >> 266       G4int Z = G4int(((*theElmTable)[i])->GetZ());
                                                   >> 267       if(Z < 1)         { Z = 1; }
                                                   >> 268       else if(Z > maxZ) { Z = maxZ; }
                                                   >> 269       //G4cout << "Z= " << Z << G4endl;
                                                   >> 270       // Initialisation 
                                                   >> 271       //      if(!data[Z]) { Initialise(Z, path); }
                                                   >> 272       if(!data.GetElementData(Z)) { Initialise(Z, path); }
303     }                                             273     }
304     l.unlock();                                << 
305   }                                            << 
306                                                << 
307   // prepare isotope selection                 << 
308   std::size_t nIso = temp.size();              << 
309   for ( auto const & elm : *table ) {          << 
310     std::size_t n = elm->GetNumberOfIsotopes() << 
311     if (n > nIso) { nIso = n; }                << 
312   }                                               274   }
313   temp.resize(nIso, 0.0);                      << 
314 }                                                 275 }
315                                                   276 
316 const G4String& G4NeutronCaptureXS::FindDirect << 277 void 
317 {                                              << 278 G4NeutronCaptureXS::Initialise(G4int Z, const char* p)
318   // build the complete string identifying the << 
319   if(gDataDirectory.empty()) {                 << 
320     std::ostringstream ost;                    << 
321     ost << G4HadronicParameters::Instance()->G << 
322     gDataDirectory = ost.str();                << 
323   }                                            << 
324   return gDataDirectory;                       << 
325 }                                              << 
326                                                << 
327 void G4NeutronCaptureXS::InitialiseOnFly(G4int << 
328 {                                                 279 {
329   G4AutoLock l(&neutronCaptureXSMutex);        << 280   if(data.GetElementData(Z)) { return; }
330   Initialise(Z);                               << 281   const char* path = p;
331   l.unlock();                                  << 
332 }                                              << 
333                                                   282 
334 void G4NeutronCaptureXS::Initialise(G4int Z)   << 283   // check environment variable 
335 {                                              << 284   if(!p) {
336   if (nullptr != data->GetElementData(Z)) { re << 285     path = getenv("G4NEUTRONXSDATA");
                                                   >> 286     if (!path) {
                                                   >> 287       throw G4HadronicException(__FILE__, __LINE__, 
                                                   >> 288         "G4NEUTRONXSDATA environment variable not defined");
                                                   >> 289       return;
                                                   >> 290     }
                                                   >> 291   }
337                                                   292 
338   // upload element data                          293   // upload element data 
339   std::ostringstream ost;                         294   std::ostringstream ost;
340   ost << FindDirectoryPath() << Z ;            << 295   ost << path << "/cap" << Z ;
341   G4PhysicsVector* v = RetrieveVector(ost, tru    296   G4PhysicsVector* v = RetrieveVector(ost, true);
342   data->InitialiseForElement(Z, v);            << 297   data.InitialiseForElement(Z, v);
343                                                   298 
344   // upload isotope data                          299   // upload isotope data
345   G4bool noComp = true;                        << 300   if(amin[Z] > 0) {
346   if (amin[Z] < amax[Z]) {                     << 301     size_t n = 0;
                                                   >> 302     size_t i = 0;
                                                   >> 303     size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
                                                   >> 304     if(work.size() < nmax) { work.resize(nmax,0); }
347     for(G4int A=amin[Z]; A<=amax[Z]; ++A) {       305     for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
348       std::ostringstream ost1;                    306       std::ostringstream ost1;
349       ost1 << gDataDirectory << Z << "_" << A; << 307       ost1 << path << "/cap" << Z << "_" << A;
350       G4PhysicsVector* v1 = RetrieveVector(ost << 308       v = RetrieveVector(ost1, false);
351       if (nullptr != v1) {                     << 309       if(v) { ++n; }
352   if (noComp) {                                << 310       work[i] = v;
353     G4int nmax = amax[Z] - A + 1;              << 311       ++i;
354     data->InitialiseForComponent(Z, nmax);     << 312     }
355     noComp = false;                            << 313     data.InitialiseForComponent(Z, n);
356   }                                            << 314     for(size_t j=0; j<i; ++j) {
357   data->AddComponent(Z, A, v1);                << 315       if(work[j]) { data.AddComponent(Z, amin[Z]+j, work[j]); }
358       }                                        << 
359     }                                             316     }
360   }                                               317   }
361   // no components case                        << 
362   if (noComp) { data->InitialiseForComponent(Z << 
363 }                                                 318 }
364                                                   319  
365 G4PhysicsVector*                                  320 G4PhysicsVector* 
366 G4NeutronCaptureXS::RetrieveVector(std::ostrin    321 G4NeutronCaptureXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
367 {                                                 322 {
368   G4PhysicsLogVector* v = nullptr;             << 323   G4PhysicsLogVector* v = 0;
369   std::ifstream filein(ost.str().c_str());        324   std::ifstream filein(ost.str().c_str());
370   if (!filein.is_open()) {                     << 325   if (!(filein)) {
371     if (warn) {                                << 326     if(!warn) { return v; }
372       G4ExceptionDescription ed;               << 327     G4cout << ost.str() << "  is not opened by G4NeutronCaptureXS" << G4endl;
373       ed << "Data file <" << ost.str().c_str() << 328     throw G4HadronicException(__FILE__, __LINE__,"NO data sets opened");
374    << "> is not opened!";                      << 329   }else{
375       G4Exception("G4NeutronCaptureXS::Retriev << 330     if(verboseLevel > 1) {
376       FatalException, ed, "Check G4PARTICLEXSD << 
377     }                                          << 
378   } else {                                     << 
379     if (verboseLevel > 1) {                    << 
380       G4cout << "File " << ost.str()              331       G4cout << "File " << ost.str() 
381        << " is opened by G4NeutronCaptureXS" <    332        << " is opened by G4NeutronCaptureXS" << G4endl;
382     }                                             333     }
383     // retrieve data from DB                      334     // retrieve data from DB
384     v = new G4PhysicsLogVector();                 335     v = new G4PhysicsLogVector();
385     if (!v->Retrieve(filein, true)) {          << 336     if(!v->Retrieve(filein, true)) {
386       G4ExceptionDescription ed;               << 337       G4cout << ost.str() << " is not retrieved in G4NeutronCaptureXS" << G4endl;
387       ed << "Data file <" << ost.str().c_str() << 338       throw G4HadronicException(__FILE__, __LINE__,"ERROR: retrieve data fail");
388    << "> is not retrieved!";                   << 
389       G4Exception("G4NeutronCaptureXS::Retriev << 
390       FatalException, ed, "Check G4PARTICLEXSD << 
391     }                                             339     }
392   }                                               340   }
393   return v;                                       341   return v;
394 }                                                 342 }
395                                                   343