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
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 #include "G4VCrossSectionHandler.hh"
00072 #include "G4VDataSetAlgorithm.hh"
00073 #include "G4LogLogInterpolation.hh"
00074 #include "G4VEMDataSet.hh"
00075 #include "G4EMDataSet.hh"
00076 #include "G4CompositeEMDataSet.hh"
00077 #include "G4ShellEMDataSet.hh"
00078 #include "G4ProductionCutsTable.hh"
00079 #include "G4Material.hh"
00080 #include "G4Element.hh"
00081 #include "Randomize.hh"
00082 #include <map>
00083 #include <vector>
00084 #include <fstream>
00085 #include <sstream>
00086
00087
00088 G4VCrossSectionHandler::G4VCrossSectionHandler()
00089 {
00090 crossSections = 0;
00091 interpolation = 0;
00092 Initialise();
00093 ActiveElements();
00094 }
00095
00096
00097 G4VCrossSectionHandler::G4VCrossSectionHandler(G4VDataSetAlgorithm* algorithm,
00098 G4double minE,
00099 G4double maxE,
00100 G4int bins,
00101 G4double unitE,
00102 G4double unitData,
00103 G4int minZ,
00104 G4int maxZ)
00105 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
00106 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
00107 {
00108 crossSections = 0;
00109 ActiveElements();
00110 }
00111
00112 G4VCrossSectionHandler::~G4VCrossSectionHandler()
00113 {
00114 delete interpolation;
00115 interpolation = 0;
00116 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
00117
00118 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00119 {
00120
00121
00122
00123
00124 G4VEMDataSet* dataSet = (*pos).second;
00125 delete dataSet;
00126 }
00127
00128 if (crossSections != 0)
00129 {
00130 size_t n = crossSections->size();
00131 for (size_t i=0; i<n; i++)
00132 {
00133 delete (*crossSections)[i];
00134 }
00135 delete crossSections;
00136 crossSections = 0;
00137 }
00138 }
00139
00140 void G4VCrossSectionHandler::Initialise(G4VDataSetAlgorithm* algorithm,
00141 G4double minE, G4double maxE,
00142 G4int numberOfBins,
00143 G4double unitE, G4double unitData,
00144 G4int minZ, G4int maxZ)
00145 {
00146 if (algorithm != 0)
00147 {
00148 delete interpolation;
00149 interpolation = algorithm;
00150 }
00151 else
00152 {
00153 delete interpolation;
00154 interpolation = CreateInterpolation();
00155 }
00156
00157 eMin = minE;
00158 eMax = maxE;
00159 nBins = numberOfBins;
00160 unit1 = unitE;
00161 unit2 = unitData;
00162 zMin = minZ;
00163 zMax = maxZ;
00164 }
00165
00166 void G4VCrossSectionHandler::PrintData() const
00167 {
00168 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00169
00170 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
00171 {
00172
00173
00174
00175
00176
00177 G4int z = (*pos).first;
00178 G4VEMDataSet* dataSet = (*pos).second;
00179 G4cout << "---- Data set for Z = "
00180 << z
00181 << G4endl;
00182 dataSet->PrintData();
00183 G4cout << "--------------------------------------------------" << G4endl;
00184 }
00185 }
00186
00187 void G4VCrossSectionHandler::LoadData(const G4String& fileName)
00188 {
00189 size_t nZ = activeZ.size();
00190 for (size_t i=0; i<nZ; i++)
00191 {
00192 G4int Z = (G4int) activeZ[i];
00193
00194
00195
00196 char* path = getenv("G4LEDATA");
00197 if (!path)
00198 {
00199 G4Exception("G4VCrossSectionHandler::LoadData",
00200 "em0006",FatalException,"G4LEDATA environment variable not set");
00201 return;
00202 }
00203
00204 std::ostringstream ost;
00205 ost << path << '/' << fileName << Z << ".dat";
00206 std::ifstream file(ost.str().c_str());
00207 std::filebuf* lsdp = file.rdbuf();
00208
00209 if (! (lsdp->is_open()) )
00210 {
00211 G4String excep = "data file: " + ost.str() + " not found";
00212 G4Exception("G4VCrossSectionHandler::LoadData",
00213 "em0003",FatalException,excep);
00214 }
00215 G4double a = 0;
00216 G4int k = 0;
00217 G4int nColumns = 2;
00218
00219 G4DataVector* orig_reg_energies = new G4DataVector;
00220 G4DataVector* orig_reg_data = new G4DataVector;
00221 G4DataVector* log_reg_energies = new G4DataVector;
00222 G4DataVector* log_reg_data = new G4DataVector;
00223
00224 do
00225 {
00226 file >> a;
00227
00228 if (a==0.) a=1e-300;
00229
00230
00231
00232
00233
00234
00235
00236 if (a != -1 && a != -2)
00237 {
00238 if (k%nColumns == 0)
00239 {
00240 orig_reg_energies->push_back(a*unit1);
00241 log_reg_energies->push_back(std::log10(a)+std::log10(unit1));
00242 }
00243 else if (k%nColumns == 1)
00244 {
00245 orig_reg_data->push_back(a*unit2);
00246 log_reg_data->push_back(std::log10(a)+std::log10(unit2));
00247 }
00248 k++;
00249 }
00250 }
00251 while (a != -2);
00252
00253 file.close();
00254 G4VDataSetAlgorithm* algo = interpolation->Clone();
00255
00256 G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,log_reg_energies,log_reg_data,algo);
00257
00258 dataMap[Z] = dataSet;
00259
00260 }
00261 }
00262
00263
00264 void G4VCrossSectionHandler::LoadNonLogData(const G4String& fileName)
00265 {
00266 size_t nZ = activeZ.size();
00267 for (size_t i=0; i<nZ; i++)
00268 {
00269 G4int Z = (G4int) activeZ[i];
00270
00271
00272
00273 char* path = getenv("G4LEDATA");
00274 if (!path)
00275 {
00276 G4Exception("G4VCrossSectionHandler::LoadNonLogData",
00277 "em0006",FatalException,"G4LEDATA environment variable not set");
00278 return;
00279 }
00280
00281 std::ostringstream ost;
00282 ost << path << '/' << fileName << Z << ".dat";
00283 std::ifstream file(ost.str().c_str());
00284 std::filebuf* lsdp = file.rdbuf();
00285
00286 if (! (lsdp->is_open()) )
00287 {
00288 G4String excep = "data file: " + ost.str() + " not found";
00289 G4Exception("G4VCrossSectionHandler::LoadNonLogData",
00290 "em0003",FatalException,excep);
00291 }
00292 G4double a = 0;
00293 G4int k = 0;
00294 G4int nColumns = 2;
00295
00296 G4DataVector* orig_reg_energies = new G4DataVector;
00297 G4DataVector* orig_reg_data = new G4DataVector;
00298
00299 do
00300 {
00301 file >> a;
00302
00303
00304
00305
00306
00307
00308
00309 if (a != -1 && a != -2)
00310 {
00311 if (k%nColumns == 0)
00312 {
00313 orig_reg_energies->push_back(a*unit1);
00314 }
00315 else if (k%nColumns == 1)
00316 {
00317 orig_reg_data->push_back(a*unit2);
00318 }
00319 k++;
00320 }
00321 }
00322 while (a != -2);
00323
00324 file.close();
00325 G4VDataSetAlgorithm* algo = interpolation->Clone();
00326
00327 G4VEMDataSet* dataSet = new G4EMDataSet(Z,orig_reg_energies,orig_reg_data,algo);
00328
00329 dataMap[Z] = dataSet;
00330
00331 }
00332 }
00333
00334 void G4VCrossSectionHandler::LoadShellData(const G4String& fileName)
00335 {
00336 size_t nZ = activeZ.size();
00337 for (size_t i=0; i<nZ; i++)
00338 {
00339 G4int Z = (G4int) activeZ[i];
00340
00341 G4VDataSetAlgorithm* algo = interpolation->Clone();
00342 G4VEMDataSet* dataSet = new G4ShellEMDataSet(Z, algo);
00343
00344 dataSet->LoadData(fileName);
00345
00346 dataMap[Z] = dataSet;
00347 }
00348 }
00349
00350
00351
00352
00353 void G4VCrossSectionHandler::Clear()
00354 {
00355
00356 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
00357
00358 if(! dataMap.empty())
00359 {
00360 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00361 {
00362
00363
00364
00365
00366 G4VEMDataSet* dataSet = (*pos).second;
00367 delete dataSet;
00368 dataSet = 0;
00369 G4int i = (*pos).first;
00370 dataMap[i] = 0;
00371 }
00372 dataMap.clear();
00373 }
00374
00375 activeZ.clear();
00376 ActiveElements();
00377 }
00378
00379 G4double G4VCrossSectionHandler::FindValue(G4int Z, G4double energy) const
00380 {
00381 G4double value = 0.;
00382
00383 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00384 pos = dataMap.find(Z);
00385 if (pos!= dataMap.end())
00386 {
00387
00388
00389
00390
00391 G4VEMDataSet* dataSet = (*pos).second;
00392 value = dataSet->FindValue(energy);
00393 }
00394 else
00395 {
00396 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
00397 << Z << G4endl;
00398 }
00399 return value;
00400 }
00401
00402 G4double G4VCrossSectionHandler::FindValue(G4int Z, G4double energy,
00403 G4int shellIndex) const
00404 {
00405 G4double value = 0.;
00406
00407 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00408 pos = dataMap.find(Z);
00409 if (pos!= dataMap.end())
00410 {
00411
00412
00413
00414
00415 G4VEMDataSet* dataSet = (*pos).second;
00416 if (shellIndex >= 0)
00417 {
00418 G4int nComponents = dataSet->NumberOfComponents();
00419 if(shellIndex < nComponents)
00420
00421 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
00422 else
00423 {
00424 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find"
00425 << " shellIndex= " << shellIndex
00426 << " for Z= "
00427 << Z << G4endl;
00428 }
00429 } else {
00430 value = dataSet->FindValue(energy);
00431 }
00432 }
00433 else
00434 {
00435 G4cout << "WARNING: G4VCrossSectionHandler::FindValue did not find Z = "
00436 << Z << G4endl;
00437 }
00438 return value;
00439 }
00440
00441
00442 G4double G4VCrossSectionHandler::ValueForMaterial(const G4Material* material,
00443 G4double energy) const
00444 {
00445 G4double value = 0.;
00446
00447 const G4ElementVector* elementVector = material->GetElementVector();
00448 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
00449 G4int nElements = material->GetNumberOfElements();
00450
00451 for (G4int i=0 ; i<nElements ; i++)
00452 {
00453 G4int Z = (G4int) (*elementVector)[i]->GetZ();
00454 G4double elementValue = FindValue(Z,energy);
00455 G4double nAtomsVol = nAtomsPerVolume[i];
00456 value += nAtomsVol * elementValue;
00457 }
00458
00459 return value;
00460 }
00461
00462
00463 G4VEMDataSet* G4VCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector* energyCuts)
00464 {
00465
00466
00467
00468 G4DataVector energyVector;
00469 G4double dBin = std::log10(eMax/eMin) / nBins;
00470
00471 for (G4int i=0; i<nBins+1; i++)
00472 {
00473 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
00474 }
00475
00476
00477
00478
00479 if (crossSections != 0)
00480 {
00481 std::vector<G4VEMDataSet*>::iterator mat;
00482 if (! crossSections->empty())
00483 {
00484 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
00485 {
00486 G4VEMDataSet* set = *mat;
00487 delete set;
00488 set = 0;
00489 }
00490 crossSections->clear();
00491 delete crossSections;
00492 crossSections = 0;
00493 }
00494 }
00495
00496 crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
00497
00498 if (crossSections == 0)
00499 {
00500 G4Exception("G4VCrossSectionHandler::BuildMeanFreePathForMaterials",
00501 "em1010",FatalException,"crossSections = 0");
00502 return 0;
00503 }
00504
00505 G4VDataSetAlgorithm* algo = CreateInterpolation();
00506 G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo);
00507
00508
00509 G4DataVector* energies;
00510 G4DataVector* data;
00511 G4DataVector* log_energies;
00512 G4DataVector* log_data;
00513
00514
00515 const G4ProductionCutsTable* theCoupleTable=
00516 G4ProductionCutsTable::GetProductionCutsTable();
00517 size_t numOfCouples = theCoupleTable->GetTableSize();
00518
00519
00520 for (size_t mLocal=0; mLocal<numOfCouples; mLocal++)
00521 {
00522 energies = new G4DataVector;
00523 data = new G4DataVector;
00524 log_energies = new G4DataVector;
00525 log_data = new G4DataVector;
00526 for (G4int bin=0; bin<nBins; bin++)
00527 {
00528 G4double energy = energyVector[bin];
00529 energies->push_back(energy);
00530 log_energies->push_back(std::log10(energy));
00531 G4VEMDataSet* matCrossSet = (*crossSections)[mLocal];
00532 G4double materialCrossSection = 0.0;
00533 G4int nElm = matCrossSet->NumberOfComponents();
00534 for(G4int j=0; j<nElm; j++) {
00535 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
00536 }
00537
00538 if (materialCrossSection > 0.)
00539 {
00540 data->push_back(1./materialCrossSection);
00541 log_data->push_back(std::log10(1./materialCrossSection));
00542 }
00543 else
00544 {
00545 data->push_back(DBL_MAX);
00546 log_data->push_back(std::log10(DBL_MAX));
00547 }
00548 }
00549 G4VDataSetAlgorithm* algoLocal = CreateInterpolation();
00550
00551
00552
00553 G4VEMDataSet* dataSet = new G4EMDataSet(mLocal,energies,data,log_energies,log_data,algoLocal,1.,1.);
00554
00555 materialSet->AddComponent(dataSet);
00556 }
00557
00558 return materialSet;
00559 }
00560
00561
00562 G4int G4VCrossSectionHandler::SelectRandomAtom(const G4MaterialCutsCouple* couple,
00563 G4double e) const
00564 {
00565
00566
00567
00568 const G4Material* material = couple->GetMaterial();
00569 G4int nElements = material->GetNumberOfElements();
00570
00571
00572 if (nElements == 1)
00573 {
00574 G4int Z = (G4int) material->GetZ();
00575 return Z;
00576 }
00577
00578
00579
00580 const G4ElementVector* elementVector = material->GetElementVector();
00581 size_t materialIndex = couple->GetIndex();
00582
00583 G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
00584 G4double materialCrossSection0 = 0.0;
00585 G4DataVector cross;
00586 cross.clear();
00587 for ( G4int i=0; i < nElements; i++ )
00588 {
00589 G4double cr = materialSet->GetComponent(i)->FindValue(e);
00590 materialCrossSection0 += cr;
00591 cross.push_back(materialCrossSection0);
00592 }
00593
00594 G4double random = G4UniformRand() * materialCrossSection0;
00595
00596 for (G4int k=0 ; k < nElements ; k++ )
00597 {
00598 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
00599 }
00600
00601 return 0;
00602 }
00603
00604 const G4Element* G4VCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
00605 G4double e) const
00606 {
00607
00608
00609
00610 const G4Material* material = couple->GetMaterial();
00611 G4Element* nullElement = 0;
00612 G4int nElements = material->GetNumberOfElements();
00613 const G4ElementVector* elementVector = material->GetElementVector();
00614
00615
00616 if (nElements == 1)
00617 {
00618 G4Element* element = (*elementVector)[0];
00619 return element;
00620 }
00621 else
00622 {
00623
00624
00625 size_t materialIndex = couple->GetIndex();
00626
00627 G4VEMDataSet* materialSet = (*crossSections)[materialIndex];
00628 G4double materialCrossSection0 = 0.0;
00629 G4DataVector cross;
00630 cross.clear();
00631 for (G4int i=0; i<nElements; i++)
00632 {
00633 G4double cr = materialSet->GetComponent(i)->FindValue(e);
00634 materialCrossSection0 += cr;
00635 cross.push_back(materialCrossSection0);
00636 }
00637
00638 G4double random = G4UniformRand() * materialCrossSection0;
00639
00640 for (G4int k=0 ; k < nElements ; k++ )
00641 {
00642 if (random <= cross[k]) return (*elementVector)[k];
00643 }
00644
00645 G4cout << "G4VCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
00646 return nullElement;
00647 }
00648 }
00649
00650 G4int G4VCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
00651 {
00652
00653
00654
00655
00656
00657
00658 G4int shell = 0;
00659
00660 G4double totCrossSection = FindValue(Z,e);
00661 G4double random = G4UniformRand() * totCrossSection;
00662 G4double partialSum = 0.;
00663
00664 G4VEMDataSet* dataSet = 0;
00665 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00666 pos = dataMap.find(Z);
00667
00668
00669
00670
00671 if (pos != dataMap.end())
00672 dataSet = (*pos).second;
00673 else
00674 {
00675 G4Exception("G4VCrossSectionHandler::SelectRandomShell",
00676 "em1011",FatalException,"unable to load the dataSet");
00677 return 0;
00678 }
00679
00680 size_t nShells = dataSet->NumberOfComponents();
00681 for (size_t i=0; i<nShells; i++)
00682 {
00683 const G4VEMDataSet* shellDataSet = dataSet->GetComponent(i);
00684 if (shellDataSet != 0)
00685 {
00686 G4double value = shellDataSet->FindValue(e);
00687 partialSum += value;
00688 if (random <= partialSum) return i;
00689 }
00690 }
00691
00692 return shell;
00693 }
00694
00695 void G4VCrossSectionHandler::ActiveElements()
00696 {
00697 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00698 if (materialTable == 0)
00699 G4Exception("G4VCrossSectionHandler::ActiveElements",
00700 "em1001",FatalException,"no MaterialTable found");
00701
00702 G4int nMaterials = G4Material::GetNumberOfMaterials();
00703
00704 for (G4int mLocal2=0; mLocal2<nMaterials; mLocal2++)
00705 {
00706 const G4Material* material= (*materialTable)[mLocal2];
00707 const G4ElementVector* elementVector = material->GetElementVector();
00708 const G4int nElements = material->GetNumberOfElements();
00709
00710 for (G4int iEl=0; iEl<nElements; iEl++)
00711 {
00712 G4Element* element = (*elementVector)[iEl];
00713 G4double Z = element->GetZ();
00714 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
00715 {
00716 activeZ.push_back(Z);
00717 }
00718 }
00719 }
00720 }
00721
00722 G4VDataSetAlgorithm* G4VCrossSectionHandler::CreateInterpolation()
00723 {
00724 G4VDataSetAlgorithm* algorithm = new G4LogLogInterpolation;
00725 return algorithm;
00726 }
00727
00728 G4int G4VCrossSectionHandler::NumberOfComponents(G4int Z) const
00729 {
00730 G4int n = 0;
00731
00732 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00733 pos = dataMap.find(Z);
00734 if (pos!= dataMap.end())
00735 {
00736 G4VEMDataSet* dataSet = (*pos).second;
00737 n = dataSet->NumberOfComponents();
00738 }
00739 else
00740 {
00741 G4cout << "WARNING: G4VCrossSectionHandler::NumberOfComponents did not "
00742 << "find Z = "
00743 << Z << G4endl;
00744 }
00745 return n;
00746 }
00747
00748