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 9.2.p3)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // -------------------------------------------    
 27 //                                                
 28 // GEANT4 Class file                              
 29 //                                                
 30 //                                                
 31 // File name:    G4NeutronInelasticXS             
 32 //                                                
 33 // Author  Ivantchenko, Geant4, 3-Aug-09          
 34 //                                                
 35                                                   
 36 #include "G4NeutronInelasticXS.hh"                
 37 #include "G4Neutron.hh"                           
 38 #include "G4DynamicParticle.hh"                   
 39 #include "G4ElementTable.hh"                      
 40 #include "G4Material.hh"                          
 41 #include "G4Element.hh"                           
 42 #include "G4PhysicsLogVector.hh"                  
 43 #include "G4CrossSectionDataSetRegistry.hh"       
 44 #include "G4ComponentGGHadronNucleusXsc.hh"       
 45 #include "G4HadronicParameters.hh"                
 46 #include "Randomize.hh"                           
 47 #include "G4Neutron.hh"                           
 48 #include "G4SystemOfUnits.hh"                     
 49 #include "G4IsotopeList.hh"                       
 50 #include "G4NuclearRadii.hh"                      
 51 #include "G4AutoLock.hh"                          
 52                                                   
 53 #include <fstream>                                
 54 #include <sstream>                                
 55 #include <thread>                                 
 56                                                   
 57 G4double G4NeutronInelasticXS::coeff[] = {1.0}    
 58 G4double G4NeutronInelasticXS::lowcoeff[] = {1    
 59 G4ElementData* G4NeutronInelasticXS::data = nu    
 60 G4String G4NeutronInelasticXS::gDataDirectory     
 61                                                   
 62 static std::once_flag applyOnce;                  
 63                                                   
 64 namespace                                         
 65 {                                                 
 66   G4Mutex nInelasticXSMutex = G4MUTEX_INITIALI    
 67 }                                                 
 68                                                   
 69 G4NeutronInelasticXS::G4NeutronInelasticXS()      
 70   : G4VCrossSectionDataSet(Default_Name()),       
 71     neutron(G4Neutron::Neutron()),                
 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 }                                                 
 93                                                   
 94 void G4NeutronInelasticXS::CrossSectionDescrip    
 95 {                                                 
 96   outFile << "G4NeutronInelasticXS calculates     
 97           << "cross section on nuclei using da    
 98           << "neutron database.  These data ar    
 99           << "the resonance region in order to    
100           << "For high energy Glauber-Gribov c    
101 }                                                 
102                                                   
103 G4bool                                            
104 G4NeutronInelasticXS::IsElementApplicable(cons    
105                                           G4in    
106 {                                                 
107   return true;                                    
108 }                                                 
109                                                   
110 G4bool                                            
111 G4NeutronInelasticXS::IsIsoApplicable(const G4    
112                                       G4int, G    
113                                       const G4    
114 {                                                 
115   return true;                                    
116 }                                                 
117                                                   
118 G4double                                          
119 G4NeutronInelasticXS::GetElementCrossSection(c    
120                                              G    
121 {                                                 
122   return ElementCrossSection(aParticle->GetKin    
123                              aParticle->GetLog    
124 }                                                 
125                                                   
126 G4double                                          
127 G4NeutronInelasticXS::ComputeCrossSectionPerEl    
128                                                   
129                                                   
130                                                   
131 {                                                 
132   return ElementCrossSection(ekin, loge, elm->    
133 }                                                 
134                                                   
135 G4double                                          
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                                                   
143   // very low energy limit                        
144   if (ekin < lowElimit) {                         
145     ekin = lowElimit;                             
146     loge = loglowElimit;                          
147   }                                               
148   // pv should exist                              
149   auto pv = GetPhysicsVector(Z);                  
150                                                   
151   const G4double e0 = pv->Energy(0);              
152   if (ekin <= e0) {                               
153     xs = (*pv)[0];                                
154     if (xs > 0.0) { xs *= std::sqrt(e0/ekin);     
155   } else if (ekin <= pv->GetMaxEnergy()) {        
156     xs = pv->LogVectorValue(ekin, loge);          
157   } else {                                        
158     xs = coeff[Z]*ggXsection->GetInelasticElem    
159                                                   
160   }                                               
161                                                   
162 #ifdef G4VERBOSE                                  
163   if(verboseLevel > 1) {                          
164     G4cout  << "G4NeutronInelasticXS::ElementC    
165             << " Ekin(MeV)= " << ekin/CLHEP::M    
166             << ", ElmXSinel(b)= " << xs/CLHEP:    
167             << G4endl;                            
168   }                                               
169 #endif                                            
170   return xs;                                      
171 }                                                 
172                                                   
173 G4double                                          
174 G4NeutronInelasticXS::ComputeIsoCrossSection(G    
175                                              c    
176                                              G    
177                                              c    
178                                              c    
179 {                                                 
180   return IsoCrossSection(ekin, loge, Z, A);       
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                                                   
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   }                                               
233                                                   
234   // use element x-section                        
235   xs = ElementCrossSection(ekin, loge, Z)*A/ae    
236                                                   
237 #ifdef G4VERBOSE                                  
238   if(verboseLevel > 1) {                          
239     G4cout  << "G4NeutronInelasticXS::IsoXS: Z    
240             << " Ekin(MeV)= " << ekin/CLHEP::M    
241             << ", ElmXS(b)= " << xs/CLHEP::bar    
242   }                                               
243 #endif                                            
244   return xs;                                      
245 }                                                 
246                                                   
247 const G4Isotope* G4NeutronInelasticXS::SelectI    
248       const G4Element* anElement, G4double kin    
249 {                                                 
250   std::size_t nIso = anElement->GetNumberOfIso    
251   const G4Isotope* iso = anElement->GetIsotope    
252   if(1 == nIso) { return iso; }                   
253                                                   
254   // more than 1 isotope                          
255   G4int Z = anElement->GetZasInt();               
256   if (nullptr == data->GetElementData(Z)) { In    
257                                                   
258   const G4double* abundVector = anElement->Get    
259   G4double q = G4UniformRand();                   
260   G4double sum = 0.0;                             
261   std::size_t j;                                  
262                                                   
263   // isotope wise cross section not available     
264   if (Z >= MAXZINEL || 0 == data->GetNumberOfC    
265     for (j=0; j<nIso; ++j) {                      
266       sum += abundVector[j];                      
267       if(q <= sum) {                              
268         iso = anElement->GetIsotope((G4int)j);    
269         break;                                    
270       }                                           
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     }                                             
292   }                                               
293   return iso;                                     
294 }                                                 
295                                                   
296 void                                              
297 G4NeutronInelasticXS::BuildPhysicsTable(const     
298 {                                                 
299   if (verboseLevel > 0) {                         
300     G4cout << "G4NeutronInelasticXS::BuildPhys    
301            << p.GetParticleName() << G4endl;      
302   }                                               
303   if (p.GetParticleName() != "neutron") {         
304     G4ExceptionDescription ed;                    
305     ed << p.GetParticleName() << " is a wrong     
306        << " only neutron is allowed";             
307     G4Exception("G4NeutronInelasticXS::BuildPh    
308                 FatalException, ed, "");          
309     return;                                       
310   }                                               
311   // it is possible re-initialisation for the     
312   const G4ElementTable* table = G4Element::Get    
313                                                   
314   // initialise static tables only once           
315   std::call_once(applyOnce, [this]() { isIniti    
316                                                   
317   if (isInitializer) {                            
318     G4AutoLock l(&nInelasticXSMutex);             
319                                                   
320     // Upload data for elements used in geomet    
321     for ( auto const & elm : *table ) {           
322       G4int Z = std::max( 1, std::min( elm->Ge    
323       if ( nullptr == data->GetElementData(Z)     
324     }                                             
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   }                                               
333   temp.resize(nIso, 0.0);                         
334 }                                                 
335                                                   
336 const G4String& G4NeutronInelasticXS::FindDire    
337 {                                                 
338   // build the complete string identifying the    
339   if (gDataDirectory.empty()) {                   
340     std::ostringstream ost;                       
341     ost << G4HadronicParameters::Instance()->G    
342     gDataDirectory = ost.str();                   
343   }                                               
344   return gDataDirectory;                          
345 }                                                 
346                                                   
347 void G4NeutronInelasticXS::InitialiseOnFly(G4i    
348 {                                                 
349   G4AutoLock l(&nInelasticXSMutex);               
350   Initialise(Z);                                  
351   l.unlock();                                     
352 }                                                 
353                                                   
354 void G4NeutronInelasticXS::Initialise(G4int Z)    
355 {                                                 
356   if (nullptr != data->GetElementData(Z)) { re    
357                                                   
358   // upload element data                          
359   std::ostringstream ost;                         
360   ost << FindDirectoryPath() << Z;                
361   G4PhysicsVector* v = RetrieveVector(ost, tru    
362   data->InitialiseForElement(Z, v);               
363   if (verboseLevel > 1) {                         
364     G4cout  << "G4NeutronInelasticXS::Initiali    
365             << " A= " << aeff[Z] << "  Amin= "    
366             << "  Amax= " << amax[Z] << G4endl    
367   }                                               
368   // upload isotope data                          
369   G4bool noComp = true;                           
370   if (amin[Z] < amax[Z]) {                        
371                                                   
372     for (G4int A=amin[Z]; A<=amax[Z]; ++A) {      
373       std::ostringstream ost1;                    
374       ost1 << gDataDirectory << Z << "_" << A;    
375       G4PhysicsVector* v1 = RetrieveVector(ost    
376       if (nullptr != 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     }                                             
385   }                                               
386   // no components case                           
387   if (noComp) { data->InitialiseForComponent(Z    
388                                                   
389   // smooth transition                            
390   G4double sig1 = (*v)[v->GetVectorLength()-1]    
391   G4double ehigh= v->GetMaxEnergy();              
392   G4double sig2 = ggXsection->GetInelasticElem    
393                               ehigh, Z, aeff[Z    
394   coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0;       
395 }                                                 
396                                                   
397 G4PhysicsVector*                                  
398 G4NeutronInelasticXS::RetrieveVector(std::ostr    
399 {                                                 
400   G4PhysicsLogVector* v = nullptr;                
401   std::ifstream filein(ost.str().c_str());        
402   if (!filein.is_open()) {                        
403     if(warn) {                                    
404       G4ExceptionDescription ed;                  
405       ed << "Data file <" << ost.str().c_str()    
406          << "> is not opened!";                   
407       G4Exception("G4NeutronInelasticXS::Retri    
408                   FatalException, ed, "Check G    
409     }                                             
410   } else {                                        
411     if(verboseLevel > 1) {                        
412       G4cout << "File " << ost.str()              
413              << " is opened by G4NeutronInelas    
414     }                                             
415     // retrieve data from DB                      
416     v = new G4PhysicsLogVector();                 
417     if(!v->Retrieve(filein, true)) {              
418       G4ExceptionDescription ed;                  
419       ed << "Data file <" << ost.str().c_str()    
420          << "> is not retrieved!";                
421       G4Exception("G4NeutronInelasticXS::Retri    
422                   FatalException, ed, "Check G    
423     }                                             
424   }                                               
425   return v;                                       
426 }                                                 
427