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 #include "G4eIonisationParameters.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4VEMDataSet.hh"
00050 #include "G4ShellEMDataSet.hh"
00051 #include "G4EMDataSet.hh"
00052 #include "G4CompositeEMDataSet.hh"
00053 #include "G4LinInterpolation.hh"
00054 #include "G4LogLogInterpolation.hh"
00055 #include "G4LinLogLogInterpolation.hh"
00056 #include "G4SemiLogInterpolation.hh"
00057 #include "G4Material.hh"
00058 #include "G4DataVector.hh"
00059 #include <fstream>
00060 #include <sstream>
00061
00062
00063 G4eIonisationParameters:: G4eIonisationParameters(G4int minZ, G4int maxZ)
00064 : zMin(minZ), zMax(maxZ),
00065 length(24)
00066 {
00067 LoadData();
00068 }
00069
00070
00071 G4eIonisationParameters::~G4eIonisationParameters()
00072 {
00073
00074 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
00075
00076 for (pos = param.begin(); pos != param.end(); ++pos)
00077 {
00078 G4VEMDataSet* dataSet = (*pos).second;
00079 delete dataSet;
00080 }
00081
00082 for (pos = excit.begin(); pos != excit.end(); ++pos)
00083 {
00084 G4VEMDataSet* dataSet = (*pos).second;
00085 delete dataSet;
00086 }
00087
00088 activeZ.clear();
00089 }
00090
00091
00092 G4double G4eIonisationParameters::Parameter(G4int Z, G4int shellIndex,
00093 G4int parameterIndex,
00094 G4double e) const
00095 {
00096 G4double value = 0.;
00097 G4int id = Z*100 + parameterIndex;
00098 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00099
00100 pos = param.find(id);
00101 if (pos!= param.end()) {
00102 G4VEMDataSet* dataSet = (*pos).second;
00103 G4int nShells = dataSet->NumberOfComponents();
00104
00105 if(shellIndex < nShells) {
00106 const G4VEMDataSet* component = dataSet->GetComponent(shellIndex);
00107 const G4DataVector ener = component->GetEnergies(0);
00108 G4double ee = std::max(ener.front(),std::min(ener.back(),e));
00109 value = component->FindValue(ee);
00110 } else {
00111 G4cout << "WARNING: G4IonisationParameters::FindParameter "
00112 << "has no parameters for shell= " << shellIndex
00113 << "; Z= " << Z
00114 << G4endl;
00115 }
00116 } else {
00117 G4cout << "WARNING: G4IonisationParameters::Parameter "
00118 << "did not find ID = "
00119 << shellIndex << G4endl;
00120 }
00121
00122 return value;
00123 }
00124
00125 G4double G4eIonisationParameters::Excitation(G4int Z, G4double e) const
00126 {
00127 G4double value = 0.;
00128 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00129
00130 pos = excit.find(Z);
00131 if (pos!= excit.end()) {
00132 G4VEMDataSet* dataSet = (*pos).second;
00133
00134 const G4DataVector ener = dataSet->GetEnergies(0);
00135 G4double ee = std::max(ener.front(),std::min(ener.back(),e));
00136 value = dataSet->FindValue(ee);
00137 } else {
00138 G4cout << "WARNING: G4IonisationParameters::Excitation "
00139 << "did not find ID = "
00140 << Z << G4endl;
00141 }
00142
00143 return value;
00144 }
00145
00146
00147 void G4eIonisationParameters::LoadData()
00148 {
00149
00150
00151
00152
00153
00154
00155 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00156 if (materialTable == 0)
00157 G4Exception("G4eIonisationParameters::LoadData",
00158 "em1001",FatalException,"Unable to find MaterialTable");
00159
00160 G4int nMaterials = G4Material::GetNumberOfMaterials();
00161
00162 for (G4int mLocal=0; mLocal<nMaterials; mLocal++) {
00163
00164 const G4Material* material= (*materialTable)[mLocal];
00165 const G4ElementVector* elementVector = material->GetElementVector();
00166 const size_t nElements = material->GetNumberOfElements();
00167
00168 for (size_t iEl=0; iEl<nElements; iEl++) {
00169 G4Element* element = (*elementVector)[iEl];
00170 G4double Z = element->GetZ();
00171 if (!(activeZ.contains(Z))) {
00172 activeZ.push_back(Z);
00173 }
00174 }
00175 }
00176
00177 char* path = getenv("G4LEDATA");
00178 if (!path)
00179 {
00180 G4Exception("G4BremsstrahlungParameters::LoadData",
00181 "em0006",FatalException,"G4LEDATA environment variable not set");
00182 return;
00183 }
00184
00185 G4String pathString(path);
00186 G4String path2("/ioni/ion-sp-");
00187 pathString += path2;
00188
00189 G4double energy, sum;
00190
00191 size_t nZ = activeZ.size();
00192
00193 for (size_t i=0; i<nZ; i++) {
00194
00195 G4int Z = (G4int)activeZ[i];
00196 std::ostringstream ost;
00197 ost << pathString << Z << ".dat";
00198 G4String name(ost.str());
00199
00200 std::ifstream file(name);
00201 std::filebuf* lsdp = file.rdbuf();
00202
00203 if (! (lsdp->is_open()) ) {
00204 G4String excep = G4String("data file: ")
00205 + name + G4String(" not found");
00206 G4Exception("G4eIonisationParameters::LoadData",
00207 "em0003",FatalException,excep);
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 std::vector<G4VEMDataSet*> p;
00222 for (size_t k=0; k<length; k++)
00223 {
00224 G4VDataSetAlgorithm* inter = new G4LinLogLogInterpolation();
00225 G4VEMDataSet* composite = new G4CompositeEMDataSet(inter,1.,1.);
00226 p.push_back(composite);
00227 }
00228
00229 G4int shell = 0;
00230 std::vector<G4DataVector*> a;
00231 a.resize(length);
00232 G4DataVector e;
00233 G4bool isReady = false;
00234 do {
00235 file >> energy >> sum;
00236
00237 if (energy < -2)
00238 {
00239 G4String excep("invalid data file");
00240 G4Exception("G4eIonisationParameters::LoadData",
00241 "em0005",FatalException,excep);
00242 return;
00243 }
00244
00245 if (energy == -2) break;
00246
00247 if (energy > -1) {
00248
00249 if(!isReady) {
00250 isReady = true;
00251 e.clear();
00252 for (size_t j=0; j<length; ++j) {
00253 a[j] = new G4DataVector();
00254 }
00255 }
00256 e.push_back(energy);
00257 a[0]->push_back(sum);
00258 for (size_t j=1; j<length; ++j) {
00259 G4double qRead;
00260 file >> qRead;
00261 a[j]->push_back(qRead);
00262 }
00263
00264 } else {
00265
00266
00267 for (size_t k=0; k<length; k++) {
00268
00269 G4VDataSetAlgorithm* interp;
00270 if(0 == k) { interp = new G4LinLogLogInterpolation(); }
00271 else { interp = new G4LogLogInterpolation(); }
00272
00273 G4DataVector* eVector = new G4DataVector;
00274 size_t eSize = e.size();
00275 for (size_t sLocal=0; sLocal<eSize; sLocal++) {
00276 eVector->push_back(e[sLocal]);
00277 }
00278 G4VEMDataSet* set = new G4EMDataSet(shell,eVector,a[k],interp,1.,1.);
00279
00280 p[k]->AddComponent(set);
00281 }
00282
00283
00284 ++shell;
00285 isReady = false;
00286 }
00287 } while (energy > -2);
00288
00289 file.close();
00290
00291 for (size_t kk=0; kk<length; kk++)
00292 {
00293 G4int id = Z*100 + kk;
00294 param[id] = p[kk];
00295 }
00296 }
00297
00298 G4String pathString_a(path);
00299 G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat");
00300 std::ifstream file_a(name_a);
00301 std::filebuf* lsdp_a = file_a.rdbuf();
00302 G4String pathString_b(path);
00303 G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat");
00304 std::ifstream file_b(name_b);
00305 std::filebuf* lsdp_b = file_b.rdbuf();
00306
00307 if (! (lsdp_a->is_open()) ) {
00308 G4String excep = G4String("cannot open file ")
00309 + name_a;
00310 G4Exception("G4eIonisationParameters::LoadData",
00311 "em0003",FatalException,excep);
00312 }
00313 if (! (lsdp_b->is_open()) ) {
00314 G4String excep = G4String("cannot open file ")
00315 + name_b;
00316 G4Exception("G4eIonisationParameters::LoadData",
00317 "em0003",FatalException,excep);
00318 }
00319
00320
00321
00322
00323
00324
00325
00326 G4double ener, ener1, sig, sig1;
00327 G4int z = 0;
00328
00329 G4DataVector e;
00330 e.clear();
00331 G4DataVector d;
00332 d.clear();
00333
00334 do {
00335 file_a >> ener >> sig;
00336 file_b >> ener1 >> sig1;
00337
00338 if(ener != ener1) {
00339 G4cout << "G4eIonisationParameters: problem in excitation data "
00340 << "ener= " << ener
00341 << " ener1= " << ener1
00342 << G4endl;
00343 }
00344
00345
00346 if (ener == -2) {
00347 break;
00348
00349
00350 } else if (ener == -1) {
00351
00352 z++;
00353 G4double Z = (G4double)z;
00354
00355
00356 if (activeZ.contains(Z)) {
00357
00358 G4VDataSetAlgorithm* inter = new G4LinInterpolation();
00359 G4DataVector* eVector = new G4DataVector;
00360 G4DataVector* dVector = new G4DataVector;
00361 size_t eSize = e.size();
00362 for (size_t sLocal2=0; sLocal2<eSize; sLocal2++) {
00363 eVector->push_back(e[sLocal2]);
00364 dVector->push_back(d[sLocal2]);
00365 }
00366 G4VEMDataSet* set = new G4EMDataSet(z,eVector,dVector,inter,1.,1.);
00367 excit[z] = set;
00368 }
00369 e.clear();
00370 d.clear();
00371
00372 } else {
00373
00374 e.push_back(ener*MeV);
00375 d.push_back(sig1*sig*barn*MeV);
00376 }
00377 } while (ener != -2);
00378
00379 file_a.close();
00380
00381 }
00382
00383
00384 void G4eIonisationParameters::PrintData() const
00385 {
00386 G4cout << G4endl;
00387 G4cout << "===== G4eIonisationParameters =====" << G4endl;
00388 G4cout << G4endl;
00389
00390 size_t nZ = activeZ.size();
00391 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00392
00393 for (size_t i=0; i<nZ; i++) {
00394 G4int Z = (G4int)activeZ[i];
00395
00396 for (size_t j=0; j<length; j++) {
00397
00398 G4int index = Z*100 + j;
00399
00400 pos = param.find(index);
00401 if (pos!= param.end()) {
00402 G4VEMDataSet* dataSet = (*pos).second;
00403 size_t nShells = dataSet->NumberOfComponents();
00404
00405 for (size_t k=0; k<nShells; k++) {
00406
00407 G4cout << "===== Z= " << Z << " shell= " << k
00408 << " parameter[" << j << "] ====="
00409 << G4endl;
00410 const G4VEMDataSet* comp = dataSet->GetComponent(k);
00411 comp->PrintData();
00412 }
00413 }
00414 }
00415 }
00416 G4cout << "====================================" << G4endl;
00417 }
00418
00419