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 // History: 26 // History: 27 // ----------- 27 // ----------- 28 // 31 Jul 2001 MGP Created 28 // 31 Jul 2001 MGP Created 29 // 29 // 30 // 15 Jul 2009 Nicolas A. Karakatsanis 30 // 15 Jul 2009 Nicolas A. Karakatsanis 31 // 31 // 32 // - New Constructor 32 // - New Constructor was added when logarithmic data are loaded as well 33 // to enhance CPU 33 // to enhance CPU performance 34 // 34 // 35 // - Destructor was 35 // - Destructor was modified accordingly 36 // 36 // 37 // - LoadNonLogData 37 // - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW 38 // dataset. It is 38 // dataset. It is essentially performing the data loading operations as in the past. 39 // 39 // 40 // - LoadData method 40 // - LoadData method was revised in order to calculate the logarithmic values of the data 41 // It retrieves th 41 // It retrieves the data values from the G4EMLOW data files but, then, calculates the 42 // respective log 42 // respective log values and loads them to seperate data structures. 43 // 43 // 44 // - SetLogEnergiesD 44 // - SetLogEnergiesData method was cretaed to set logarithmic values to G4 data vectors. 45 // The EM data set 45 // The EM data sets, initialized this way, contain both non-log and log values. 46 // These initializ 46 // These initialized data sets can enhance the computing performance of data interpolation 47 // operations 47 // operations 48 // 48 // 49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war 49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings and cleanup logic 50 // 50 // 51 // ------------------------------------------- 51 // ------------------------------------------------------------------- 52 52 53 #include "G4EMDataSet.hh" 53 #include "G4EMDataSet.hh" 54 #include "G4VDataSetAlgorithm.hh" 54 #include "G4VDataSetAlgorithm.hh" 55 #include <fstream> 55 #include <fstream> 56 #include <sstream> 56 #include <sstream> 57 #include "G4Integrator.hh" 57 #include "G4Integrator.hh" 58 #include "Randomize.hh" 58 #include "Randomize.hh" 59 #include "G4LinInterpolation.hh" 59 #include "G4LinInterpolation.hh" 60 60 >> 61 61 G4EMDataSet::G4EMDataSet(G4int Z, 62 G4EMDataSet::G4EMDataSet(G4int Z, 62 G4VDataSetAlgorithm* algo, 63 G4VDataSetAlgorithm* algo, 63 G4double xUnit, 64 G4double xUnit, 64 G4double yUnit, 65 G4double yUnit, 65 G4bool random): << 66 G4bool random): 66 energies(nullptr), << 67 z(Z), 67 data(nullptr), << 68 energies(0), 68 log_energies(nullptr), << 69 data(0), 69 log_data(nullptr), << 70 log_energies(0), >> 71 log_data(0), 70 algorithm(algo), 72 algorithm(algo), 71 pdf(nullptr), << 72 unitEnergies(xUnit), 73 unitEnergies(xUnit), 73 unitData(yUnit), << 74 unitData(yUnit), 74 z(Z), << 75 pdf(0), 75 randomSet(random) 76 randomSet(random) 76 { 77 { 77 if (algorithm == nullptr) { << 78 if (algorithm == 0) { 78 G4Exception("G4EMDataSet::G4EMDataSet", 79 G4Exception("G4EMDataSet::G4EMDataSet", 79 "em1012",FatalException,"interpolation 80 "em1012",FatalException,"interpolation == 0"); 80 } else if (randomSet) { BuildPdf(); } 81 } else if (randomSet) { BuildPdf(); } 81 } 82 } 82 83 83 //....oooOO0OOooo........oooOO0OOooo........oo << 84 << 85 G4EMDataSet::G4EMDataSet(G4int argZ, 84 G4EMDataSet::G4EMDataSet(G4int argZ, 86 G4DataVector* dataX, 85 G4DataVector* dataX, 87 G4DataVector* dataY, 86 G4DataVector* dataY, 88 G4VDataSetAlgorithm* algo, 87 G4VDataSetAlgorithm* algo, 89 G4double xUnit, 88 G4double xUnit, 90 G4double yUnit, 89 G4double yUnit, 91 G4bool random): 90 G4bool random): >> 91 z(argZ), 92 energies(dataX), 92 energies(dataX), 93 data(dataY), 93 data(dataY), 94 log_energies(nullptr), << 94 log_energies(0), 95 log_data(nullptr), << 95 log_data(0), 96 algorithm(algo), 96 algorithm(algo), 97 pdf(nullptr), << 98 unitEnergies(xUnit), 97 unitEnergies(xUnit), 99 unitData(yUnit), << 98 unitData(yUnit), 100 z(argZ), << 99 pdf(0), 101 randomSet(random) 100 randomSet(random) 102 { 101 { 103 if (!algorithm || !energies || !data) { 102 if (!algorithm || !energies || !data) { 104 G4Exception("G4EMDataSet::G4EMDataSet", 103 G4Exception("G4EMDataSet::G4EMDataSet", 105 "em1012",FatalException,"interpolation 104 "em1012",FatalException,"interpolation == 0"); 106 } else { 105 } else { 107 if (energies->size() != data->size()) { 106 if (energies->size() != data->size()) { 108 G4Exception("G4EMDataSet::G4EMDataSet", 107 G4Exception("G4EMDataSet::G4EMDataSet", 109 "em1012",FatalException,"different size 108 "em1012",FatalException,"different size for energies and data"); 110 } else if (randomSet) { 109 } else if (randomSet) { 111 BuildPdf(); 110 BuildPdf(); 112 } 111 } 113 } 112 } 114 } 113 } 115 114 116 //....oooOO0OOooo........oooOO0OOooo........oo << 117 << 118 G4EMDataSet::G4EMDataSet(G4int argZ, 115 G4EMDataSet::G4EMDataSet(G4int argZ, 119 G4DataVector* dataX, 116 G4DataVector* dataX, 120 G4DataVector* dataY, 117 G4DataVector* dataY, 121 G4DataVector* dataLog 118 G4DataVector* dataLogX, 122 G4DataVector* dataLog 119 G4DataVector* dataLogY, 123 G4VDataSetAlgorithm* algo, 120 G4VDataSetAlgorithm* algo, 124 G4double xUnit, 121 G4double xUnit, 125 G4double yUnit, 122 G4double yUnit, 126 G4bool random): 123 G4bool random): >> 124 z(argZ), 127 energies(dataX), 125 energies(dataX), 128 data(dataY), 126 data(dataY), 129 log_energies(dataLogX), 127 log_energies(dataLogX), 130 log_data(dataLogY), 128 log_data(dataLogY), 131 algorithm(algo), 129 algorithm(algo), 132 pdf(nullptr), << 133 unitEnergies(xUnit), 130 unitEnergies(xUnit), 134 unitData(yUnit), 131 unitData(yUnit), 135 z(argZ), << 132 pdf(0), 136 randomSet(random) 133 randomSet(random) 137 { 134 { 138 if (!algorithm || !energies || !data || !log 135 if (!algorithm || !energies || !data || !log_energies || !log_data) { 139 G4Exception("G4EMDataSet::G4EMDataSet", 136 G4Exception("G4EMDataSet::G4EMDataSet", 140 "em1012",FatalException,"interpolation 137 "em1012",FatalException,"interpolation == 0"); 141 } else { 138 } else { 142 if ((energies->size() != data->size()) || 139 if ((energies->size() != data->size()) || 143 (energies->size() != log_energies->size()) | 140 (energies->size() != log_energies->size()) || 144 (energies->size() != log_data->size())) { 141 (energies->size() != log_data->size())) { 145 G4Exception("G4EMDataSet::G4EMDataSet", 142 G4Exception("G4EMDataSet::G4EMDataSet", 146 "em1012",FatalException,"different size 143 "em1012",FatalException,"different size for energies and data"); 147 } else if (randomSet) { 144 } else if (randomSet) { 148 BuildPdf(); 145 BuildPdf(); 149 } 146 } 150 } 147 } 151 } 148 } 152 149 153 //....oooOO0OOooo........oooOO0OOooo........oo << 154 << 155 G4EMDataSet::~G4EMDataSet() 150 G4EMDataSet::~G4EMDataSet() 156 { 151 { 157 delete algorithm; 152 delete algorithm; 158 delete energies; 153 delete energies; 159 delete data; 154 delete data; 160 delete pdf; 155 delete pdf; 161 delete log_energies; 156 delete log_energies; 162 delete log_data; 157 delete log_data; 163 } 158 } 164 159 165 //....oooOO0OOooo........oooOO0OOooo........oo << 166 << 167 G4double G4EMDataSet::FindValue(G4double energ 160 G4double G4EMDataSet::FindValue(G4double energy, G4int /* componentId */) const 168 { 161 { 169 if (energy <= (*energies)[0]) { 162 if (energy <= (*energies)[0]) { 170 return (*data)[0]; 163 return (*data)[0]; 171 } 164 } 172 std::size_t i = energies->size()-1; << 165 size_t i = energies->size()-1; 173 if (energy >= (*energies)[i]) { return (*dat 166 if (energy >= (*energies)[i]) { return (*data)[i]; } 174 167 175 //Nicolas A. Karakatsanis: Check if the loga 168 //Nicolas A. Karakatsanis: Check if the logarithmic data have been loaded to decide 176 // which Interpolati 169 // which Interpolation-Calculation method will be applied 177 if (log_energies != nullptr) << 170 if (log_energies != 0) 178 { 171 { 179 return algorithm->Calculate(energy, (G4in << 172 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data, *log_energies, *log_data); 180 } 173 } 181 174 182 return algorithm->Calculate(energy, (G4int)F << 175 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data); 183 } 176 } 184 177 185 //....oooOO0OOooo........oooOO0OOooo........oo << 186 << 187 void G4EMDataSet::PrintData(void) const 178 void G4EMDataSet::PrintData(void) const 188 { 179 { 189 std::size_t size = energies->size(); << 180 size_t size = energies->size(); 190 for (std::size_t i(0); i<size; ++i) << 181 for (size_t i(0); i<size; i++) 191 { 182 { 192 G4cout << "Point: " << ((*energies)[i]/u 183 G4cout << "Point: " << ((*energies)[i]/unitEnergies) 193 << " - Data value: " << ((*data)[i]/uni 184 << " - Data value: " << ((*data)[i]/unitData); 194 if (pdf != 0) G4cout << " - PDF : " << ( 185 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i]; 195 G4cout << G4endl; 186 G4cout << G4endl; 196 } 187 } 197 } 188 } 198 189 199 //....oooOO0OOooo........oooOO0OOooo........oo << 200 << 201 void G4EMDataSet::SetEnergiesData(G4DataVector 190 void G4EMDataSet::SetEnergiesData(G4DataVector* dataX, 202 G4DataVector* dataY, 191 G4DataVector* dataY, 203 G4int /* componentId */) 192 G4int /* componentId */) 204 { 193 { 205 if(!dataX || !dataY) { 194 if(!dataX || !dataY) { 206 G4Exception("G4EMDataSet::SetEnergiesData" 195 G4Exception("G4EMDataSet::SetEnergiesData", 207 "em1012",FatalException,"new interpola 196 "em1012",FatalException,"new interpolation == 0"); 208 } else { 197 } else { 209 if (dataX->size() != dataY->size()) { 198 if (dataX->size() != dataY->size()) { 210 G4Exception("G4EMDataSet::SetEnergiesDat 199 G4Exception("G4EMDataSet::SetEnergiesData", 211 "em1012",FatalException,"different size 200 "em1012",FatalException,"different size for energies and data"); 212 } else { 201 } else { 213 202 214 delete energies; 203 delete energies; 215 energies = dataX; 204 energies = dataX; 216 205 217 delete data; 206 delete data; 218 data = dataY; 207 data = dataY; 219 //G4cout << "Size of energies: " << ener 208 //G4cout << "Size of energies: " << energies->size() << G4endl 220 //<< "Size of data: " << data->size() << 209 //<< "Size of data: " << data->size() << G4endl; 221 } 210 } 222 } 211 } 223 } 212 } 224 213 225 //....oooOO0OOooo........oooOO0OOooo........oo << 226 << 227 void G4EMDataSet::SetLogEnergiesData(G4DataVec 214 void G4EMDataSet::SetLogEnergiesData(G4DataVector* dataX, 228 G4DataVec 215 G4DataVector* dataY, 229 G4DataVec 216 G4DataVector* data_logX, 230 G4DataVector* data_logY, 217 G4DataVector* data_logY, 231 G4int /* 218 G4int /* componentId */) 232 { 219 { 233 if(!dataX || !dataY || !data_logX || !data_l 220 if(!dataX || !dataY || !data_logX || !data_logY) { 234 G4Exception("G4EMDataSet::SetEnergiesData" 221 G4Exception("G4EMDataSet::SetEnergiesData", 235 "em1012",FatalException,"new interpola 222 "em1012",FatalException,"new interpolation == 0"); 236 } else { 223 } else { 237 if (dataX->size() != dataY->size() || 224 if (dataX->size() != dataY->size() || 238 dataX->size() != data_logX->size() || 225 dataX->size() != data_logX->size() || 239 dataX->size() != data_logY->size()) { 226 dataX->size() != data_logY->size()) { 240 G4Exception("G4EMDataSet::SetEnergiesDat 227 G4Exception("G4EMDataSet::SetEnergiesData", 241 "em1012",FatalException,"different size 228 "em1012",FatalException,"different size for energies and data"); 242 } else { 229 } else { 243 230 244 delete energies; 231 delete energies; 245 energies = dataX; 232 energies = dataX; 246 233 247 delete data; 234 delete data; 248 data = dataY; 235 data = dataY; 249 236 250 delete log_energies; 237 delete log_energies; 251 log_energies = data_logX; 238 log_energies = data_logX; 252 239 253 delete log_data; 240 delete log_data; 254 log_data = data_logY; 241 log_data = data_logY; 255 //G4cout << "Size of energies: " << ener 242 //G4cout << "Size of energies: " << energies->size() << G4endl 256 //<< "Size of data: " << data->size() << 243 //<< "Size of data: " << data->size() << G4endl; 257 } 244 } 258 } 245 } 259 } 246 } 260 247 261 //....oooOO0OOooo........oooOO0OOooo........oo << 262 << 263 G4bool G4EMDataSet::LoadData(const G4String& f 248 G4bool G4EMDataSet::LoadData(const G4String& fileName) 264 { 249 { 265 // The file is organized into four columns: 250 // The file is organized into four columns: 266 // 1st column contains the values of energy 251 // 1st column contains the values of energy 267 // 2nd column contains the corresponding data 252 // 2nd column contains the corresponding data value 268 // The file terminates with the pattern: -1 253 // The file terminates with the pattern: -1 -1 269 // -2 254 // -2 -2 270 255 271 G4String fullFileName(FullFileName(fileName) 256 G4String fullFileName(FullFileName(fileName)); 272 std::ifstream in(fullFileName); 257 std::ifstream in(fullFileName); 273 258 274 if (!in.is_open()) 259 if (!in.is_open()) 275 { 260 { 276 G4String message("data file \""); 261 G4String message("data file \""); 277 message += fullFileName; 262 message += fullFileName; 278 message += "\" not found"; 263 message += "\" not found"; 279 G4Exception("G4EMDataSet::LoadData", 264 G4Exception("G4EMDataSet::LoadData", 280 "em1012",FatalException,message); 265 "em1012",FatalException,message); 281 return false; 266 return false; 282 } 267 } 283 268 284 delete energies; 269 delete energies; 285 delete data; 270 delete data; 286 delete log_energies; 271 delete log_energies; 287 delete log_data; 272 delete log_data; 288 energies = new G4DataVector; 273 energies = new G4DataVector; 289 data = new G4DataVector; 274 data = new G4DataVector; 290 log_energies = new G4DataVector; 275 log_energies = new G4DataVector; 291 log_data = new G4DataVector; 276 log_data = new G4DataVector; 292 277 293 G4double a, b; 278 G4double a, b; >> 279 //G4int k = 0; >> 280 //G4int nColumns = 2; >> 281 294 do 282 do 295 { 283 { 296 in >> a >> b; 284 in >> a >> b; 297 285 298 if (a != -1 && a != -2) 286 if (a != -1 && a != -2) 299 { 287 { 300 if (a==0.) { a=1e-300; } 288 if (a==0.) { a=1e-300; } 301 if (b==0.) { b=1e-300; } 289 if (b==0.) { b=1e-300; } 302 a *= unitEnergies; 290 a *= unitEnergies; 303 b *= unitData; 291 b *= unitData; 304 energies->push_back(a); 292 energies->push_back(a); 305 log_energies->push_back(std::log10(a)); 293 log_energies->push_back(std::log10(a)); 306 data->push_back(b); 294 data->push_back(b); 307 log_data->push_back(std::log10(b)); 295 log_data->push_back(std::log10(b)); 308 } 296 } 309 } 297 } 310 while (a != -2); 298 while (a != -2); 311 299 312 if (randomSet) { BuildPdf(); } 300 if (randomSet) { BuildPdf(); } 313 301 314 return true; 302 return true; 315 } 303 } 316 304 317 //....oooOO0OOooo........oooOO0OOooo........oo << 318 305 319 G4bool G4EMDataSet::LoadNonLogData(const G4Str 306 G4bool G4EMDataSet::LoadNonLogData(const G4String& fileName) 320 { 307 { 321 // The file is organized into four columns: 308 // The file is organized into four columns: 322 // 1st column contains the values of energy 309 // 1st column contains the values of energy 323 // 2nd column contains the corresponding data 310 // 2nd column contains the corresponding data value 324 // The file terminates with the pattern: -1 311 // The file terminates with the pattern: -1 -1 325 // -2 312 // -2 -2 326 313 327 G4String fullFileName(FullFileName(fileName) 314 G4String fullFileName(FullFileName(fileName)); 328 std::ifstream in(fullFileName); 315 std::ifstream in(fullFileName); 329 if (!in.is_open()) 316 if (!in.is_open()) 330 { 317 { 331 G4String message("data file \""); 318 G4String message("data file \""); 332 message += fullFileName; 319 message += fullFileName; 333 message += "\" not found"; 320 message += "\" not found"; 334 G4Exception("G4EMDataSet::LoadNonLogData 321 G4Exception("G4EMDataSet::LoadNonLogData", 335 "em1012",FatalException,message); 322 "em1012",FatalException,message); 336 } 323 } 337 324 338 G4DataVector* argEnergies=new G4DataVector; 325 G4DataVector* argEnergies=new G4DataVector; 339 G4DataVector* argData=new G4DataVector; 326 G4DataVector* argData=new G4DataVector; 340 327 341 G4double a; 328 G4double a; 342 G4int k = 0; 329 G4int k = 0; 343 G4int nColumns = 2; 330 G4int nColumns = 2; 344 331 345 do 332 do 346 { 333 { 347 in >> a; 334 in >> a; 348 335 349 if (a != -1 && a != -2) 336 if (a != -1 && a != -2) 350 { 337 { 351 if (k%nColumns == 0) 338 if (k%nColumns == 0) 352 { 339 { 353 argEnergies->push_back(a*unitEnergies); 340 argEnergies->push_back(a*unitEnergies); 354 } 341 } 355 else if (k%nColumns == 1) 342 else if (k%nColumns == 1) 356 { 343 { 357 argData->push_back(a*unitData); 344 argData->push_back(a*unitData); 358 } 345 } 359 k++; 346 k++; 360 } 347 } 361 } 348 } 362 while (a != -2); 349 while (a != -2); 363 350 364 SetEnergiesData(argEnergies, argData, 0); 351 SetEnergiesData(argEnergies, argData, 0); 365 352 366 if (randomSet) BuildPdf(); 353 if (randomSet) BuildPdf(); 367 354 368 return true; 355 return true; 369 } 356 } 370 357 371 //....oooOO0OOooo........oooOO0OOooo........oo << 358 372 359 373 G4bool G4EMDataSet::SaveData(const G4String& n 360 G4bool G4EMDataSet::SaveData(const G4String& name) const 374 { 361 { 375 // The file is organized into two columns: 362 // The file is organized into two columns: 376 // 1st column is the energy 363 // 1st column is the energy 377 // 2nd column is the corresponding value 364 // 2nd column is the corresponding value 378 // The file terminates with the pattern: -1 365 // The file terminates with the pattern: -1 -1 379 // -2 366 // -2 -2 380 367 381 G4String fullFileName(FullFileName(name)); 368 G4String fullFileName(FullFileName(name)); 382 std::ofstream out(fullFileName); 369 std::ofstream out(fullFileName); 383 370 384 if (!out.is_open()) 371 if (!out.is_open()) 385 { 372 { 386 G4String message("cannot open \""); 373 G4String message("cannot open \""); 387 message+=fullFileName; 374 message+=fullFileName; 388 message+="\""; 375 message+="\""; 389 G4Exception("G4EMDataSet::SaveData", 376 G4Exception("G4EMDataSet::SaveData", 390 "em1012",FatalException,message); 377 "em1012",FatalException,message); 391 } 378 } 392 379 393 out.precision(10); 380 out.precision(10); 394 out.width(15); 381 out.width(15); 395 out.setf(std::ofstream::left); 382 out.setf(std::ofstream::left); 396 383 397 if (energies!=0 && data!=0) 384 if (energies!=0 && data!=0) 398 { 385 { 399 G4DataVector::const_iterator i(energies- 386 G4DataVector::const_iterator i(energies->begin()); 400 G4DataVector::const_iterator endI(energi 387 G4DataVector::const_iterator endI(energies->end()); 401 G4DataVector::const_iterator j(data->beg 388 G4DataVector::const_iterator j(data->begin()); 402 389 403 while (i!=endI) 390 while (i!=endI) 404 { 391 { 405 out.precision(10); 392 out.precision(10); 406 out.width(15); 393 out.width(15); 407 out.setf(std::ofstream::left); 394 out.setf(std::ofstream::left); 408 out << ((*i)/unitEnergies) << ' '; 395 out << ((*i)/unitEnergies) << ' '; 409 396 410 out.precision(10); 397 out.precision(10); 411 out.width(15); 398 out.width(15); 412 out.setf(std::ofstream::left); 399 out.setf(std::ofstream::left); 413 out << ((*j)/unitData) << std::endl; 400 out << ((*j)/unitData) << std::endl; 414 401 415 i++; 402 i++; 416 j++; 403 j++; 417 } 404 } 418 } 405 } 419 406 420 out.precision(10); 407 out.precision(10); 421 out.width(15); 408 out.width(15); 422 out.setf(std::ofstream::left); 409 out.setf(std::ofstream::left); 423 out << -1.f << ' '; 410 out << -1.f << ' '; 424 411 425 out.precision(10); 412 out.precision(10); 426 out.width(15); 413 out.width(15); 427 out.setf(std::ofstream::left); 414 out.setf(std::ofstream::left); 428 out << -1.f << std::endl; 415 out << -1.f << std::endl; 429 416 430 out.precision(10); 417 out.precision(10); 431 out.width(15); 418 out.width(15); 432 out.setf(std::ofstream::left); 419 out.setf(std::ofstream::left); 433 out << -2.f << ' '; 420 out << -2.f << ' '; 434 421 435 out.precision(10); 422 out.precision(10); 436 out.width(15); 423 out.width(15); 437 out.setf(std::ofstream::left); 424 out.setf(std::ofstream::left); 438 out << -2.f << std::endl; 425 out << -2.f << std::endl; 439 426 440 return true; 427 return true; 441 } 428 } 442 429 443 //....oooOO0OOooo........oooOO0OOooo........oo << 444 430 445 std::size_t G4EMDataSet::FindLowerBound(G4doub << 431 >> 432 size_t G4EMDataSet::FindLowerBound(G4double x) const 446 { 433 { 447 std::size_t lowerBound = 0; << 434 size_t lowerBound = 0; 448 std::size_t upperBound(energies->size() - 1) << 435 size_t upperBound(energies->size() - 1); 449 436 450 while (lowerBound <= upperBound) 437 while (lowerBound <= upperBound) 451 { 438 { 452 std::size_t midBin((lowerBound + upperBo << 439 size_t midBin((lowerBound + upperBound) / 2); 453 440 454 if (x < (*energies)[midBin]) upperBound 441 if (x < (*energies)[midBin]) upperBound = midBin - 1; 455 else lowerBound = midBin + 1; 442 else lowerBound = midBin + 1; 456 } 443 } 457 444 458 return upperBound; 445 return upperBound; 459 } 446 } 460 447 461 //....oooOO0OOooo........oooOO0OOooo........oo << 462 448 463 std::size_t G4EMDataSet::FindLowerBound(G4doub << 449 size_t G4EMDataSet::FindLowerBound(G4double x, G4DataVector* values) const 464 { 450 { 465 std::size_t lowerBound = 0;; << 451 size_t lowerBound = 0;; 466 std::size_t upperBound(values->size() - 1); << 452 size_t upperBound(values->size() - 1); 467 453 468 while (lowerBound <= upperBound) 454 while (lowerBound <= upperBound) 469 { 455 { 470 std::size_t midBin((lowerBound + upperBo << 456 size_t midBin((lowerBound + upperBound) / 2); 471 457 472 if (x < (*values)[midBin]) upperBound = 458 if (x < (*values)[midBin]) upperBound = midBin - 1; 473 else lowerBound = midBin + 1; 459 else lowerBound = midBin + 1; 474 } 460 } 475 461 476 return upperBound; 462 return upperBound; 477 } 463 } 478 464 479 //....oooOO0OOooo........oooOO0OOooo........oo << 480 465 481 G4String G4EMDataSet::FullFileName(const G4Str 466 G4String G4EMDataSet::FullFileName(const G4String& name) const 482 { 467 { 483 const char* path = G4FindDataDir("G4LEDATA") << 468 char* path = getenv("G4LEDATA"); 484 if (!path) { 469 if (!path) { 485 G4Exception("G4EMDataSet::FullFileName", 470 G4Exception("G4EMDataSet::FullFileName", 486 "em0006",FatalException,"G4LEDATA envi 471 "em0006",FatalException,"G4LEDATA environment variable not set"); 487 return ""; 472 return ""; 488 } 473 } 489 std::ostringstream fullFileName; 474 std::ostringstream fullFileName; 490 fullFileName << path << '/' << name << z << 475 fullFileName << path << '/' << name << z << ".dat"; 491 476 492 return G4String(fullFileName.str().c_str()); 477 return G4String(fullFileName.str().c_str()); 493 } 478 } 494 479 495 //....oooOO0OOooo........oooOO0OOooo........oo << 480 496 481 497 void G4EMDataSet::BuildPdf() 482 void G4EMDataSet::BuildPdf() 498 { 483 { 499 pdf = new G4DataVector; 484 pdf = new G4DataVector; 500 G4Integrator <G4EMDataSet, G4double(G4EMData 485 G4Integrator <G4EMDataSet, G4double(G4EMDataSet::*)(G4double)> integrator; 501 486 502 std::size_t nData = data->size(); << 487 G4int nData = data->size(); 503 pdf->push_back(0.); 488 pdf->push_back(0.); 504 489 505 // Integrate the data distribution 490 // Integrate the data distribution 506 std::size_t i; << 491 G4int i; 507 G4double totalSum = 0.; 492 G4double totalSum = 0.; 508 for (i=1; i<nData; ++i) << 493 for (i=1; i<nData; i++) 509 { 494 { 510 G4double xLow = (*energies)[i-1]; 495 G4double xLow = (*energies)[i-1]; 511 G4double xHigh = (*energies)[i]; 496 G4double xHigh = (*energies)[i]; 512 G4double sum = integrator.Legendre96(thi 497 G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh); 513 totalSum = totalSum + sum; 498 totalSum = totalSum + sum; 514 pdf->push_back(totalSum); 499 pdf->push_back(totalSum); 515 } 500 } 516 501 517 // Normalize to the last bin 502 // Normalize to the last bin 518 G4double tot = 0.; 503 G4double tot = 0.; 519 if (totalSum > 0.) tot = 1. / totalSum; 504 if (totalSum > 0.) tot = 1. / totalSum; 520 for (i=1; i<nData; ++i) << 505 for (i=1; i<nData; i++) 521 { 506 { 522 (*pdf)[i] = (*pdf)[i] * tot; 507 (*pdf)[i] = (*pdf)[i] * tot; 523 } 508 } 524 } 509 } 525 510 526 //....oooOO0OOooo........oooOO0OOooo........oo << 527 << 528 G4double G4EMDataSet::RandomSelect(G4int /* co 511 G4double G4EMDataSet::RandomSelect(G4int /* componentId */) const 529 { 512 { 530 G4double value = 0.; 513 G4double value = 0.; 531 // Random select a X value according to the 514 // Random select a X value according to the cumulative probability distribution 532 // derived from the data 515 // derived from the data 533 516 534 if (!pdf) { 517 if (!pdf) { 535 G4Exception("G4EMDataSet::RandomSelect", 518 G4Exception("G4EMDataSet::RandomSelect", 536 "em1012",FatalException,"PDF has not b 519 "em1012",FatalException,"PDF has not been created for this data set"); 537 return value; 520 return value; 538 } 521 } 539 522 540 G4double x = G4UniformRand(); 523 G4double x = G4UniformRand(); 541 524 542 // Locate the random value in the X vector b 525 // Locate the random value in the X vector based on the PDF 543 G4int bin = (G4int)FindLowerBound(x,pdf); << 526 size_t bin = FindLowerBound(x,pdf); 544 527 545 // Interpolate the PDF to calculate the X va 528 // Interpolate the PDF to calculate the X value: 546 // linear interpolation in the first bin (to 529 // linear interpolation in the first bin (to avoid problem with 0), 547 // interpolation with associated data set al 530 // interpolation with associated data set algorithm in other bins >> 531 548 G4LinInterpolation linearAlgo; 532 G4LinInterpolation linearAlgo; 549 if (bin == 0) value = linearAlgo.Calculate(x 533 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies); 550 else value = algorithm->Calculate(x, bin, *p 534 else value = algorithm->Calculate(x, bin, *pdf, *energies); 551 535 >> 536 // G4cout << x << " random bin "<< bin << " - " << value << G4endl; 552 return value; 537 return value; 553 } 538 } 554 539 555 //....oooOO0OOooo........oooOO0OOooo........oo << 556 << 557 G4double G4EMDataSet::IntegrationFunction(G4do 540 G4double G4EMDataSet::IntegrationFunction(G4double x) 558 { 541 { 559 // This function is needed by G4Integrator t 542 // This function is needed by G4Integrator to calculate the integral of the data distribution >> 543 560 G4double y = 0; 544 G4double y = 0; 561 545 562 // Locate the random value in the X vector b << 546 // Locate the random value in the X vector based on the PDF 563 G4int bin = (G4int)FindLowerBound(x); << 547 size_t bin = FindLowerBound(x); 564 548 565 // Interpolate to calculate the X value: 549 // Interpolate to calculate the X value: 566 // linear interpolation in the first bin (to 550 // linear interpolation in the first bin (to avoid problem with 0), 567 // interpolation with associated algorithm i 551 // interpolation with associated algorithm in other bins 568 552 569 G4LinInterpolation linearAlgo; 553 G4LinInterpolation linearAlgo; 570 554 571 if (bin == 0) y = linearAlgo.Calculate(x, bi 555 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data); 572 else y = algorithm->Calculate(x, bin, *energ 556 else y = algorithm->Calculate(x, bin, *energies, *data); 573 557 574 return y; 558 return y; 575 } 559 } 576 560 577 561