Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPJENDLHEData.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 // Class Description
 27 // Cross-section data set for a high precision (based on JENDL_HE evaluated data
 28 // libraries) description of elastic scattering 20 MeV ~ 3 GeV;
 29 // Class Description - End
 30 
 31 // 15-Nov-06 First Implementation is done by T. Koi (SLAC/SCCS)
 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
 33 //
 34 #include "G4ParticleHPJENDLHEData.hh"
 35 
 36 #include "G4ElementTable.hh"
 37 #include "G4ParticleHPData.hh"
 38 #include "G4PhysicsFreeVector.hh"
 39 #include "G4Pow.hh"
 40 #include "G4SystemOfUnits.hh"
 41 
 42 G4bool G4ParticleHPJENDLHEData::IsApplicable(const G4DynamicParticle* aP, const G4Element* anE)
 43 {
 44   G4bool result = true;
 45   G4double eKin = aP->GetKineticEnergy();
 46   // if(eKin>20*MeV||aP->GetDefinition()!=G4Neutron::Neutron()) result = false;
 47   if (eKin < 20 * MeV || 3 * GeV < eKin || aP->GetDefinition() != G4Neutron::Neutron()) {
 48     result = false;
 49   }
 50   // Element Check
 51   else if (!(vElement[anE->GetIndex()]))
 52     result = false;
 53 
 54   return result;
 55 }
 56 
 57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData()
 58 {
 59   for (auto& itZ : mIsotope) {
 60     std::map<G4int, G4PhysicsVector*>* pointer_map = itZ.second;
 61     if (pointer_map != nullptr) {
 62       for (auto& itA : *pointer_map) {
 63         G4PhysicsVector* pointerPhysicsVector = itA.second;
 64         if (pointerPhysicsVector != nullptr) {
 65           delete pointerPhysicsVector;
 66           itA.second = NULL;
 67         }
 68       }
 69       delete pointer_map;
 70       itZ.second = NULL;
 71     }
 72   }
 73   mIsotope.clear();
 74 }
 75 
 76 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEData(const G4String& reaction, G4ParticleDefinition* pd)
 77   : G4VCrossSectionDataSet("JENDLHE" + reaction + "CrossSection")
 78 {
 79   reactionName = reaction;
 80   BuildPhysicsTable(*pd);
 81 }
 82 
 83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHEData() = default;
 84 
 85 void G4ParticleHPJENDLHEData::BuildPhysicsTable(const G4ParticleDefinition& aP)
 86 {
 87   particleName = aP.GetParticleName();
 88 
 89   const G4String& baseName = G4FindDataDir("G4NEUTRONHPDATA");
 90   const G4String& dirName = baseName + "/JENDL_HE/" + particleName + "/" + reactionName;
 91   const G4String& aFSType = "/CrossSection/";
 92   G4ParticleHPNames theNames;
 93 
 94   G4String filename;
 95 
 96   // Create JENDL_HE data
 97   // Create map element or isotope
 98 
 99   std::size_t numberOfElements = G4Element::GetNumberOfElements();
100 
101   // make a PhysicsVector for each element
102 
103   auto theElementTable = G4Element::GetElementTable();
104   vElement.clear();
105   vElement.resize(numberOfElements);
106   for (std::size_t i = 0; i < numberOfElements; ++i) {
107     G4Element* theElement = (*theElementTable)[i];
108     vElement[i] = false;
109 
110     // isotope
111     auto nIso = (G4int)(*theElementTable)[i]->GetNumberOfIsotopes();
112     auto Z = (G4int)(*theElementTable)[i]->GetZ();
113     for (G4int i1 = 0; i1 < nIso; ++i1) {
114       G4int A = theElement->GetIsotope(i1)->GetN();
115 
116       if (isThisNewIsotope(Z, A)) {
117   std::stringstream ss;
118   ss << dirName << aFSType << Z << "_" << A << "_" << theNames.GetName(Z - 1);
119   filename = ss.str();
120   std::fstream file;
121   file.open(filename, std::fstream::in);
122   G4int dummy;
123   file >> dummy;
124   if (file.good()) {
125     vElement[i] = true;
126 
127     // read the file
128     G4PhysicsVector* aPhysVec = readAFile(&file);
129     registAPhysicsVector(Z, A, aPhysVec);
130   }
131   file.close();
132       }
133     }
134   }
135 }
136 
137 void G4ParticleHPJENDLHEData::DumpPhysicsTable(const G4ParticleDefinition&)
138 {}
139 
140 G4double G4ParticleHPJENDLHEData::GetCrossSection(const G4DynamicParticle* aP,
141                                                   const G4Element* anE, G4double)
142 {
143   // Primary energy >20MeV
144   // Thus not taking into account of Doppler broadening
145   // also not taking into account of Target thermal motions
146 
147   G4double result = 0;
148 
149   G4double ek = aP->GetKineticEnergy();
150 
151   auto nIso = (G4int)anE->GetNumberOfIsotopes();
152   auto Z = (G4int)anE->GetZ();
153   for (G4int i1 = 0; i1 < nIso; ++i1) {
154     G4int A = anE->GetIsotope(i1)->GetN();
155     G4double frac = anE->GetRelativeAbundanceVector()[i1];
156     // This case does NOT request "*perCent".
157     result += frac * getXSfromThisIsotope(Z, A, ek);
158   }
159   return result;
160 }
161 
162 G4PhysicsVector* G4ParticleHPJENDLHEData::readAFile(std::fstream* file)
163 {
164   G4int dummy;
165   G4int len;
166   *file >> dummy;
167   *file >> len;
168 
169   std::vector<G4double> v_e;
170   std::vector<G4double> v_xs;
171 
172   for (G4int i = 0; i < len; ++i) {
173     G4double e;
174     G4double xs;
175 
176     *file >> e;
177     *file >> xs;
178     // data are written in eV and barn.
179     v_e.push_back(e * eV);
180     v_xs.push_back(xs * barn);
181   }
182 
183   auto aPhysVec = new G4PhysicsFreeVector(static_cast<std::size_t>(len), v_e.front(), v_e.back());
184 
185   for (G4int i = 0; i < len; ++i) {
186     aPhysVec->PutValues(static_cast<std::size_t>(i), v_e[i], v_xs[i]);
187   }
188 
189   return aPhysVec;
190 }
191 
192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4int z, G4int a)
193 {
194   if (mIsotope.find(z) == mIsotope.end()) return false;
195   if (mIsotope.find(z)->second->find(a) == mIsotope.find(z)->second->end()) return false;
196   return true;
197 }
198 
199 void G4ParticleHPJENDLHEData::registAPhysicsVector(G4int Z, G4int A, G4PhysicsVector* aPhysVec)
200 {
201   std::pair<G4int, G4PhysicsVector*> aPair = std::pair<G4int, G4PhysicsVector*>(A, aPhysVec);
202   auto itm = mIsotope.find(Z);
203   if (itm != mIsotope.cend()) {
204     itm->second->insert(aPair);
205   }
206   else {
207     auto aMap = new std::map<G4int, G4PhysicsVector*>;
208     aMap->insert(aPair);
209     mIsotope.insert(std::pair<G4int, std::map<G4int, G4PhysicsVector*>*>(Z, aMap));
210   }
211 }
212 
213 G4double G4ParticleHPJENDLHEData::getXSfromThisIsotope(G4int Z, G4int A, G4double ek)
214 {
215   G4double aXSection = 0.0;
216 
217   G4PhysicsVector* aPhysVec;
218   if (mIsotope.find(Z)->second->find(A) != mIsotope.find(Z)->second->end()) {
219     aPhysVec = mIsotope.find(Z)->second->find(A)->second;
220     aXSection = aPhysVec->Value(ek);
221   }
222   else {
223     // Select closest one in the same Z
224     G4int delta0 = 99;  // no mean for 99
225     for (auto it = mIsotope.find(Z)->second->cbegin(); it != mIsotope.find(Z)->second->cend(); ++it)
226     {
227       G4int delta = std::abs(A - it->first);
228       if (delta < delta0) delta0 = delta;
229     }
230 
231     // Randomize of selection larger or smaller than A
232     if (G4UniformRand() < 0.5) delta0 *= -1;
233     G4int A1 = A + delta0;
234     if (mIsotope.find(Z)->second->find(A1) != mIsotope.find(Z)->second->cend()) {
235       aPhysVec = mIsotope.find(Z)->second->find(A1)->second;
236     }
237     else {
238       A1 = A - delta0;
239       aPhysVec = mIsotope.find(Z)->second->find(A1)->second;
240     }
241 
242     aXSection = aPhysVec->Value(ek);
243     // X^(2/3) factor
244     aXSection *= G4Pow::GetInstance()->A23(1.0 * A / A1);
245   }
246 
247   return aXSection;
248 }
249