00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include "G4ShellData.hh"
00039 #include "G4DataVector.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include <fstream>
00042 #include <sstream>
00043 #include <numeric>
00044 #include <algorithm>
00045 #include <valarray>
00046 #include <functional>
00047 #include "Randomize.hh"
00048
00049
00050
00051 G4ShellData::G4ShellData(G4int minZ, G4int maxZ, G4bool isOccupancy)
00052 : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
00053 { }
00054
00055
00056 G4ShellData::~G4ShellData()
00057 {
00058 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos;
00059 for (pos = idMap.begin(); pos != idMap.end(); ++pos)
00060 {
00061 std::vector<G4double>* dataSet = (*pos).second;
00062 delete dataSet;
00063 }
00064
00065 std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos2;
00066 for (pos2 = bindingMap.begin(); pos2 != bindingMap.end(); ++pos2)
00067 {
00068 G4DataVector* dataSet = (*pos2).second;
00069 delete dataSet;
00070 }
00071
00072 if (occupancyData)
00073 {
00074 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos3;
00075 for (pos3 = occupancyPdfMap.begin(); pos3 != occupancyPdfMap.end(); ++pos3)
00076 {
00077 std::vector<G4double>* dataSet = (*pos3).second;
00078 delete dataSet;
00079 }
00080 }
00081 }
00082
00083
00084 size_t G4ShellData::NumberOfShells(G4int Z) const
00085 {
00086 G4int z = Z - 1;
00087 G4int n = 0;
00088
00089 if (Z>= zMin && Z <= zMax)
00090 {
00091 n = nShells[z];
00092 }
00093 return n;
00094 }
00095
00096
00097 const std::vector<G4double>& G4ShellData::ShellIdVector(G4int Z) const
00098 {
00099 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
00100 if (Z < zMin || Z > zMax) {
00101
00102 G4Exception("G4ShellData::ShellIdVector","de0001",FatalErrorInArgument, "Z outside boundaries");
00103
00104
00105 } pos = idMap.find(Z);
00106 std::vector<G4double>* dataSet = (*pos).second;
00107 return *dataSet;
00108 }
00109
00110
00111 const std::vector<G4double>& G4ShellData::ShellVector(G4int Z) const
00112 {
00113 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
00114 if (Z < zMin || Z > zMax) G4Exception("G4ShellData::ShellVector()","de0001",JustWarning,"Z outside boundaries");
00115 pos = occupancyPdfMap.find(Z);
00116 std::vector<G4double>* dataSet = (*pos).second;
00117 return *dataSet;
00118 }
00119
00120
00121 G4int G4ShellData::ShellId(G4int Z, G4int shellIndex) const
00122 {
00123 G4int n = -1;
00124
00125 if (Z >= zMin && Z <= zMax)
00126 {
00127 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
00128 pos = idMap.find(Z);
00129 if (pos!= idMap.end())
00130 {
00131 std::vector<G4double> dataSet = *((*pos).second);
00132 G4int nData = dataSet.size();
00133 if (shellIndex >= 0 && shellIndex < nData)
00134 {
00135 n = (G4int) dataSet[shellIndex];
00136 }
00137 }
00138 }
00139 return n;
00140 }
00141
00142
00143 G4double G4ShellData::ShellOccupancyProbability(G4int Z, G4int shellIndex) const
00144 {
00145 G4double prob = -1.;
00146
00147 if (Z >= zMin && Z <= zMax)
00148 {
00149 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
00150 pos = idMap.find(Z);
00151 if (pos!= idMap.end())
00152 {
00153 std::vector<G4double> dataSet = *((*pos).second);
00154 G4int nData = dataSet.size();
00155 if (shellIndex >= 0 && shellIndex < nData)
00156 {
00157 prob = dataSet[shellIndex];
00158 }
00159 }
00160 }
00161 return prob;
00162 }
00163
00164
00165
00166 G4double G4ShellData::BindingEnergy(G4int Z, G4int shellIndex) const
00167 {
00168 G4double value = 0.;
00169
00170 if (Z >= zMin && Z <= zMax)
00171 {
00172 std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
00173 pos = bindingMap.find(Z);
00174 if (pos!= bindingMap.end())
00175 {
00176 G4DataVector dataSet = *((*pos).second);
00177 G4int nData = dataSet.size();
00178 if (shellIndex >= 0 && shellIndex < nData)
00179 {
00180 value = dataSet[shellIndex];
00181 }
00182 }
00183 }
00184 return value;
00185 }
00186
00187 void G4ShellData::PrintData() const
00188 {
00189 for (G4int Z = zMin; Z <= zMax; Z++)
00190 {
00191 G4cout << "---- Shell data for Z = "
00192 << Z
00193 << " ---- "
00194 << G4endl;
00195 G4int nSh = nShells[Z-1];
00196 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posId;
00197 posId = idMap.find(Z);
00198 std::vector<G4double>* ids = (*posId).second;
00199 std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posE;
00200 posE = bindingMap.find(Z);
00201 G4DataVector* energies = (*posE).second;
00202 for (G4int i=0; i<nSh; i++)
00203 {
00204 G4int id = (G4int) (*ids)[i];
00205 G4double e = (*energies)[i] / keV;
00206 G4cout << i << ") ";
00207
00208 if (occupancyData)
00209 {
00210 G4cout << " Occupancy: ";
00211 }
00212 else
00213 {
00214 G4cout << " Shell id: ";
00215 }
00216 G4cout << id << " - Binding energy = "
00217 << e << " keV ";
00218 if (occupancyData)
00219 {
00220 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posOcc;
00221 posOcc = occupancyPdfMap.find(Z);
00222 std::vector<G4double> probs = *((*posOcc).second);
00223 G4double prob = probs[i];
00224 G4cout << "- Probability = " << prob;
00225 }
00226 G4cout << G4endl;
00227 }
00228 G4cout << "-------------------------------------------------"
00229 << G4endl;
00230 }
00231 }
00232
00233
00234 void G4ShellData::LoadData(const G4String& fileName)
00235 {
00236
00237
00238 std::ostringstream ost;
00239
00240 ost << fileName << ".dat";
00241
00242 G4String name(ost.str());
00243
00244 char* path = getenv("G4LEDATA");
00245 if (!path)
00246 {
00247 G4String excep("G4ShellData::LoadData()");
00248 G4Exception(excep,"em0006",FatalException,"Please set G4LEDATA");
00249 return;
00250 }
00251
00252 G4String pathString(path);
00253 G4String dirFile = pathString + name;
00254 std::ifstream file(dirFile);
00255 std::filebuf* lsdp = file.rdbuf();
00256
00257 if (! (lsdp->is_open()) )
00258 {
00259
00260 G4String excep = "G4ShellData::LoadData()";
00261 G4String msg = "data file: " + dirFile + " not found";
00262 G4Exception(excep, "em0003",FatalException, msg );
00263 return;
00264 }
00265
00266 G4double a = 0;
00267 G4int k = 1;
00268 G4int sLocal = 0;
00269
00270 G4int Z = 1;
00271 G4DataVector* energies = new G4DataVector;
00272 std::vector<G4double>* ids = new std::vector<G4double>;
00273
00274 do {
00275 file >> a;
00276 G4int nColumns = 2;
00277 if (a == -1)
00278 {
00279 if (sLocal == 0)
00280 {
00281
00282 idMap[Z] = ids;
00283 bindingMap[Z] = energies;
00284 G4int n = ids->size();
00285 nShells.push_back(n);
00286
00287
00288 ids = new std::vector<G4double>;
00289 energies = new G4DataVector;
00290 Z++;
00291 }
00292 sLocal++;
00293 if (sLocal == nColumns)
00294 {
00295 sLocal = 0;
00296 }
00297 }
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 else
00308 {
00309
00310 if(k%nColumns != 0)
00311 {
00312 ids->push_back(a);
00313 k++;
00314 }
00315 else if (k%nColumns == 0)
00316 {
00317
00318 G4double e = a * MeV;
00319 energies->push_back(e);
00320 k = 1;
00321 }
00322 }
00323 } while (a != -2);
00324 file.close();
00325 delete energies;
00326 delete ids;
00327
00328
00329
00330
00331 if (occupancyData)
00332 {
00333
00334
00335 for (G4int ZLocal=zMin; ZLocal <= zMax; ZLocal++)
00336 {
00337 std::vector<G4double> occupancy = ShellIdVector(ZLocal);
00338
00339 std::vector<G4double>* prob = new std::vector<G4double>;
00340 G4double scale = 1. / G4double(ZLocal);
00341
00342 prob->push_back(occupancy[0] * scale);
00343 for (size_t i=1; i<occupancy.size(); i++)
00344 {
00345 prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
00346 }
00347 occupancyPdfMap[ZLocal] = prob;
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358 }
00359 }
00360 }
00361
00362
00363 G4int G4ShellData::SelectRandomShell(G4int Z) const
00364 {
00365 if (Z < zMin || Z > zMax) {
00366
00367 G4Exception("G4ShellData::SelectrandomShell","de0001",FatalErrorInArgument, "Z outside boundaries");
00368
00369 }
00370 G4int shellIndex = 0;
00371 std::vector<G4double> prob = ShellVector(Z);
00372 G4double random = G4UniformRand();
00373
00374
00375
00376
00377
00378
00379 G4int nShellsLocal = NumberOfShells(Z);
00380 G4int upperBound = nShellsLocal;
00381
00382 while (shellIndex <= upperBound)
00383 {
00384 G4int midShell = (shellIndex + upperBound) / 2;
00385 if ( random < prob[midShell] )
00386 upperBound = midShell - 1;
00387 else
00388 shellIndex = midShell + 1;
00389 }
00390 if (shellIndex >= nShellsLocal) shellIndex = nShellsLocal - 1;
00391
00392 return shellIndex;
00393 }