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 #include "G4DataSet.hh"
00039 #include "G4IInterpolator.hh"
00040 #include <fstream>
00041 #include <sstream>
00042 #include "G4Integrator.hh"
00043 #include "Randomize.hh"
00044 #include "G4LinInterpolation.hh"
00045
00046
00047 G4DataSet::G4DataSet(G4int Z,
00048 G4IInterpolator* algo,
00049 G4double xUnit,
00050 G4double yUnit,
00051 G4bool random):
00052 z(Z),
00053 energies(0),
00054 data(0),
00055 algorithm(algo),
00056 unitEnergies(xUnit),
00057 unitData(yUnit),
00058 pdf(0),
00059 randomSet(random)
00060 {
00061 if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
00062 "pii00000101",
00063 FatalException,
00064 "Interpolation == 0");
00065 if (randomSet) BuildPdf();
00066 }
00067
00068 G4DataSet::G4DataSet(G4int argZ,
00069 G4DataVector* dataX,
00070 G4DataVector* dataY,
00071 G4IInterpolator* algo,
00072 G4double xUnit,
00073 G4double yUnit,
00074 G4bool random):
00075 z(argZ),
00076 energies(dataX),
00077 data(dataY),
00078 algorithm(algo),
00079 unitEnergies(xUnit),
00080 unitData(yUnit),
00081 pdf(0),
00082 randomSet(random)
00083 {
00084 if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
00085 "pii00000110",
00086 FatalException,
00087 "Interpolation == 0");
00088
00089 if ((energies == 0) ^ (data == 0))
00090 G4Exception("G4DataSet::G4DataSet",
00091 "pii00000111-",
00092 FatalException,
00093 "different size for energies and data (zero case)");
00094
00095 if (energies == 0) return;
00096
00097 if (energies->size() != data->size())
00098 G4Exception("G4DataSet::G4DataSet",
00099 "pii00000112",
00100 FatalException,
00101 "different size for energies and data");
00102
00103 if (randomSet) BuildPdf();
00104 }
00105
00106 G4DataSet::~G4DataSet()
00107 {
00108 delete algorithm;
00109 if (energies) delete energies;
00110 if (data) delete data;
00111 if (pdf) delete pdf;
00112 }
00113
00114 G4double G4DataSet::FindValue(G4double energy, G4int ) const
00115 {
00116 if (!energies) G4Exception("G4DataSet::FindValue",
00117 "pii00000120",
00118 FatalException,
00119 "energies == 0");
00120 if (energies->empty()) return 0;
00121 if (energy <= (*energies)[0]) return (*data)[0];
00122
00123 size_t i = energies->size()-1;
00124 if (energy >= (*energies)[i]) return (*data)[i];
00125
00126 G4double interpolated = algorithm->Calculate(energy,FindLowerBound(energy),*energies,*data);
00127 return interpolated;
00128 }
00129
00130
00131 void G4DataSet::PrintData(void) const
00132 {
00133 if (!energies)
00134 {
00135 G4cout << "Data not available." << G4endl;
00136 }
00137 else
00138 {
00139 size_t size = energies->size();
00140 for (size_t i(0); i<size; i++)
00141 {
00142 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
00143 << " - Data value: " << ((*data)[i]/unitData);
00144 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
00145 G4cout << G4endl;
00146 }
00147 }
00148 }
00149
00150
00151 void G4DataSet::SetEnergiesData(G4DataVector* dataX,
00152 G4DataVector* dataY,
00153 G4int )
00154 {
00155 if (energies) delete energies;
00156 energies = dataX;
00157
00158 if (data) delete data;
00159 data = dataY;
00160
00161 if ((energies == 0) ^ (data==0))
00162 G4Exception("G4DataSet::SetEnergiesData",
00163 "pii00000130",
00164 FatalException,
00165 "different size for energies and data (zero case)");
00166
00167 if (energies == 0) return;
00168
00169 if (energies->size() != data->size())
00170 G4Exception("G4DataSet::SetEnergiesData",
00171 "pii00000131",
00172 FatalException,
00173 "different size for energies and data");
00174 }
00175
00176 G4bool G4DataSet::LoadData(const G4String& fileName)
00177 {
00178
00179
00180
00181
00182
00183
00184 G4String fullFileName(FullFileName(fileName));
00185 std::ifstream in(fullFileName);
00186
00187 if (!in.is_open())
00188 {
00189
00190 std::ostringstream message;
00191 message << "G4DataSet::LoadData - data file " << fullFileName << " not found";
00192
00193 G4Exception("G4CompositeDataSet::LoadData",
00194 "pii00000140",
00195 FatalException,
00196 message.str().c_str());
00197 }
00198
00199 G4DataVector* argEnergies=new G4DataVector;
00200 G4DataVector* argData=new G4DataVector;
00201
00202 G4double a;
00203 bool energyColumn(true);
00204
00205 do
00206 {
00207 in >> a;
00208
00209 if (a!=-1 && a!=-2)
00210 {
00211 if (energyColumn)
00212 {
00213
00214 argEnergies->push_back(a*unitEnergies);
00215 }
00216 else
00217 argData->push_back(a*unitData);
00218 energyColumn=(!energyColumn);
00219 }
00220 }
00221 while (a != -2);
00222
00223 SetEnergiesData(argEnergies, argData, 0);
00224 if (randomSet) BuildPdf();
00225
00226 return true;
00227 }
00228
00229 G4bool G4DataSet::SaveData(const G4String& name) const
00230 {
00231
00232
00233
00234
00235
00236
00237 G4String fullFileName(FullFileName(name));
00238 std::ofstream out(fullFileName);
00239
00240 if (!out.is_open())
00241 {
00242
00243 std::ostringstream message;
00244 message << "G4DataSet:: SaveData - cannot open " << fullFileName;
00245
00246 G4Exception("G4CompositeDataSet::SaveData",
00247 "pii00000150",
00248 FatalException,
00249 message.str().c_str());
00250
00251 }
00252
00253 out.precision(10);
00254 out.width(15);
00255 out.setf(std::ofstream::left);
00256
00257 if (energies!=0 && data!=0)
00258 {
00259 G4DataVector::const_iterator i(energies->begin());
00260 G4DataVector::const_iterator endI(energies->end());
00261 G4DataVector::const_iterator j(data->begin());
00262
00263 while (i!=endI)
00264 {
00265 out.precision(10);
00266 out.width(15);
00267 out.setf(std::ofstream::left);
00268 out << ((*i)/unitEnergies) << ' ';
00269
00270 out.precision(10);
00271 out.width(15);
00272 out.setf(std::ofstream::left);
00273 out << ((*j)/unitData) << std::endl;
00274
00275 i++;
00276 j++;
00277 }
00278 }
00279
00280 out.precision(10);
00281 out.width(15);
00282 out.setf(std::ofstream::left);
00283 out << -1.f << ' ';
00284
00285 out.precision(10);
00286 out.width(15);
00287 out.setf(std::ofstream::left);
00288 out << -1.f << std::endl;
00289
00290 out.precision(10);
00291 out.width(15);
00292 out.setf(std::ofstream::left);
00293 out << -2.f << ' ';
00294
00295 out.precision(10);
00296 out.width(15);
00297 out.setf(std::ofstream::left);
00298 out << -2.f << std::endl;
00299
00300 return true;
00301 }
00302
00303 size_t G4DataSet::FindLowerBound(G4double x) const
00304 {
00305 size_t lowerBound = 0;
00306 size_t upperBound(energies->size() - 1);
00307
00308 while (lowerBound <= upperBound)
00309 {
00310 size_t midBin((lowerBound + upperBound) / 2);
00311
00312 if (x < (*energies)[midBin]) upperBound = midBin - 1;
00313 else lowerBound = midBin + 1;
00314 }
00315
00316 return upperBound;
00317 }
00318
00319
00320 size_t G4DataSet::FindLowerBound(G4double x, G4DataVector* values) const
00321 {
00322 size_t lowerBound = 0;;
00323 size_t upperBound(values->size() - 1);
00324
00325 while (lowerBound <= upperBound)
00326 {
00327 size_t midBin((lowerBound + upperBound) / 2);
00328
00329 if (x < (*values)[midBin]) upperBound = midBin - 1;
00330 else lowerBound = midBin + 1;
00331 }
00332
00333 return upperBound;
00334 }
00335
00336
00337 G4String G4DataSet::FullFileName(const G4String& name) const
00338 {
00339 char* path = getenv("G4PIIDATA");
00340 if (!path)
00341 G4Exception("G4DataSet::FullFileName",
00342 "pii00000160",
00343 FatalException,
00344 "G4PIIDATA environment variable not set");
00345
00346 std::ostringstream fullFileName;
00347 fullFileName << path << '/' << name << z << ".dat";
00348
00349 return G4String(fullFileName.str().c_str());
00350 }
00351
00352
00353 void G4DataSet::BuildPdf()
00354 {
00355 pdf = new G4DataVector;
00356 G4Integrator <G4DataSet, G4double(G4DataSet::*)(G4double)> integrator;
00357
00358 G4int nData = data->size();
00359 pdf->push_back(0.);
00360
00361
00362 G4int i;
00363 G4double totalSum = 0.;
00364 for (i=1; i<nData; i++)
00365 {
00366 G4double xLow = (*energies)[i-1];
00367 G4double xHigh = (*energies)[i];
00368 G4double sum = integrator.Legendre96(this, &G4DataSet::IntegrationFunction, xLow, xHigh);
00369 totalSum = totalSum + sum;
00370 pdf->push_back(totalSum);
00371 }
00372
00373
00374 G4double tot = 0.;
00375 if (totalSum > 0.) tot = 1. / totalSum;
00376 for (i=1; i<nData; i++)
00377 {
00378 (*pdf)[i] = (*pdf)[i] * tot;
00379 }
00380 }
00381
00382
00383 G4double G4DataSet::RandomSelect(G4int ) const
00384 {
00385
00386
00387
00388 if (!pdf) G4Exception("G4DataSet::RandomSelect",
00389 "pii00000170",
00390 FatalException,
00391 "PDF has not been created for this data set");
00392
00393 G4double value = 0.;
00394 G4double x = G4UniformRand();
00395
00396
00397 size_t bin = FindLowerBound(x,pdf);
00398
00399
00400
00401
00402
00403 G4LinInterpolation linearAlgo;
00404 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
00405 else value = algorithm->Calculate(x, bin, *pdf, *energies);
00406
00407
00408 return value;
00409 }
00410
00411 G4double G4DataSet::IntegrationFunction(G4double x)
00412 {
00413
00414
00415 G4double y = 0;
00416
00417
00418 size_t bin = FindLowerBound(x);
00419
00420
00421
00422
00423
00424 G4LinInterpolation linearAlgo;
00425
00426 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
00427 else y = algorithm->Calculate(x, bin, *energies, *data);
00428
00429 return y;
00430 }
00431