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