Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4ShellData.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/electromagnetic/lowenergy/src/G4ShellData.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4ShellData.cc (Version 9.1.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 //                                                 26 //
                                                   >>  27 // $Id: G4ShellData.cc,v 1.8 2006/06/29 19:41:21 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
 27 //                                                 29 //
 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@     30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
 29 //                                                 31 //
 30 // History:                                        32 // History:
 31 // -----------                                     33 // -----------
 32 // 31 Jul 2001   MGP        Created                34 // 31 Jul 2001   MGP        Created
 33 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war << 
 34 //                                                 35 //
 35 // -------------------------------------------     36 // -------------------------------------------------------------------
 36                                                    37 
 37 #include "G4ShellData.hh"                          38 #include "G4ShellData.hh"
 38 #include "G4DataVector.hh"                         39 #include "G4DataVector.hh"
 39 #include "G4SystemOfUnits.hh"                  << 
 40 #include <fstream>                                 40 #include <fstream>
 41 #include <sstream>                                 41 #include <sstream>
 42 #include <numeric>                             << 
 43 #include <algorithm>                           << 
 44 #include <valarray>                            << 
 45 #include <functional>                          << 
 46 #include "Randomize.hh"                        << 
 47                                                << 
 48 //....oooOO0OOooo........oooOO0OOooo........oo << 
 49                                                << 
 50 G4ShellData::G4ShellData(G4int minZ, G4int max << 
 51   : zMin(minZ), zMax(maxZ), occupancyData(isOc << 
 52 {}                                             << 
 53                                                    42 
 54 //....oooOO0OOooo........oooOO0OOooo........oo <<  43 // Constructor
 55                                                    44 
                                                   >>  45 G4ShellData::G4ShellData(G4int minZ, G4int maxZ)
                                                   >>  46   : zMin(minZ), zMax(maxZ)
                                                   >>  47 { }
                                                   >>  48 
                                                   >>  49 // Destructor
 56 G4ShellData::~G4ShellData()                        50 G4ShellData::~G4ShellData()
 57 {                                                  51 {
 58   for (auto& pos : idMap)                      <<  52   std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos;
                                                   >>  53 
                                                   >>  54   for (pos = idMap.begin(); pos != idMap.end(); ++pos)
 59     {                                              55     {
 60       std::vector<G4double>* dataSet = pos.sec <<  56       G4DataVector* dataSet = (*pos).second;
 61       delete dataSet;                              57       delete dataSet;
 62     }                                              58     }
 63   for (auto& pos2 : bindingMap)                <<  59   for (pos = bindingMap.begin(); pos != bindingMap.end(); ++pos)
 64     {                                              60     {
 65       G4DataVector* dataSet = pos2.second;     <<  61       G4DataVector* dataSet = (*pos).second;
 66       delete dataSet;                              62       delete dataSet;
 67     }                                              63     }
 68                                                << 
 69   if (occupancyData)                           << 
 70     {                                          << 
 71       for (auto& pos3 : occupancyPdfMap)       << 
 72   {                                            << 
 73     std::vector<G4double>* dataSet = pos3.seco << 
 74     delete dataSet;                            << 
 75   }                                            << 
 76     }                                          << 
 77 }                                                  64 }
 78                                                    65 
 79 //....oooOO0OOooo........oooOO0OOooo........oo << 
 80                                                    66 
 81 std::size_t G4ShellData::NumberOfShells(G4int  <<  67 size_t G4ShellData::NumberOfShells(G4int Z) const
 82 {                                                  68 {
 83   G4int z = Z - 1;                                 69   G4int z = Z - 1;
 84   G4int n = 0;                                     70   G4int n = 0;
 85                                                    71 
 86   if (Z>= zMin && Z <= zMax)                       72   if (Z>= zMin && Z <= zMax)
 87     {                                              73     {
 88       n = nShells[z];                              74       n = nShells[z];
 89     }                                              75     }
 90   return n;                                        76   return n;
 91 }                                                  77 }
 92                                                    78 
 93 //....oooOO0OOooo........oooOO0OOooo........oo << 
 94                                                    79 
 95 const std::vector<G4double>& G4ShellData::Shel <<  80 const G4DataVector& G4ShellData::ShellIdVector(G4int Z) const
 96 {                                                  81 {
 97   if (Z < zMin || Z > zMax) {                  <<  82   std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
 98     G4Exception("G4ShellData::ShellIdVector"," <<  83   if (Z < zMin || Z > zMax)
 99   }                                            <<  84     G4Exception("G4ShellData::ShellIdVector - Z outside boundaries");
100   auto pos = idMap.find(Z);                    <<  85   pos = idMap.find(Z);
101   std::vector<G4double>* dataSet = (*pos).seco <<  86   G4DataVector* dataSet = (*pos).second;
102   return *dataSet;                                 87   return *dataSet;
103 }                                                  88 }
104                                                    89 
105 //....oooOO0OOooo........oooOO0OOooo........oo << 
106                                                << 
107 const std::vector<G4double>& G4ShellData::Shel << 
108 {                                              << 
109   if (Z < zMin || Z > zMax)                    << 
110     G4Exception("G4ShellData::ShellVector()"," << 
111   auto pos = occupancyPdfMap.find(Z);          << 
112   std::vector<G4double>* dataSet = (*pos).seco << 
113   return *dataSet;                             << 
114 }                                              << 
115                                                << 
116 //....oooOO0OOooo........oooOO0OOooo........oo << 
117                                                << 
118 G4int G4ShellData::ShellId(G4int Z, G4int shel     90 G4int G4ShellData::ShellId(G4int Z, G4int shellIndex) const
119 {                                                  91 {
120   G4int n = -1;                                    92   G4int n = -1;
121                                                    93 
122   if (Z >= zMin && Z <= zMax)                      94   if (Z >= zMin && Z <= zMax)
123     {                                              95     {
124       auto pos = idMap.find(Z);                <<  96       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
                                                   >>  97       pos = idMap.find(Z);
125       if (pos!= idMap.end())                       98       if (pos!= idMap.end())
126   {                                                99   {
127     std::vector<G4double> dataSet = *((*pos).s << 100     G4DataVector dataSet = *((*pos).second);
128     G4int nData = (G4int)dataSet.size();       << 101     G4int nData = dataSet.size();
129     if (shellIndex >= 0 && shellIndex < nData)    102     if (shellIndex >= 0 && shellIndex < nData)
130       {                                           103       {
131         n = (G4int) dataSet[shellIndex];          104         n = (G4int) dataSet[shellIndex];
132       }                                           105       }
133   }                                               106   }
134     }                                             107     }
135   return n;                                       108   return n;
136 }                                                 109 }
137                                                   110 
138 //....oooOO0OOooo........oooOO0OOooo........oo << 
139                                                << 
140 G4double G4ShellData::ShellOccupancyProbabilit << 
141 {                                              << 
142   G4double prob = -1.;                         << 
143   if (Z >= zMin && Z <= zMax)                  << 
144     {                                          << 
145       auto pos = idMap.find(Z);                << 
146       if (pos!= idMap.end())                   << 
147   {                                            << 
148     std::vector<G4double> dataSet = *((*pos).s << 
149     G4int nData = (G4int)dataSet.size();       << 
150     if (shellIndex >= 0 && shellIndex < nData) << 
151       {                                        << 
152         prob = dataSet[shellIndex];            << 
153       }                                        << 
154   }                                            << 
155     }                                          << 
156   return prob;                                 << 
157 }                                              << 
158                                                << 
159 //....oooOO0OOooo........oooOO0OOooo........oo << 
160                                                   111 
161 G4double G4ShellData::BindingEnergy(G4int Z, G    112 G4double G4ShellData::BindingEnergy(G4int Z, G4int shellIndex)  const
162 {                                                 113 {
163   G4double value = 0.;                            114   G4double value = 0.;
164                                                   115 
165   if (Z >= zMin && Z <= zMax)                     116   if (Z >= zMin && Z <= zMax)
166     {                                             117     {
167       auto pos = bindingMap.find(Z);           << 118       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
                                                   >> 119       pos = bindingMap.find(Z);
168       if (pos!= bindingMap.end())                 120       if (pos!= bindingMap.end())
169   {                                               121   {
170     G4DataVector dataSet = *((*pos).second);      122     G4DataVector dataSet = *((*pos).second);
171     G4int nData = (G4int)dataSet.size();       << 123     G4int nData = dataSet.size();
172     if (shellIndex >= 0 && shellIndex < nData)    124     if (shellIndex >= 0 && shellIndex < nData)
173       {                                           125       {
174         value = dataSet[shellIndex];              126         value = dataSet[shellIndex];
175       }                                           127       }
176   }                                               128   }
177     }                                             129     }
178   return value;                                   130   return value;
179 }                                                 131 }
180                                                   132 
181 //....oooOO0OOooo........oooOO0OOooo........oo << 
182                                                << 
183 void G4ShellData::PrintData() const               133 void G4ShellData::PrintData() const
184 {                                                 134 {
185   for (G4int Z = zMin; Z <= zMax; Z++)            135   for (G4int Z = zMin; Z <= zMax; Z++)
186     {                                             136     {
187       G4cout << "---- Shell data for Z = "        137       G4cout << "---- Shell data for Z = "
188        << Z                                       138        << Z
189        << " ---- "                                139        << " ---- "
190        << G4endl;                                 140        << G4endl;
191       G4int nSh = nShells[Z-1];                   141       G4int nSh = nShells[Z-1];
192       auto posId = idMap.find(Z);              << 142       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posId;
193       std::vector<G4double>* ids = (*posId).se << 143       posId = idMap.find(Z);
194       auto posE = bindingMap.find(Z);          << 144       G4DataVector* ids = (*posId).second;
                                                   >> 145       std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posE;
                                                   >> 146       posE = bindingMap.find(Z);
195       G4DataVector* energies = (*posE).second;    147       G4DataVector* energies = (*posE).second;
196       for (G4int i=0; i<nSh; ++i)              << 148       for (G4int i=0; i<nSh; i++)
197   {                                               149   {
198     G4int id = (G4int) (*ids)[i];                 150     G4int id = (G4int) (*ids)[i];
199     G4double e = (*energies)[i] / keV;         << 151     G4double e = (*energies)[i] / MeV;
200     G4cout << i << ") ";                       << 152     G4cout << i <<") Shell id: " << id 
201                                                << 153      << " - Binding energy = "
202     if (occupancyData)                         << 154      << e << " MeV " << G4endl;
203       {                                        << 
204         G4cout << " Occupancy: ";              << 
205       }                                        << 
206     else                                       << 
207       {                                        << 
208         G4cout << " Shell id: ";               << 
209       }                                        << 
210     G4cout << id << " - Binding energy = "     << 
211      << e << " keV ";                          << 
212       if (occupancyData)                       << 
213         {                                      << 
214     auto posOcc = occupancyPdfMap.find(Z);     << 
215                 std::vector<G4double> probs =  << 
216                 G4double prob = probs[i];      << 
217     G4cout << "- Probability = " << prob;      << 
218         }                                      << 
219       G4cout << G4endl;                        << 
220   }                                               155   }
221       G4cout << "-----------------------------    156       G4cout << "-------------------------------------------------" 
222        << G4endl;                                 157        << G4endl;
223     }                                             158     }
224 }                                                 159 }
225                                                   160 
226 //....oooOO0OOooo........oooOO0OOooo........oo << 
227                                                   161 
228 void G4ShellData::LoadData(const G4String& fil    162 void G4ShellData::LoadData(const G4String& fileName)
229 {                                                 163 { 
230   // Build the complete string identifying the    164   // Build the complete string identifying the file with the data set
                                                   >> 165   
231   std::ostringstream ost;                         166   std::ostringstream ost;
                                                   >> 167   
232   ost << fileName << ".dat";                      168   ost << fileName << ".dat";
                                                   >> 169   
233   G4String name(ost.str());                       170   G4String name(ost.str());
234                                                   171   
235   const char* path = G4FindDataDir("G4LEDATA") << 172   char* path = getenv("G4LEDATA");
236   if (!path)                                      173   if (!path)
237     {                                             174     { 
238       G4String excep("G4ShellData::LoadData()" << 175       G4String excep("G4EMDataSet - G4LEDATA environment variable not set");
239       G4Exception(excep,"em0006",FatalExceptio << 176       G4Exception(excep);
240       return;                                  << 
241     }                                             177     }
242                                                   178   
243   G4String pathString(path);                      179   G4String pathString(path);
244   G4String dirFile = pathString + name;           180   G4String dirFile = pathString + name;
245   std::ifstream file(dirFile);                    181   std::ifstream file(dirFile);
246   std::filebuf* lsdp = file.rdbuf();              182   std::filebuf* lsdp = file.rdbuf();
247                                                   183 
248   if (! (lsdp->is_open()) )                       184   if (! (lsdp->is_open()) )
249     {                                             185     {
250                                                << 186       G4String s1("G4ShellData - data file: ");
251       G4String excep = "G4ShellData::LoadData( << 187       G4String s2(" not found");
252       G4String msg = "data file: " + dirFile + << 188       G4String excep = s1 + dirFile + s2;
253       G4Exception(excep, "em0003",FatalExcepti << 189       G4Exception(excep);
254       return;                                  << 
255     }                                             190     }
256                                                   191 
257   G4double a = 0;                                 192   G4double a = 0;
258   G4int k = 1;                                    193   G4int k = 1;
259   G4int sLocal = 0;                            << 194   G4int s = 0;
260                                                   195   
261   G4int Z = 1;                                    196   G4int Z = 1;
262   G4DataVector* energies = new G4DataVector;      197   G4DataVector* energies = new G4DataVector;
263   std::vector<G4double>* ids = new std::vector << 198   G4DataVector* ids = new G4DataVector;
264                                                   199 
265   do {                                            200   do {
266     file >> a;                                    201     file >> a;
267     G4int nColumns = 2;                           202     G4int nColumns = 2;
268     if (a == -1)                                  203     if (a == -1)
269       {                                           204       {
270   if (sLocal == 0)                             << 205   if (s == 0)
271     {                                             206     {
272       // End of a shell data set                  207       // End of a shell data set
273       idMap[Z] = ids;                             208       idMap[Z] = ids;
274             bindingMap[Z] = energies;             209             bindingMap[Z] = energies;
275             G4int n = (G4int)ids->size();      << 210             G4int n = ids->size();
276       nShells.push_back(n);                       211       nShells.push_back(n);
277       // Start of new shell data set              212       // Start of new shell data set
278                                                << 213       ids = new G4DataVector;
279       ids = new std::vector<G4double>;         << 
280             energies = new G4DataVector;          214             energies = new G4DataVector;
281             Z++;                                  215             Z++;      
282     }                                             216     }      
283   sLocal++;                                    << 217   s++;
284   if (sLocal == nColumns)                      << 218   if (s == nColumns)
285   {                                               219   {
286     sLocal = 0;                                << 220     s = 0;
287   }                                               221   }
288       }                                           222       }
289                                                << 223     else if (a == -2)
290     // moved out of the do-while since might g << 224       {
291     //    else if (a == -2)                    << 225   // End of file; delete the empty vectors created when encountering the last -1 -1 row
292     //      {                                  << 226   delete energies;
293     // End of file; delete the empty vectors c << 227   delete ids;
294     //  delete energies;                       << 228   //nComponents = components.size();
295     //  delete ids;                            << 229       }
296     //nComponents = components.size();         << 
297     //      }                                  << 
298     else                                          230     else
299       {                                           231       {
300   // 1st column is shell id                       232   // 1st column is shell id
301   if(k%nColumns != 0)                             233   if(k%nColumns != 0)
302     {                                             234     {     
303       ids->push_back(a);                          235       ids->push_back(a);
304       k++;                                        236       k++;
305     }                                             237     }
306   else if (k%nColumns == 0)                       238   else if (k%nColumns == 0)
307     {                                             239     {
308       // 2nd column is binding energy             240       // 2nd column is binding energy
309       G4double e = a * MeV;                       241       G4double e = a * MeV;
310       energies->push_back(e);                     242       energies->push_back(e);
311       k = 1;                                      243       k = 1;
312     }                                             244     }
313       }                                           245       }
314   } while (a != -2); // end of file               246   } while (a != -2); // end of file
315   file.close();                                   247   file.close();    
316   delete energies;                             << 
317   delete ids;                                  << 
318                                                << 
319   // For Doppler broadening: the data set cont << 
320   // Build additional map with probability for << 
321                                                << 
322   if (occupancyData)                           << 
323     {                                          << 
324       // Build cumulative from raw shell occup << 
325       for (G4int ZLocal=zMin; ZLocal <= zMax;  << 
326   {                                            << 
327     std::vector<G4double> occupancy = ShellIdV << 
328                                                << 
329     std::vector<G4double>* prob = new std::vec << 
330     G4double scale = 1. / G4double(ZLocal);    << 
331                                                << 
332     prob->push_back(occupancy[0] * scale);     << 
333     for (std::size_t i=1; i<occupancy.size();  << 
334       {                                        << 
335         prob->push_back(occupancy[i]*scale + ( << 
336       }                                        << 
337     occupancyPdfMap[ZLocal] = prob;            << 
338   }                                            << 
339     }                                          << 
340 }                                                 248 }
341                                                   249 
342 //....oooOO0OOooo........oooOO0OOooo........oo << 
343                                                << 
344 G4int G4ShellData::SelectRandomShell(G4int Z)  << 
345 {                                              << 
346   if (Z < zMin || Z > zMax)                    << 
347     G4Exception("G4ShellData::SelectrandomShel << 
348                                                << 
349   G4int shellIndex = 0;                        << 
350   std::vector<G4double> prob = ShellVector(Z); << 
351   G4double random = G4UniformRand();           << 
352                                                << 
353   // Binary search the shell with probability  << 
354   G4int nShellsLocal = (G4int)NumberOfShells(Z << 
355   G4int upperBound = nShellsLocal;             << 
356                                                << 
357   while (shellIndex <= upperBound)             << 
358     {                                          << 
359       G4int midShell = (shellIndex + upperBoun << 
360       if ( random < prob[midShell] )           << 
361   upperBound = midShell - 1;                   << 
362       else                                     << 
363   shellIndex = midShell + 1;                   << 
364     }                                          << 
365   if (shellIndex >= nShellsLocal) shellIndex = << 
366                                                << 
367   return shellIndex;                           << 
368 }                                              << 
369                                                   250