Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4GammaNuclearXS.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // -------------------------------------------------------------------
 27 //
 28 // GEANT4 Class file
 29 //
 30 //
 31 // File name:    G4GammaNuclearXS
 32 //
 33 // Authors  V.Ivantchenko, Geant4, 20 October 2020
 34 //          B.Kutsenko, BINP/NSU, 10 August 2021
 35 //
 36 // Modifications:
 37 //
 38 
 39 #include "G4GammaNuclearXS.hh"
 40 #include "G4DynamicParticle.hh"
 41 #include "G4ThreeVector.hh"
 42 #include "G4ElementTable.hh"
 43 #include "G4Material.hh"
 44 #include "G4Element.hh"
 45 #include "G4ElementData.hh"
 46 #include "G4PhysicsVector.hh"
 47 #include "G4PhysicsLinearVector.hh"
 48 #include "G4PhysicsFreeVector.hh"
 49 #include "G4CrossSectionDataSetRegistry.hh"
 50 #include "G4PhotoNuclearCrossSection.hh"
 51 #include "G4HadronicParameters.hh"
 52 #include "Randomize.hh"
 53 #include "G4SystemOfUnits.hh"
 54 #include "G4Gamma.hh"
 55 #include "G4IsotopeList.hh"
 56 #include "G4AutoLock.hh"
 57 
 58 #include <fstream>
 59 #include <sstream>
 60 #include <vector>
 61 
 62 G4ElementData* G4GammaNuclearXS::data = nullptr;
 63 
 64 G4double G4GammaNuclearXS::coeff[3][3];
 65 G4double G4GammaNuclearXS::xs150[] = {0.0};
 66 const G4int G4GammaNuclearXS::freeVectorException[] = {
 67 4, 6, 7, 8, 27, 39, 45, 65, 67, 69, 73};
 68 G4String G4GammaNuclearXS::gDataDirectory = "";
 69 
 70 namespace
 71 {
 72   // Upper limit of the linear transition between IAEA database and CHIPS model
 73   const G4double eTransitionBound = 150.*CLHEP::MeV;
 74   // A limit energy to correct CHIPS parameterisation for light isotopes 
 75   const G4double ehigh = 10*CLHEP::GeV;
 76 }
 77 
 78 G4GammaNuclearXS::G4GammaNuclearXS() 
 79   : G4VCrossSectionDataSet(Default_Name()), gamma(G4Gamma::Gamma())
 80 {
 81   verboseLevel = 0;
 82   if (verboseLevel > 0) {
 83     G4cout << "G4GammaNuclearXS::G4GammaNuclearXS Initialise for Z < " 
 84      << MAXZGAMMAXS << G4endl;
 85   }
 86   ggXsection = dynamic_cast<G4PhotoNuclearCrossSection*>
 87     (G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet("PhotoNuclearXS"));
 88   if (ggXsection == nullptr) {
 89     ggXsection = new G4PhotoNuclearCrossSection();
 90   }
 91   SetForAllAtomsAndEnergies(true);
 92 
 93   // full data set is uploaded once
 94   if (nullptr == data) { 
 95     data = new G4ElementData(MAXZGAMMAXS);
 96     data->SetName("gNuclear");
 97     for (G4int Z=1; Z<MAXZGAMMAXS; ++Z) {
 98       Initialise(Z);
 99     }
100   }
101 }
102 
103 void G4GammaNuclearXS::CrossSectionDescription(std::ostream& outFile) const
104 {
105   outFile << "G4GammaNuclearXS calculates the gamma nuclear\n"
106           << "cross-section for GDR energy region on nuclei using "
107     << "data from the high precision\n"
108           << "IAEA photonuclear database (2019). Then liniear connection\n"
109     << "implemented with previous CHIPS photonuclear model." << G4endl;
110 }
111 
112 G4bool 
113 G4GammaNuclearXS::IsElementApplicable(const G4DynamicParticle*, 
114                                       G4int, const G4Material*)
115 {
116   return true;
117 }
118 
119 G4bool G4GammaNuclearXS::IsIsoApplicable(const G4DynamicParticle*,
120                                          G4int, G4int,
121                                          const G4Element*, const G4Material*)
122 {
123   return true;
124 }
125 
126 G4double 
127 G4GammaNuclearXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
128                                          G4int Z, const G4Material*)
129 {
130   return ElementCrossSection(aParticle->GetKineticEnergy(), Z); 
131 }
132 
133 G4double
134 G4GammaNuclearXS::ElementCrossSection(const G4double ekin, const G4int ZZ)
135 {
136   // check cache
137   const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
138   if(Z == fZ && ekin == fEkin) { return fXS; }
139   fZ = Z;
140   fEkin = ekin;
141 
142   auto pv = data->GetElementData(Z);
143   const G4double limCHIPS1 = 25*CLHEP::MeV;
144   const G4double limCHIPS2 = 16*CLHEP::MeV;
145   if (pv == nullptr || 1 == Z || Z == 40 || Z == 74 ||
146      (Z == 24 && ekin >= limCHIPS1) ||
147      (Z == 39 && ekin >= limCHIPS1) ||
148      (Z == 50 && ekin >= limCHIPS2) ||
149      (Z == 64 && ekin >= limCHIPS2)
150      ) {
151     fXS = ggXsection->ComputeElementXSection(ekin, Z);
152     return fXS;
153   }
154   const G4double emax = pv->GetMaxEnergy();
155 
156   // low energy based on data
157   if(ekin <= emax) {
158     fXS = pv->Value(ekin);
159     // high energy CHIPS parameterisation
160   } else if(ekin >= eTransitionBound) {
161     fXS = ggXsection->ComputeElementXSection(ekin, Z);
162     // linear interpolation
163   } else {
164     const G4double rxs = xs150[Z];
165     const G4double lxs = pv->Value(emax);
166     fXS = lxs + (ekin - emax)*(rxs - lxs)/(eTransitionBound - emax);
167   }
168 
169 #ifdef G4VERBOSE
170   if(verboseLevel > 1) {
171     G4cout  << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV 
172       << ",  nElmXS(b)= " << fXS/CLHEP::barn 
173       << G4endl;
174   }
175 #endif
176   return fXS;
177 }
178 
179 G4double G4GammaNuclearXS::LowEnergyCrossSection(G4double ekin, G4int ZZ)
180 {    
181   const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
182   auto pv = data->GetElementData(Z);
183   return pv->Value(ekin);
184 }
185 
186 G4double G4GammaNuclearXS::GetIsoCrossSection(
187          const G4DynamicParticle* aParticle,
188    G4int Z, G4int A,
189    const G4Isotope*, const G4Element*, const G4Material*)
190 {
191   return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
192 }
193 
194 G4double 
195 G4GammaNuclearXS::IsoCrossSection(const G4double ekin, const G4int ZZ, const G4int A)
196 {
197   const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
198   // cross section per element
199   G4double xs = ElementCrossSection(ekin, Z);
200 
201   if (Z > 2) {
202     xs *= A/aeff[Z];
203   } else {
204     G4int AA = A - amin[Z];
205     if(ekin >= ehigh && AA >=0 && AA <=2) { 
206       xs *= coeff[Z][AA];
207     } else {
208       xs = ggXsection->ComputeIsoXSection(ekin, Z, A);
209     }
210   }
211 
212 #ifdef G4VERBOSE
213   if(verboseLevel > 1) {
214     G4cout  << "G4GammaNuclearXS::IsoXS: Z= " << Z << " A= " << A 
215       << " Ekin(MeV)= " << ekin/CLHEP::MeV 
216       << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
217   }
218 #endif
219   return xs;
220 }
221 
222 const G4Isotope* G4GammaNuclearXS::SelectIsotope(
223        const G4Element* anElement, G4double kinEnergy, G4double)
224 {
225   std::size_t nIso = anElement->GetNumberOfIsotopes();
226   const G4Isotope* iso = anElement->GetIsotope(0);
227 
228   if(1 == nIso) { return iso; }
229 
230   const G4double* abundVector = anElement->GetRelativeAbundanceVector();
231   G4double sum = 0.0;
232   G4int Z = anElement->GetZasInt();
233 
234   // use isotope cross sections
235   std::size_t nn = temp.size();
236   if(nn < nIso) { temp.resize(nIso, 0.); }
237   
238   for (std::size_t j=0; j<nIso; ++j) {
239     //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN() 
240     //       <<  " abund= " << abundVector[j] << G4endl;
241     sum += abundVector[j]*
242       IsoCrossSection(kinEnergy, Z, anElement->GetIsotope((G4int)j)->GetN());
243     temp[j] = sum;
244   }
245   sum *= G4UniformRand();
246   for (std::size_t j = 0; j<nIso; ++j) {
247     if(temp[j] >= sum) {
248       iso = anElement->GetIsotope((G4int)j);
249       break;
250     }
251   }
252   return iso;
253 }
254 
255 void G4GammaNuclearXS::BuildPhysicsTable(const G4ParticleDefinition& p)
256 {
257   if(verboseLevel > 1) {
258     G4cout << "G4GammaNuclearXS::BuildPhysicsTable for " 
259      << p.GetParticleName() << G4endl;
260   }
261   if(p.GetParticleName() != "gamma") { 
262     G4ExceptionDescription ed;
263     ed << p.GetParticleName() << " is a wrong particle type -"
264        << " only gamma is allowed";
265     G4Exception("G4GammaNuclearXS::BuildPhysicsTable(..)","had012",
266     FatalException, ed, "");
267     return; 
268   }
269 
270   // prepare isotope selection
271   const G4ElementTable* table = G4Element::GetElementTable();
272   std::size_t nIso = temp.size();
273   for (auto const & elm : *table ) {
274     std::size_t n = elm->GetNumberOfIsotopes();
275     if (n > nIso) { nIso = n; }
276   }
277   temp.resize(nIso, 0.0);   
278 }
279 
280 const G4String& G4GammaNuclearXS::FindDirectoryPath()
281 {
282   // build the complete string identifying the file with the data set
283   if(gDataDirectory.empty()) {
284     std::ostringstream ost;
285     ost << G4HadronicParameters::Instance()->GetDirPARTICLEXS() << "/gamma/inel";
286     gDataDirectory = ost.str();
287   }
288   return gDataDirectory;
289 }
290 
291 void G4GammaNuclearXS::Initialise(G4int Z)
292 {
293   // upload data from file
294   std::ostringstream ost;
295   ost << FindDirectoryPath() << Z ;
296   G4PhysicsVector* v = RetrieveVector(ost, true, Z);
297   
298   data->InitialiseForElement(Z, v);
299   /*
300   G4cout << "G4GammaNuclearXS::Initialise for Z= " << Z 
301    << " A= " << Amean << "  Amin= " << amin[Z] 
302    << "  Amax= " << amax[Z] << G4endl;
303   */
304   xs150[Z] = ggXsection->ComputeElementXSection(eTransitionBound, Z);
305 
306   // compute corrections for low Z data
307   if(Z <= 2){
308     if(amax[Z] > amin[Z]) {
309       for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
310         G4int AA = A - amin[Z];
311   if(AA >= 0 && AA <= 2) {
312     G4double sig1 = ggXsection->ComputeIsoXSection(ehigh, Z, A);
313     G4double sig2 = ggXsection->ComputeElementXSection(ehigh, Z);
314     if(sig2 > 0.) { coeff[Z][AA] = (sig1/sig2); }
315     else { coeff[Z][AA] = 1.; }
316   }
317       }
318     }
319   }
320 }
321 
322 G4PhysicsVector* 
323 G4GammaNuclearXS::RetrieveVector(std::ostringstream& ost, G4bool warn, G4int Z)
324 {
325   G4PhysicsVector* v = nullptr;
326 
327   std::ifstream filein(ost.str().c_str());
328   if (!filein.is_open()) {
329     if(warn) {
330       G4ExceptionDescription ed;
331       ed << "Data file <" << ost.str().c_str()
332    << "> is not opened!";
333       G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had014",
334       FatalException, ed, "Check G4PARTICLEXSDATA");
335     }
336   } else {
337     if(verboseLevel > 1) {
338       G4cout << "File " << ost.str() 
339        << " is opened by G4GammaNuclearXS" << G4endl;
340     }
341     // retrieve data from DB
342     if(std::find(std::begin(freeVectorException), std::end(freeVectorException), Z)
343        == std::end(freeVectorException)) {
344       v = new G4PhysicsLinearVector(false);
345     } else {
346       v = new G4PhysicsFreeVector(false);
347     }
348     if(!v->Retrieve(filein, true)) {
349       G4ExceptionDescription ed;
350       ed << "Data file <" << ost.str().c_str()
351    << "> is not retrieved!";
352       G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had015",
353       FatalException, ed, "Check G4PARTICLEXSDATA");
354     }
355   }
356   return v;
357 }
358 
359