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 ]

Diff markup

Differences between /processes/hadronic/models/particle_hp/src/G4ParticleHPJENDLHEData.cc (Version 11.3.0) and /processes/hadronic/models/particle_hp/src/G4ParticleHPJENDLHEData.cc (Version 8.1)


  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 // Class Description                              
 27 // Cross-section data set for a high precision    
 28 // libraries) description of elastic scatterin    
 29 // Class Description - End                        
 30                                                   
 31 // 15-Nov-06 First Implementation is done by T    
 32 // P. Arce, June-2014 Conversion neutron_hp to    
 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(c    
 43 {                                                 
 44   G4bool result = true;                           
 45   G4double eKin = aP->GetKineticEnergy();         
 46   // if(eKin>20*MeV||aP->GetDefinition()!=G4Ne    
 47   if (eKin < 20 * MeV || 3 * GeV < eKin || aP-    
 48     result = false;                               
 49   }                                               
 50   // Element Check                                
 51   else if (!(vElement[anE->GetIndex()]))          
 52     result = false;                               
 53                                                   
 54   return result;                                  
 55 }                                                 
 56                                                   
 57 G4ParticleHPJENDLHEData::G4ParticleHPJENDLHEDa    
 58 {                                                 
 59   for (auto& itZ : mIsotope) {                    
 60     std::map<G4int, G4PhysicsVector*>* pointer    
 61     if (pointer_map != nullptr) {                 
 62       for (auto& itA : *pointer_map) {            
 63         G4PhysicsVector* pointerPhysicsVector     
 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::G4ParticleHPJENDLHEDa    
 77   : G4VCrossSectionDataSet("JENDLHE" + reactio    
 78 {                                                 
 79   reactionName = reaction;                        
 80   BuildPhysicsTable(*pd);                         
 81 }                                                 
 82                                                   
 83 G4ParticleHPJENDLHEData::~G4ParticleHPJENDLHED    
 84                                                   
 85 void G4ParticleHPJENDLHEData::BuildPhysicsTabl    
 86 {                                                 
 87   particleName = aP.GetParticleName();            
 88                                                   
 89   const G4String& baseName = G4FindDataDir("G4    
 90   const G4String& dirName = baseName + "/JENDL    
 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::Ge    
100                                                   
101   // make a PhysicsVector for each element        
102                                                   
103   auto theElementTable = G4Element::GetElement    
104   vElement.clear();                               
105   vElement.resize(numberOfElements);              
106   for (std::size_t i = 0; i < numberOfElements    
107     G4Element* theElement = (*theElementTable)    
108     vElement[i] = false;                          
109                                                   
110     // isotope                                    
111     auto nIso = (G4int)(*theElementTable)[i]->    
112     auto Z = (G4int)(*theElementTable)[i]->Get    
113     for (G4int i1 = 0; i1 < nIso; ++i1) {         
114       G4int A = theElement->GetIsotope(i1)->Ge    
115                                                   
116       if (isThisNewIsotope(Z, A)) {               
117   std::stringstream ss;                           
118   ss << dirName << aFSType << Z << "_" << A <<    
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(&fil    
129     registAPhysicsVector(Z, A, aPhysVec);         
130   }                                               
131   file.close();                                   
132       }                                           
133     }                                             
134   }                                               
135 }                                                 
136                                                   
137 void G4ParticleHPJENDLHEData::DumpPhysicsTable    
138 {}                                                
139                                                   
140 G4double G4ParticleHPJENDLHEData::GetCrossSect    
141                                                   
142 {                                                 
143   // Primary energy >20MeV                        
144   // Thus not taking into account of Doppler b    
145   // also not taking into account of Target th    
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->GetRelativeAbundanceV    
156     // This case does NOT request "*perCent".     
157     result += frac * getXSfromThisIsotope(Z, A    
158   }                                               
159   return result;                                  
160 }                                                 
161                                                   
162 G4PhysicsVector* G4ParticleHPJENDLHEData::read    
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(stat    
184                                                   
185   for (G4int i = 0; i < len; ++i) {               
186     aPhysVec->PutValues(static_cast<std::size_    
187   }                                               
188                                                   
189   return aPhysVec;                                
190 }                                                 
191                                                   
192 G4bool G4ParticleHPJENDLHEData::isThisInMap(G4    
193 {                                                 
194   if (mIsotope.find(z) == mIsotope.end()) retu    
195   if (mIsotope.find(z)->second->find(a) == mIs    
196   return true;                                    
197 }                                                 
198                                                   
199 void G4ParticleHPJENDLHEData::registAPhysicsVe    
200 {                                                 
201   std::pair<G4int, G4PhysicsVector*> aPair = s    
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, G4PhysicsV    
208     aMap->insert(aPair);                          
209     mIsotope.insert(std::pair<G4int, std::map<    
210   }                                               
211 }                                                 
212                                                   
213 G4double G4ParticleHPJENDLHEData::getXSfromThi    
214 {                                                 
215   G4double aXSection = 0.0;                       
216                                                   
217   G4PhysicsVector* aPhysVec;                      
218   if (mIsotope.find(Z)->second->find(A) != mIs    
219     aPhysVec = mIsotope.find(Z)->second->find(    
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->c    
226     {                                             
227       G4int delta = std::abs(A - it->first);      
228       if (delta < delta0) delta0 = delta;         
229     }                                             
230                                                   
231     // Randomize of selection larger or smalle    
232     if (G4UniformRand() < 0.5) delta0 *= -1;      
233     G4int A1 = A + delta0;                        
234     if (mIsotope.find(Z)->second->find(A1) !=     
235       aPhysVec = mIsotope.find(Z)->second->fin    
236     }                                             
237     else {                                        
238       A1 = A - delta0;                            
239       aPhysVec = mIsotope.find(Z)->second->fin    
240     }                                             
241                                                   
242     aXSection = aPhysVec->Value(ek);              
243     // X^(2/3) factor                             
244     aXSection *= G4Pow::GetInstance()->A23(1.0    
245   }                                               
246                                                   
247   return aXSection;                               
248 }                                                 
249