00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "G4EMDataSet.hh"
00054 #include "G4VDataSetAlgorithm.hh"
00055 #include <fstream>
00056 #include <sstream>
00057 #include "G4Integrator.hh"
00058 #include "Randomize.hh"
00059 #include "G4LinInterpolation.hh"
00060
00061
00062 G4EMDataSet::G4EMDataSet(G4int Z,
00063 G4VDataSetAlgorithm* algo,
00064 G4double xUnit,
00065 G4double yUnit,
00066 G4bool random):
00067 z(Z),
00068 energies(0),
00069 data(0),
00070 log_energies(0),
00071 log_data(0),
00072 algorithm(algo),
00073 unitEnergies(xUnit),
00074 unitData(yUnit),
00075 pdf(0),
00076 randomSet(random)
00077 {
00078 if (algorithm == 0) {
00079 G4Exception("G4EMDataSet::G4EMDataSet",
00080 "em1012",FatalException,"interpolation == 0");
00081 } else if (randomSet) { BuildPdf(); }
00082 }
00083
00084 G4EMDataSet::G4EMDataSet(G4int argZ,
00085 G4DataVector* dataX,
00086 G4DataVector* dataY,
00087 G4VDataSetAlgorithm* algo,
00088 G4double xUnit,
00089 G4double yUnit,
00090 G4bool random):
00091 z(argZ),
00092 energies(dataX),
00093 data(dataY),
00094 log_energies(0),
00095 log_data(0),
00096 algorithm(algo),
00097 unitEnergies(xUnit),
00098 unitData(yUnit),
00099 pdf(0),
00100 randomSet(random)
00101 {
00102 if (algorithm == 0) {
00103 G4Exception("G4EMDataSet::G4EMDataSet",
00104 "em1012",FatalException,"interpolation == 0");
00105 } else if ((energies == 0) ^ (data == 0)) {
00106 G4Exception("G4EMDataSet::G4EMDataSet",
00107 "em1012",FatalException,"different size for energies and data (zero case)");
00108 } else if (energies->size() != data->size()) {
00109 G4Exception("G4EMDataSet::G4EMDataSet",
00110 "em1012",FatalException,"different size for energies and data");
00111 } else if (randomSet) {
00112 BuildPdf();
00113 }
00114 }
00115
00116 G4EMDataSet::G4EMDataSet(G4int argZ,
00117 G4DataVector* dataX,
00118 G4DataVector* dataY,
00119 G4DataVector* dataLogX,
00120 G4DataVector* dataLogY,
00121 G4VDataSetAlgorithm* algo,
00122 G4double xUnit,
00123 G4double yUnit,
00124 G4bool random):
00125 z(argZ),
00126 energies(dataX),
00127 data(dataY),
00128 log_energies(dataLogX),
00129 log_data(dataLogY),
00130 algorithm(algo),
00131 unitEnergies(xUnit),
00132 unitData(yUnit),
00133 pdf(0),
00134 randomSet(random)
00135 {
00136 if (algorithm == 0) {
00137 G4Exception("G4EMDataSet::G4EMDataSet",
00138 "em1012",FatalException,"interpolation == 0");
00139 } else if ((energies == 0) ^ (data == 0)) {
00140 G4Exception("G4EMDataSet::G4EMDataSet",
00141 "em1012",FatalException,"different size for energies and data (zero case)");
00142 } else if (energies->size() != data->size()) {
00143 G4Exception("G4EMDataSet::G4EMDataSet",
00144 "em1012",FatalException,"different size for energies and data");
00145 } else if ((log_energies == 0) ^ (log_data == 0)) {
00146 G4Exception("G4EMDataSet::G4EMDataSet",
00147 "em1012",FatalException,"different size for log energies and log data (zero case)");
00148 } else if (log_energies->size() != log_data->size()) {
00149 G4Exception("G4EMDataSet::G4EMDataSet",
00150 "em1012",FatalException,"different size for log energies and log data");
00151 } else if (randomSet) {
00152 BuildPdf();
00153 }
00154 }
00155
00156
00157 G4EMDataSet::~G4EMDataSet()
00158 {
00159 delete algorithm;
00160 delete energies;
00161 delete data;
00162 delete pdf;
00163 delete log_energies;
00164 delete log_data;
00165 }
00166
00167
00168 G4double G4EMDataSet::FindValue(G4double energy, G4int ) const
00169 {
00170 if (!energies) {
00171 G4Exception("G4EMDataSet::FindValue",
00172 "em1012",FatalException,"energies == 0");
00173 return 0.0;
00174 } else if (energies->empty()) {
00175 return 0.0;
00176 } else if (energy <= (*energies)[0]) {
00177 return (*data)[0];
00178 }
00179 size_t i = energies->size()-1;
00180 if (energy >= (*energies)[i]) { return (*data)[i]; }
00181
00182
00183
00184 if (log_energies != 0)
00185 {
00186 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data, *log_energies, *log_data);
00187 }
00188
00189 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
00190 }
00191
00192
00193 void G4EMDataSet::PrintData(void) const
00194 {
00195 if (!energies)
00196 {
00197 G4cout << "Data not available." << G4endl;
00198 }
00199 else
00200 {
00201 size_t size = energies->size();
00202 for (size_t i(0); i<size; i++)
00203 {
00204 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
00205 << " - Data value: " << ((*data)[i]/unitData);
00206 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
00207 G4cout << G4endl;
00208 }
00209 }
00210 }
00211
00212
00213 void G4EMDataSet::SetEnergiesData(G4DataVector* dataX,
00214 G4DataVector* dataY,
00215 G4int )
00216 {
00217 if (energies) { delete energies; }
00218 energies = dataX;
00219
00220 if (data) { delete data; }
00221 data = dataY;
00222
00223 if ((energies == 0) ^ (data==0)) {
00224 G4Exception("G4EMDataSet::SetEnergiesData",
00225 "em1012",FatalException,"different size for energies and data (zero case)");
00226 return;
00227 } else if (energies == 0) { return; }
00228
00229
00230 if (energies->size() != data->size()) {
00231 G4Exception("G4EMDataSet::SetEnergiesData",
00232 "em1012",FatalException,"different size for energies and data");
00233 }
00234 }
00235
00236 void G4EMDataSet::SetLogEnergiesData(G4DataVector* dataX,
00237 G4DataVector* dataY,
00238 G4DataVector* data_logX,
00239 G4DataVector* data_logY,
00240 G4int )
00241 {
00242
00243 if (energies) { delete energies; }
00244 energies = dataX;
00245 if (data) { delete data; }
00246 data = dataY;
00247
00248 if (log_energies) { delete log_energies; }
00249 log_energies = data_logX;
00250 if (log_data) { delete log_data; }
00251 log_data = data_logY;
00252
00253
00254 if ( !energies ) {
00255 if(data || log_energies || log_data ) {
00256 G4Exception("G4EMDataSet::SetLogEnergiesData",
00257 "em1012",FatalException,"inconsistent data");
00258 }
00259 return;
00260 } else {
00261 if ( !data ) {
00262 G4Exception("G4EMDataSet::SetLogEnergiesData",
00263 "em1012",FatalException,"only energy, no data");
00264 return;
00265 } else if (energies->size() != data->size()) {
00266 G4Exception("G4EMDataSet::SetLogEnergiesData",
00267 "em1012",FatalException,"different size for energies and data");
00268 return;
00269 }
00270
00271
00272
00273 if ( !log_energies ) {
00274 if(log_data) {
00275 G4Exception("G4EMDataSet::SetLogEnergiesData",
00276 "em1012",FatalException,"inconsistence of log_data");
00277 }
00278 return;
00279 } else {
00280 if ( !log_data ) {
00281 G4Exception("G4EMDataSet::SetLogEnergiesData",
00282 "em1012",FatalException,"only log_energies, no data");
00283 } else if ((log_energies->size() != log_data->size()) || (log_energies->size() != data->size())) {
00284 G4Exception("G4EMDataSet::SetLogEnergiesData",
00285 "em1012",FatalException,"different size for log energies and data");
00286 }
00287 }
00288 }
00289
00290 }
00291
00292 G4bool G4EMDataSet::LoadData(const G4String& fileName)
00293 {
00294
00295
00296
00297
00298
00299
00300 G4String fullFileName(FullFileName(fileName));
00301 std::ifstream in(fullFileName);
00302
00303 if (!in.is_open())
00304 {
00305 G4String message("data file \"");
00306 message += fullFileName;
00307 message += "\" not found";
00308 G4Exception("G4EMDataSet::LoadData",
00309 "em1012",FatalException,message);
00310 return false;
00311 }
00312
00313 delete energies;
00314 delete data;
00315 delete log_energies;
00316 delete log_data;
00317 energies = new G4DataVector;
00318 data = new G4DataVector;
00319 log_energies = new G4DataVector;
00320 log_data = new G4DataVector;
00321
00322 G4double a, b;
00323
00324
00325
00326 do
00327 {
00328 in >> a >> b;
00329
00330 if (a != -1 && a != -2)
00331 {
00332 if (a==0.) { a=1e-300; }
00333 if (b==0.) { b=1e-300; }
00334 a *= unitEnergies;
00335 b *= unitData;
00336 energies->push_back(a);
00337 log_energies->push_back(std::log10(a));
00338 data->push_back(b);
00339 log_data->push_back(std::log10(b));
00340 }
00341 }
00342 while (a != -2);
00343
00344 if (randomSet) { BuildPdf(); }
00345
00346 return true;
00347 }
00348
00349
00350 G4bool G4EMDataSet::LoadNonLogData(const G4String& fileName)
00351 {
00352
00353
00354
00355
00356
00357
00358 G4String fullFileName(FullFileName(fileName));
00359 std::ifstream in(fullFileName);
00360 if (!in.is_open())
00361 {
00362 G4String message("data file \"");
00363 message += fullFileName;
00364 message += "\" not found";
00365 G4Exception("G4EMDataSet::LoadNonLogData",
00366 "em1012",FatalException,message);
00367 }
00368
00369 G4DataVector* argEnergies=new G4DataVector;
00370 G4DataVector* argData=new G4DataVector;
00371
00372 G4double a;
00373 G4int k = 0;
00374 G4int nColumns = 2;
00375
00376 do
00377 {
00378 in >> a;
00379
00380 if (a != -1 && a != -2)
00381 {
00382 if (k%nColumns == 0)
00383 {
00384 argEnergies->push_back(a*unitEnergies);
00385 }
00386 else if (k%nColumns == 1)
00387 {
00388 argData->push_back(a*unitData);
00389 }
00390 k++;
00391 }
00392 }
00393 while (a != -2);
00394
00395 SetEnergiesData(argEnergies, argData, 0);
00396
00397 if (randomSet) BuildPdf();
00398
00399 return true;
00400 }
00401
00402
00403
00404 G4bool G4EMDataSet::SaveData(const G4String& name) const
00405 {
00406
00407
00408
00409
00410
00411
00412 G4String fullFileName(FullFileName(name));
00413 std::ofstream out(fullFileName);
00414
00415 if (!out.is_open())
00416 {
00417 G4String message("cannot open \"");
00418 message+=fullFileName;
00419 message+="\"";
00420 G4Exception("G4EMDataSet::SaveData",
00421 "em1012",FatalException,message);
00422 }
00423
00424 out.precision(10);
00425 out.width(15);
00426 out.setf(std::ofstream::left);
00427
00428 if (energies!=0 && data!=0)
00429 {
00430 G4DataVector::const_iterator i(energies->begin());
00431 G4DataVector::const_iterator endI(energies->end());
00432 G4DataVector::const_iterator j(data->begin());
00433
00434 while (i!=endI)
00435 {
00436 out.precision(10);
00437 out.width(15);
00438 out.setf(std::ofstream::left);
00439 out << ((*i)/unitEnergies) << ' ';
00440
00441 out.precision(10);
00442 out.width(15);
00443 out.setf(std::ofstream::left);
00444 out << ((*j)/unitData) << std::endl;
00445
00446 i++;
00447 j++;
00448 }
00449 }
00450
00451 out.precision(10);
00452 out.width(15);
00453 out.setf(std::ofstream::left);
00454 out << -1.f << ' ';
00455
00456 out.precision(10);
00457 out.width(15);
00458 out.setf(std::ofstream::left);
00459 out << -1.f << std::endl;
00460
00461 out.precision(10);
00462 out.width(15);
00463 out.setf(std::ofstream::left);
00464 out << -2.f << ' ';
00465
00466 out.precision(10);
00467 out.width(15);
00468 out.setf(std::ofstream::left);
00469 out << -2.f << std::endl;
00470
00471 return true;
00472 }
00473
00474
00475
00476 size_t G4EMDataSet::FindLowerBound(G4double x) const
00477 {
00478 size_t lowerBound = 0;
00479 size_t upperBound(energies->size() - 1);
00480
00481 while (lowerBound <= upperBound)
00482 {
00483 size_t midBin((lowerBound + upperBound) / 2);
00484
00485 if (x < (*energies)[midBin]) upperBound = midBin - 1;
00486 else lowerBound = midBin + 1;
00487 }
00488
00489 return upperBound;
00490 }
00491
00492
00493 size_t G4EMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
00494 {
00495 size_t lowerBound = 0;;
00496 size_t upperBound(values->size() - 1);
00497
00498 while (lowerBound <= upperBound)
00499 {
00500 size_t midBin((lowerBound + upperBound) / 2);
00501
00502 if (x < (*values)[midBin]) upperBound = midBin - 1;
00503 else lowerBound = midBin + 1;
00504 }
00505
00506 return upperBound;
00507 }
00508
00509
00510 G4String G4EMDataSet::FullFileName(const G4String& name) const
00511 {
00512 char* path = getenv("G4LEDATA");
00513 if (!path) {
00514 G4Exception("G4EMDataSet::FullFileName",
00515 "em0006",FatalException,"G4LEDATA environment variable not set");
00516 return "";
00517 }
00518 std::ostringstream fullFileName;
00519 fullFileName << path << '/' << name << z << ".dat";
00520
00521 return G4String(fullFileName.str().c_str());
00522 }
00523
00524
00525
00526 void G4EMDataSet::BuildPdf()
00527 {
00528 pdf = new G4DataVector;
00529 G4Integrator <G4EMDataSet, G4double(G4EMDataSet::*)(G4double)> integrator;
00530
00531 G4int nData = data->size();
00532 pdf->push_back(0.);
00533
00534
00535 G4int i;
00536 G4double totalSum = 0.;
00537 for (i=1; i<nData; i++)
00538 {
00539 G4double xLow = (*energies)[i-1];
00540 G4double xHigh = (*energies)[i];
00541 G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh);
00542 totalSum = totalSum + sum;
00543 pdf->push_back(totalSum);
00544 }
00545
00546
00547 G4double tot = 0.;
00548 if (totalSum > 0.) tot = 1. / totalSum;
00549 for (i=1; i<nData; i++)
00550 {
00551 (*pdf)[i] = (*pdf)[i] * tot;
00552 }
00553 }
00554
00555 G4double G4EMDataSet::RandomSelect(G4int ) const
00556 {
00557 G4double value = 0.;
00558
00559
00560
00561 if (!pdf) {
00562 G4Exception("G4EMDataSet::RandomSelect",
00563 "em1012",FatalException,"PDF has not been created for this data set");
00564 return value;
00565 }
00566
00567 G4double x = G4UniformRand();
00568
00569
00570 size_t bin = FindLowerBound(x,pdf);
00571
00572
00573
00574
00575
00576 G4LinInterpolation linearAlgo;
00577 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
00578 else value = algorithm->Calculate(x, bin, *pdf, *energies);
00579
00580
00581 return value;
00582 }
00583
00584 G4double G4EMDataSet::IntegrationFunction(G4double x)
00585 {
00586
00587
00588 G4double y = 0;
00589
00590
00591 size_t bin = FindLowerBound(x);
00592
00593
00594
00595
00596
00597 G4LinInterpolation linearAlgo;
00598
00599 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
00600 else y = algorithm->Calculate(x, bin, *energies, *data);
00601
00602 return y;
00603 }
00604