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