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