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
00054
00055
00056
00057 #include "G4CrossSectionDataSet.hh"
00058 #include "G4VDataSetAlgorithm.hh"
00059 #include "G4EMDataSet.hh"
00060 #include <vector>
00061 #include <fstream>
00062 #include <sstream>
00063
00064
00065 G4CrossSectionDataSet::G4CrossSectionDataSet(G4VDataSetAlgorithm* argAlgorithm,
00066 G4double argUnitEnergies,
00067 G4double argUnitData)
00068 :
00069 algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
00070 {
00071 z = 0;
00072
00073 }
00074
00075
00076 G4CrossSectionDataSet::~G4CrossSectionDataSet()
00077 {
00078 CleanUpComponents();
00079
00080 if (algorithm)
00081 delete algorithm;
00082 }
00083
00084 G4bool G4CrossSectionDataSet::LoadData(const G4String & argFileName)
00085 {
00086 CleanUpComponents();
00087
00088 G4String fullFileName(FullFileName(argFileName));
00089 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
00090
00091 if (!in.is_open())
00092 {
00093 G4String message("data file \"");
00094 message+=fullFileName;
00095 message+="\" not found";
00096 G4Exception("G4CrossSectionDataSet::LoadData",
00097 "em0003",FatalException,message);
00098 return false;
00099 }
00100
00101 std::vector<G4DataVector *> columns;
00102 std::vector<G4DataVector *> log_columns;
00103
00104 std::stringstream *stream(new std::stringstream);
00105 char c;
00106 G4bool comment(false);
00107 G4bool space(true);
00108 G4bool first(true);
00109
00110 try
00111 {
00112 while (!in.eof())
00113 {
00114 in.get(c);
00115
00116 switch (c)
00117 {
00118 case '\r':
00119 case '\n':
00120 if (!first)
00121 {
00122 unsigned long i(0);
00123 G4double value;
00124
00125 while (!stream->eof())
00126 {
00127 (*stream) >> value;
00128
00129 while (i>=columns.size())
00130 {
00131 columns.push_back(new G4DataVector);
00132 log_columns.push_back(new G4DataVector);
00133 }
00134
00135 columns[i]->push_back(value);
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149 if (value <=0.) value = 1e-300;
00150 log_columns[i]->push_back(std::log10(value));
00151
00152 i++;
00153 }
00154
00155 delete stream;
00156 stream=new std::stringstream;
00157 }
00158
00159 first=true;
00160 comment=false;
00161 space=true;
00162 break;
00163
00164 case '#':
00165 comment=true;
00166 break;
00167
00168 case '\t':
00169 c=' ';
00170 case ' ':
00171 if (space)
00172 break;
00173 default:
00174 if (comment)
00175 break;
00176
00177 if (c==' ')
00178 space=true;
00179 else
00180 {
00181 if (space && (!first))
00182 (*stream) << ' ';
00183
00184 first=false;
00185 (*stream) << c;
00186 space=false;
00187 }
00188 }
00189 }
00190 }
00191 catch(const std::ios::failure &e)
00192 {
00193
00194
00195 }
00196
00197 delete stream;
00198
00199 std::vector<G4DataVector *>::size_type maxI(columns.size());
00200
00201 if (maxI<2)
00202 {
00203 G4String message("data file \"");
00204 message+=fullFileName;
00205 message+="\" should have at least two columns";
00206 G4Exception("G4CrossSectionDataSet::LoadData",
00207 "em0005",FatalException,message);
00208 return false;
00209 }
00210
00211 std::vector<G4DataVector*>::size_type i(1);
00212 while (i<maxI)
00213 {
00214 G4DataVector::size_type maxJ(columns[i]->size());
00215
00216 if (maxJ!=columns[0]->size())
00217 {
00218 G4String message("data file \"");
00219 message+=fullFileName;
00220 message+="\" has lines with a different number of columns";
00221 G4Exception("G4CrossSectionDataSet::LoadData",
00222 "em0005",FatalException,message);
00223 return false;
00224 }
00225
00226 G4DataVector::size_type j(0);
00227
00228 G4DataVector *argEnergies=new G4DataVector;
00229 G4DataVector *argData=new G4DataVector;
00230 G4DataVector *argLogEnergies=new G4DataVector;
00231 G4DataVector *argLogData=new G4DataVector;
00232
00233 while(j<maxJ)
00234 {
00235 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
00236 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
00237 argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
00238 argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
00239 j++;
00240 }
00241
00242 AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
00243
00244 i++;
00245 }
00246
00247 i=maxI;
00248 while (i>0)
00249 {
00250 i--;
00251 delete columns[i];
00252 delete log_columns[i];
00253 }
00254
00255 return true;
00256 }
00257
00258
00259 G4bool G4CrossSectionDataSet::LoadNonLogData(const G4String & argFileName)
00260 {
00261 CleanUpComponents();
00262
00263 G4String fullFileName(FullFileName(argFileName));
00264 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
00265
00266 if (!in.is_open())
00267 {
00268 G4String message("data file \"");
00269 message+=fullFileName;
00270 message+="\" not found";
00271 G4Exception("G4CrossSectionDataSet::LoadNonLogData",
00272 "em0003",FatalException,message);
00273 return false;
00274 }
00275
00276 std::vector<G4DataVector *> columns;
00277
00278 std::stringstream *stream(new std::stringstream);
00279 char c;
00280 G4bool comment(false);
00281 G4bool space(true);
00282 G4bool first(true);
00283
00284 try
00285 {
00286 while (!in.eof())
00287 {
00288 in.get(c);
00289
00290 switch (c)
00291 {
00292 case '\r':
00293 case '\n':
00294 if (!first)
00295 {
00296 unsigned long i(0);
00297 G4double value;
00298
00299 while (!stream->eof())
00300 {
00301 (*stream) >> value;
00302
00303 while (i>=columns.size())
00304 {
00305 columns.push_back(new G4DataVector);
00306 }
00307
00308 columns[i]->push_back(value);
00309
00310 i++;
00311 }
00312
00313 delete stream;
00314 stream=new std::stringstream;
00315 }
00316
00317 first=true;
00318 comment=false;
00319 space=true;
00320 break;
00321
00322 case '#':
00323 comment=true;
00324 break;
00325
00326 case '\t':
00327 c=' ';
00328 case ' ':
00329 if (space)
00330 break;
00331 default:
00332 if (comment)
00333 break;
00334
00335 if (c==' ')
00336 space=true;
00337 else
00338 {
00339 if (space && (!first))
00340 (*stream) << ' ';
00341
00342 first=false;
00343 (*stream) << c;
00344 space=false;
00345 }
00346 }
00347 }
00348 }
00349 catch(const std::ios::failure &e)
00350 {
00351
00352
00353 }
00354
00355 delete stream;
00356
00357 std::vector<G4DataVector *>::size_type maxI(columns.size());
00358
00359 if (maxI<2)
00360 {
00361 G4String message("data file \"");
00362 message+=fullFileName;
00363 message+="\" should have at least two columns";
00364 G4Exception("G4CrossSectionDataSet::LoadNonLogData",
00365 "em0005",FatalException,message);
00366 return false;
00367 }
00368
00369 std::vector<G4DataVector*>::size_type i(1);
00370 while (i<maxI)
00371 {
00372 G4DataVector::size_type maxJ(columns[i]->size());
00373
00374 if (maxJ!=columns[0]->size())
00375 {
00376 G4String message("data file \"");
00377 message+=fullFileName;
00378 message+="\" has lines with a different number of columns";
00379 G4Exception("G4CrossSectionDataSet::LoadNonLogData",
00380 "em0005",FatalException,message);
00381 return false;
00382 }
00383
00384 G4DataVector::size_type j(0);
00385
00386 G4DataVector *argEnergies=new G4DataVector;
00387 G4DataVector *argData=new G4DataVector;
00388
00389 while(j<maxJ)
00390 {
00391 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
00392 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
00393 j++;
00394 }
00395
00396 AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
00397
00398 i++;
00399 }
00400
00401 i=maxI;
00402 while (i>0)
00403 {
00404 i--;
00405 delete columns[i];
00406 }
00407
00408 return true;
00409 }
00410
00411
00412 G4bool G4CrossSectionDataSet::SaveData(const G4String & argFileName) const
00413 {
00414 const size_t n(NumberOfComponents());
00415
00416 if (n==0)
00417 {
00418 G4Exception("G4CrossSectionDataSet::SaveData",
00419 "em0005",FatalException,"expected at least one component");
00420 return false;
00421 }
00422
00423 G4String fullFileName(FullFileName(argFileName));
00424 std::ofstream out(fullFileName);
00425
00426 if (!out.is_open())
00427 {
00428 G4String message("cannot open \"");
00429 message+=fullFileName;
00430 message+="\"";
00431 G4Exception("G4CrossSectionDataSet::SaveData",
00432 "em0003",FatalException,message);
00433 return false;
00434 }
00435
00436 G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
00437 G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
00438 G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
00439
00440 size_t k(n);
00441
00442 while (k>0)
00443 {
00444 k--;
00445 iData[k]=GetComponent(k)->GetData(0).begin();
00446 }
00447
00448 while (iEnergies!=iEnergiesEnd)
00449 {
00450 out.precision(10);
00451 out.width(15);
00452 out.setf(std::ofstream::left);
00453 out << ((*iEnergies)/GetUnitEnergies());
00454
00455 k=0;
00456
00457 while (k<n)
00458 {
00459 out << ' ';
00460 out.precision(10);
00461 out.width(15);
00462 out.setf(std::ofstream::left);
00463 out << ((*(iData[k]))/GetUnitData());
00464
00465 iData[k]++;
00466 k++;
00467 }
00468
00469 out << std::endl;
00470
00471 iEnergies++;
00472 }
00473
00474 delete[] iData;
00475
00476 return true;
00477 }
00478
00479
00480 G4String G4CrossSectionDataSet::FullFileName(const G4String& argFileName) const
00481 {
00482 char* path = getenv("G4LEDATA");
00483 if (!path)
00484 {
00485 G4Exception("G4CrossSectionDataSet::FullFileName",
00486 "em0006",FatalException,"G4LEDATA environment variable not set");
00487 return "NULL";
00488 }
00489
00490 std::ostringstream fullFileName;
00491
00492 fullFileName << path << "/" << argFileName << ".dat";
00493
00494 return G4String(fullFileName.str().c_str());
00495 }
00496
00497
00498 G4double G4CrossSectionDataSet::FindValue(G4double argEnergy, G4int ) const
00499 {
00500
00501 G4double value = 0.;
00502
00503 std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
00504 std::vector<G4VEMDataSet *>::const_iterator end(components.end());
00505
00506 while (i!=end)
00507 {
00508 value+=(*i)->FindValue(argEnergy);
00509 i++;
00510 }
00511
00512 return value;
00513 }
00514
00515
00516 void G4CrossSectionDataSet::PrintData(void) const
00517 {
00518 const size_t n(NumberOfComponents());
00519
00520 G4cout << "The data set has " << n << " components" << G4endl;
00521 G4cout << G4endl;
00522
00523 size_t i(0);
00524
00525 while (i<n)
00526 {
00527 G4cout << "--- Component " << i << " ---" << G4endl;
00528 GetComponent(i)->PrintData();
00529 i++;
00530 }
00531 }
00532
00533
00534 void G4CrossSectionDataSet::SetEnergiesData(G4DataVector* argEnergies,
00535 G4DataVector* argData,
00536 G4int argComponentId)
00537 {
00538 G4VEMDataSet * component(components[argComponentId]);
00539
00540 if (component)
00541 {
00542 component->SetEnergiesData(argEnergies, argData, 0);
00543 return;
00544 }
00545
00546 std::ostringstream message;
00547 message << "component " << argComponentId << " not found";
00548
00549 G4Exception("G4CrossSectionDataSet::SetEnergiesData",
00550 "em0005",FatalException,message.str().c_str());
00551 }
00552
00553
00554 void G4CrossSectionDataSet::SetLogEnergiesData(G4DataVector* argEnergies,
00555 G4DataVector* argData,
00556 G4DataVector* argLogEnergies,
00557 G4DataVector* argLogData,
00558 G4int argComponentId)
00559 {
00560 G4VEMDataSet * component(components[argComponentId]);
00561
00562 if (component)
00563 {
00564 component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
00565 return;
00566 }
00567
00568 std::ostringstream message;
00569 message << "component " << argComponentId << " not found";
00570
00571 G4Exception("G4CrossSectionDataSet::SetLogEnergiesData",
00572 "em0005",FatalException,message.str().c_str());
00573 }
00574
00575
00576 void G4CrossSectionDataSet::CleanUpComponents()
00577 {
00578 while (!components.empty())
00579 {
00580 if (components.back()) delete components.back();
00581 components.pop_back();
00582 }
00583 }
00584
00585