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 // 27 // 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch) 29 // 29 // 30 // History: 30 // History: 31 // ----------- 31 // ----------- 32 // 31 Jul 2001 MGP Created 32 // 31 Jul 2001 MGP Created 33 // 33 // 34 // ------------------------------------------- 34 // ------------------------------------------------------------------- 35 35 36 #include "G4RDEMDataSet.hh" 36 #include "G4RDEMDataSet.hh" 37 #include "G4RDVDataSetAlgorithm.hh" 37 #include "G4RDVDataSetAlgorithm.hh" 38 #include <fstream> 38 #include <fstream> 39 #include <sstream> 39 #include <sstream> 40 #include "G4Integrator.hh" 40 #include "G4Integrator.hh" 41 #include "Randomize.hh" 41 #include "Randomize.hh" 42 #include "G4RDLinInterpolation.hh" 42 #include "G4RDLinInterpolation.hh" 43 43 44 44 45 G4RDEMDataSet::G4RDEMDataSet(G4int Z, 45 G4RDEMDataSet::G4RDEMDataSet(G4int Z, 46 G4RDVDataSetAlgorithm* algo, 46 G4RDVDataSetAlgorithm* algo, 47 G4double xUnit, 47 G4double xUnit, 48 G4double yUnit, 48 G4double yUnit, 49 G4bool random): 49 G4bool random): 50 z(Z), 50 z(Z), 51 energies(0), 51 energies(0), 52 data(0), 52 data(0), 53 algorithm(algo), 53 algorithm(algo), 54 unitEnergies(xUnit), 54 unitEnergies(xUnit), 55 unitData(yUnit), 55 unitData(yUnit), 56 pdf(0), 56 pdf(0), 57 randomSet(random) 57 randomSet(random) 58 { 58 { 59 if (algorithm == 0) 59 if (algorithm == 0) 60 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 60 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", 61 "InvalidSetup", FatalException 61 "InvalidSetup", FatalException, "Interpolation == 0!"); 62 if (randomSet) BuildPdf(); 62 if (randomSet) BuildPdf(); 63 } 63 } 64 64 65 G4RDEMDataSet::G4RDEMDataSet(G4int argZ, 65 G4RDEMDataSet::G4RDEMDataSet(G4int argZ, 66 G4DataVector* dataX, 66 G4DataVector* dataX, 67 G4DataVector* dataY, 67 G4DataVector* dataY, 68 G4RDVDataSetAlgorithm* algo, 68 G4RDVDataSetAlgorithm* algo, 69 G4double xUnit, 69 G4double xUnit, 70 G4double yUnit, 70 G4double yUnit, 71 G4bool random): 71 G4bool random): 72 z(argZ), 72 z(argZ), 73 energies(dataX), 73 energies(dataX), 74 data(dataY), 74 data(dataY), 75 algorithm(algo), 75 algorithm(algo), 76 unitEnergies(xUnit), 76 unitEnergies(xUnit), 77 unitData(yUnit), 77 unitData(yUnit), 78 pdf(0), 78 pdf(0), 79 randomSet(random) 79 randomSet(random) 80 { 80 { 81 if (algorithm == 0) 81 if (algorithm == 0) 82 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 82 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", 83 "InvalidSetup", FatalException 83 "InvalidSetup", FatalException, "Interpolation == 0!"); 84 84 85 if ((energies == 0) ^ (data == 0)) 85 if ((energies == 0) ^ (data == 0)) 86 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 86 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup", 87 FatalException, "Different size for 87 FatalException, "Different size for energies and data (zero case)!"); 88 88 89 if (energies == 0) return; 89 if (energies == 0) return; 90 90 91 if (energies->size() != data->size()) 91 if (energies->size() != data->size()) 92 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 92 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup", 93 FatalException, "Different size for 93 FatalException, "Different size for energies and data!"); 94 94 95 if (randomSet) BuildPdf(); 95 if (randomSet) BuildPdf(); 96 } 96 } 97 97 98 G4RDEMDataSet::~G4RDEMDataSet() 98 G4RDEMDataSet::~G4RDEMDataSet() 99 { 99 { 100 delete algorithm; 100 delete algorithm; 101 if (energies) delete energies; 101 if (energies) delete energies; 102 if (data) delete data; 102 if (data) delete data; 103 if (pdf) delete pdf; 103 if (pdf) delete pdf; 104 } 104 } 105 105 106 G4double G4RDEMDataSet::FindValue(G4double ene 106 G4double G4RDEMDataSet::FindValue(G4double energy, G4int /* componentId */) const 107 { 107 { 108 if (!energies) 108 if (!energies) 109 G4Exception("G4RDEMDataSet::FindValue()", 109 G4Exception("G4RDEMDataSet::FindValue()", "InvalidSetup", 110 FatalException, "Energies == 0 110 FatalException, "Energies == 0!"); 111 if (energies->empty()) return 0; 111 if (energies->empty()) return 0; 112 if (energy <= (*energies)[0]) return (*data) 112 if (energy <= (*energies)[0]) return (*data)[0]; 113 113 114 size_t i = energies->size()-1; 114 size_t i = energies->size()-1; 115 if (energy >= (*energies)[i]) return (*data) 115 if (energy >= (*energies)[i]) return (*data)[i]; 116 116 117 return algorithm->Calculate(energy, FindLowe 117 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data); 118 } 118 } 119 119 120 120 121 void G4RDEMDataSet::PrintData(void) const 121 void G4RDEMDataSet::PrintData(void) const 122 { 122 { 123 if (!energies) 123 if (!energies) 124 { 124 { 125 G4cout << "Data not available." << G4end 125 G4cout << "Data not available." << G4endl; 126 } 126 } 127 else 127 else 128 { 128 { 129 size_t size = energies->size(); 129 size_t size = energies->size(); 130 for (size_t i(0); i<size; i++) 130 for (size_t i(0); i<size; i++) 131 { 131 { 132 G4cout << "Point: " << ((*energies)[i]/uni 132 G4cout << "Point: " << ((*energies)[i]/unitEnergies) 133 << " - Data value: " << ((*data)[i]/unitD 133 << " - Data value: " << ((*data)[i]/unitData); 134 if (pdf != 0) G4cout << " - PDF : " << (*p 134 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i]; 135 G4cout << G4endl; 135 G4cout << G4endl; 136 } 136 } 137 } 137 } 138 } 138 } 139 139 140 140 141 void G4RDEMDataSet::SetEnergiesData(G4DataVect 141 void G4RDEMDataSet::SetEnergiesData(G4DataVector* dataX, 142 G4DataVector* dataY, 142 G4DataVector* dataY, 143 G4int /* componentId */) 143 G4int /* componentId */) 144 { 144 { 145 if (energies) delete energies; 145 if (energies) delete energies; 146 energies = dataX; 146 energies = dataX; 147 147 148 if (data) delete data; 148 if (data) delete data; 149 data = dataY; 149 data = dataY; 150 150 151 if ((energies == 0) ^ (data==0)) 151 if ((energies == 0) ^ (data==0)) 152 G4Exception("G4RDEMDataSet::SetEnergiesDat 152 G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup", 153 FatalException, "Different size for ene 153 FatalException, "Different size for energies and data (zero case)!"); 154 154 155 if (energies == 0) return; 155 if (energies == 0) return; 156 156 157 if (energies->size() != data->size()) 157 if (energies->size() != data->size()) 158 G4Exception("G4RDEMDataSet::SetEnergiesDat 158 G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup", 159 FatalException, "Different size for ene 159 FatalException, "Different size for energies and data!"); 160 } 160 } 161 161 162 G4bool G4RDEMDataSet::LoadData(const G4String& 162 G4bool G4RDEMDataSet::LoadData(const G4String& fileName) 163 { 163 { 164 // The file is organized into two columns: 164 // The file is organized into two columns: 165 // 1st column is the energy 165 // 1st column is the energy 166 // 2nd column is the corresponding value 166 // 2nd column is the corresponding value 167 // The file terminates with the pattern: -1 167 // The file terminates with the pattern: -1 -1 168 // -2 168 // -2 -2 169 169 170 G4String fullFileName(FullFileName(fileName) 170 G4String fullFileName(FullFileName(fileName)); 171 std::ifstream in(fullFileName); 171 std::ifstream in(fullFileName); 172 172 173 if (!in.is_open()) 173 if (!in.is_open()) 174 { 174 { 175 G4String message("Data file \""); 175 G4String message("Data file \""); 176 message += fullFileName; 176 message += fullFileName; 177 message += "\" not found"; 177 message += "\" not found"; 178 G4Exception("G4RDEMDataSet::LoadData()", 178 G4Exception("G4RDEMDataSet::LoadData()", "DataNotFound", 179 FatalException, message); 179 FatalException, message); 180 } 180 } 181 181 182 G4DataVector* argEnergies=new G4DataVector; 182 G4DataVector* argEnergies=new G4DataVector; 183 G4DataVector* argData=new G4DataVector; 183 G4DataVector* argData=new G4DataVector; 184 184 185 G4double a; 185 G4double a; 186 bool energyColumn(true); 186 bool energyColumn(true); 187 187 188 do 188 do 189 { 189 { 190 in >> a; 190 in >> a; 191 191 192 if (a!=-1 && a!=-2) 192 if (a!=-1 && a!=-2) 193 { 193 { 194 if (energyColumn) 194 if (energyColumn) 195 argEnergies->push_back(a*unitEnergies); 195 argEnergies->push_back(a*unitEnergies); 196 else 196 else 197 argData->push_back(a*unitData); 197 argData->push_back(a*unitData); 198 energyColumn=(!energyColumn); 198 energyColumn=(!energyColumn); 199 } 199 } 200 } 200 } 201 while (a != -2); 201 while (a != -2); 202 202 203 SetEnergiesData(argEnergies, argData, 0); 203 SetEnergiesData(argEnergies, argData, 0); 204 if (randomSet) BuildPdf(); 204 if (randomSet) BuildPdf(); 205 205 206 return true; 206 return true; 207 } 207 } 208 208 209 G4bool G4RDEMDataSet::SaveData(const G4String& 209 G4bool G4RDEMDataSet::SaveData(const G4String& name) const 210 { 210 { 211 // The file is organized into two columns: 211 // The file is organized into two columns: 212 // 1st column is the energy 212 // 1st column is the energy 213 // 2nd column is the corresponding value 213 // 2nd column is the corresponding value 214 // The file terminates with the pattern: -1 214 // The file terminates with the pattern: -1 -1 215 // -2 215 // -2 -2 216 216 217 G4String fullFileName(FullFileName(name)); 217 G4String fullFileName(FullFileName(name)); 218 std::ofstream out(fullFileName); 218 std::ofstream out(fullFileName); 219 219 220 if (!out.is_open()) 220 if (!out.is_open()) 221 { 221 { 222 G4String message("Cannot open \""); 222 G4String message("Cannot open \""); 223 message+=fullFileName; 223 message+=fullFileName; 224 message+="\""; 224 message+="\""; 225 G4Exception("G4RDEMDataSet::SaveData()", 225 G4Exception("G4RDEMDataSet::SaveData()", "CannotOpenFile", 226 FatalException, message); 226 FatalException, message); 227 } 227 } 228 228 229 out.precision(10); 229 out.precision(10); 230 out.width(15); 230 out.width(15); 231 out.setf(std::ofstream::left); 231 out.setf(std::ofstream::left); 232 232 233 if (energies!=0 && data!=0) 233 if (energies!=0 && data!=0) 234 { 234 { 235 G4DataVector::const_iterator i(energies- 235 G4DataVector::const_iterator i(energies->begin()); 236 G4DataVector::const_iterator endI(energi 236 G4DataVector::const_iterator endI(energies->end()); 237 G4DataVector::const_iterator j(data->beg 237 G4DataVector::const_iterator j(data->begin()); 238 238 239 while (i!=endI) 239 while (i!=endI) 240 { 240 { 241 out.precision(10); 241 out.precision(10); 242 out.width(15); 242 out.width(15); 243 out.setf(std::ofstream::left); 243 out.setf(std::ofstream::left); 244 out << ((*i)/unitEnergies) << ' '; 244 out << ((*i)/unitEnergies) << ' '; 245 245 246 out.precision(10); 246 out.precision(10); 247 out.width(15); 247 out.width(15); 248 out.setf(std::ofstream::left); 248 out.setf(std::ofstream::left); 249 out << ((*j)/unitData) << std::endl; 249 out << ((*j)/unitData) << std::endl; 250 250 251 i++; 251 i++; 252 j++; 252 j++; 253 } 253 } 254 } 254 } 255 255 256 out.precision(10); 256 out.precision(10); 257 out.width(15); 257 out.width(15); 258 out.setf(std::ofstream::left); 258 out.setf(std::ofstream::left); 259 out << -1.f << ' '; 259 out << -1.f << ' '; 260 260 261 out.precision(10); 261 out.precision(10); 262 out.width(15); 262 out.width(15); 263 out.setf(std::ofstream::left); 263 out.setf(std::ofstream::left); 264 out << -1.f << std::endl; 264 out << -1.f << std::endl; 265 265 266 out.precision(10); 266 out.precision(10); 267 out.width(15); 267 out.width(15); 268 out.setf(std::ofstream::left); 268 out.setf(std::ofstream::left); 269 out << -2.f << ' '; 269 out << -2.f << ' '; 270 270 271 out.precision(10); 271 out.precision(10); 272 out.width(15); 272 out.width(15); 273 out.setf(std::ofstream::left); 273 out.setf(std::ofstream::left); 274 out << -2.f << std::endl; 274 out << -2.f << std::endl; 275 275 276 return true; 276 return true; 277 } 277 } 278 278 279 size_t G4RDEMDataSet::FindLowerBound(G4double 279 size_t G4RDEMDataSet::FindLowerBound(G4double x) const 280 { 280 { 281 size_t lowerBound = 0; 281 size_t lowerBound = 0; 282 size_t upperBound(energies->size() - 1); 282 size_t upperBound(energies->size() - 1); 283 283 284 while (lowerBound <= upperBound) 284 while (lowerBound <= upperBound) 285 { 285 { 286 size_t midBin((lowerBound + upperBound) 286 size_t midBin((lowerBound + upperBound) / 2); 287 287 288 if (x < (*energies)[midBin]) upperBound 288 if (x < (*energies)[midBin]) upperBound = midBin - 1; 289 else lowerBound = midBin + 1; 289 else lowerBound = midBin + 1; 290 } 290 } 291 291 292 return upperBound; 292 return upperBound; 293 } 293 } 294 294 295 295 296 size_t G4RDEMDataSet::FindLowerBound(G4double 296 size_t G4RDEMDataSet::FindLowerBound(G4double x, G4DataVector* values) const 297 { 297 { 298 size_t lowerBound = 0;; 298 size_t lowerBound = 0;; 299 size_t upperBound(values->size() - 1); 299 size_t upperBound(values->size() - 1); 300 300 301 while (lowerBound <= upperBound) 301 while (lowerBound <= upperBound) 302 { 302 { 303 size_t midBin((lowerBound + upperBound) 303 size_t midBin((lowerBound + upperBound) / 2); 304 304 305 if (x < (*values)[midBin]) upperBound = 305 if (x < (*values)[midBin]) upperBound = midBin - 1; 306 else lowerBound = midBin + 1; 306 else lowerBound = midBin + 1; 307 } 307 } 308 308 309 return upperBound; 309 return upperBound; 310 } 310 } 311 311 312 312 313 G4String G4RDEMDataSet::FullFileName(const G4S 313 G4String G4RDEMDataSet::FullFileName(const G4String& name) const 314 { 314 { 315 const char* path = G4FindDataDir("G4LEDATA") << 315 char* path = std::getenv("G4LEDATA"); 316 if (!path) 316 if (!path) 317 G4Exception("G4RDEMDataSet::FullFileName() 317 G4Exception("G4RDEMDataSet::FullFileName()", "InvalidSetup", 318 FatalException, "G4LEDATA envi 318 FatalException, "G4LEDATA environment variable not set!"); 319 319 320 std::ostringstream fullFileName; 320 std::ostringstream fullFileName; 321 fullFileName << path << '/' << name << z << 321 fullFileName << path << '/' << name << z << ".dat"; 322 322 323 return G4String(fullFileName.str().c_str()); 323 return G4String(fullFileName.str().c_str()); 324 } 324 } 325 325 326 326 327 void G4RDEMDataSet::BuildPdf() 327 void G4RDEMDataSet::BuildPdf() 328 { 328 { 329 pdf = new G4DataVector; 329 pdf = new G4DataVector; 330 G4Integrator <G4RDEMDataSet, G4double(G4RDEM 330 G4Integrator <G4RDEMDataSet, G4double(G4RDEMDataSet::*)(G4double)> integrator; 331 331 332 G4int nData = data->size(); 332 G4int nData = data->size(); 333 pdf->push_back(0.); 333 pdf->push_back(0.); 334 334 335 // Integrate the data distribution 335 // Integrate the data distribution 336 G4int i; 336 G4int i; 337 G4double totalSum = 0.; 337 G4double totalSum = 0.; 338 for (i=1; i<nData; i++) 338 for (i=1; i<nData; i++) 339 { 339 { 340 G4double xLow = (*energies)[i-1]; 340 G4double xLow = (*energies)[i-1]; 341 G4double xHigh = (*energies)[i]; 341 G4double xHigh = (*energies)[i]; 342 G4double sum = integrator.Legendre96(thi 342 G4double sum = integrator.Legendre96(this, &G4RDEMDataSet::IntegrationFunction, xLow, xHigh); 343 totalSum = totalSum + sum; 343 totalSum = totalSum + sum; 344 pdf->push_back(totalSum); 344 pdf->push_back(totalSum); 345 } 345 } 346 346 347 // Normalize to the last bin 347 // Normalize to the last bin 348 G4double tot = 0.; 348 G4double tot = 0.; 349 if (totalSum > 0.) tot = 1. / totalSum; 349 if (totalSum > 0.) tot = 1. / totalSum; 350 for (i=1; i<nData; i++) 350 for (i=1; i<nData; i++) 351 { 351 { 352 (*pdf)[i] = (*pdf)[i] * tot; 352 (*pdf)[i] = (*pdf)[i] * tot; 353 } 353 } 354 } 354 } 355 355 356 356 357 G4double G4RDEMDataSet::RandomSelect(G4int /* 357 G4double G4RDEMDataSet::RandomSelect(G4int /* componentId */) const 358 { 358 { 359 // Random select a X value according to the 359 // Random select a X value according to the cumulative probability distribution 360 // derived from the data 360 // derived from the data 361 361 362 if (!pdf) G4Exception("G4RDEMDataSet::Random 362 if (!pdf) G4Exception("G4RDEMDataSet::RandomSelect()", "InvalidSetup", 363 FatalException, "PDF has not been c 363 FatalException, "PDF has not been created for this data set"); 364 364 365 G4double value = 0.; 365 G4double value = 0.; 366 G4double x = G4UniformRand(); 366 G4double x = G4UniformRand(); 367 367 368 // Locate the random value in the X vector b 368 // Locate the random value in the X vector based on the PDF 369 size_t bin = FindLowerBound(x,pdf); 369 size_t bin = FindLowerBound(x,pdf); 370 370 371 // Interpolate the PDF to calculate the X va 371 // Interpolate the PDF to calculate the X value: 372 // linear interpolation in the first bin (to 372 // linear interpolation in the first bin (to avoid problem with 0), 373 // interpolation with associated data set al 373 // interpolation with associated data set algorithm in other bins 374 374 375 G4RDLinInterpolation linearAlgo; 375 G4RDLinInterpolation linearAlgo; 376 if (bin == 0) value = linearAlgo.Calculate(x 376 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies); 377 else value = algorithm->Calculate(x, bin, *p 377 else value = algorithm->Calculate(x, bin, *pdf, *energies); 378 378 379 // G4cout << x << " random bin "<< bin << " 379 // G4cout << x << " random bin "<< bin << " - " << value << G4endl; 380 return value; 380 return value; 381 } 381 } 382 382 383 G4double G4RDEMDataSet::IntegrationFunction(G4 383 G4double G4RDEMDataSet::IntegrationFunction(G4double x) 384 { 384 { 385 // This function is needed by G4Integrator t 385 // This function is needed by G4Integrator to calculate the integral of the data distribution 386 386 387 G4double y = 0; 387 G4double y = 0; 388 388 389 // Locate the random value in the X vector ba 389 // Locate the random value in the X vector based on the PDF 390 size_t bin = FindLowerBound(x); 390 size_t bin = FindLowerBound(x); 391 391 392 // Interpolate to calculate the X value: 392 // Interpolate to calculate the X value: 393 // linear interpolation in the first bin (to 393 // linear interpolation in the first bin (to avoid problem with 0), 394 // interpolation with associated algorithm i 394 // interpolation with associated algorithm in other bins 395 395 396 G4RDLinInterpolation linearAlgo; 396 G4RDLinInterpolation linearAlgo; 397 397 398 if (bin == 0) y = linearAlgo.Calculate(x, bi 398 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data); 399 else y = algorithm->Calculate(x, bin, *energ 399 else y = algorithm->Calculate(x, bin, *energies, *data); 400 400 401 return y; 401 return y; 402 } 402 } 403 403 404 404