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