Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch) 29 // 30 // History: 31 // ----------- 32 // 31 Jul 2001 MGP Created, with dummy implementation 33 // 12.09.01 V.Ivanchenko Add param and interpolation of parameters 34 // 04.10.01 V.Ivanchenko Add BindingEnergy method 35 // 25.10.01 MGP Many bug fixes, mostly related to the 36 // management of pointers 37 // 29.11.01 V.Ivanchenko New parametrisation + Excitation 38 // 30.05.02 V.Ivanchenko Format and names of the data files were 39 // chenged to "ion-..." 40 // 17.02.04 V.Ivanchenko Increase buffer size 41 // 42 // ------------------------------------------------------------------- 43 44 #include <fstream> 45 #include <sstream> 46 47 #include "G4RDeIonisationParameters.hh" 48 #include "G4RDVEMDataSet.hh" 49 #include "G4RDShellEMDataSet.hh" 50 #include "G4RDEMDataSet.hh" 51 #include "G4RDCompositeEMDataSet.hh" 52 #include "G4RDLinInterpolation.hh" 53 #include "G4RDLogLogInterpolation.hh" 54 #include "G4RDLinLogLogInterpolation.hh" 55 #include "G4RDSemiLogInterpolation.hh" 56 #include "G4SystemOfUnits.hh" 57 #include "G4Material.hh" 58 #include "G4DataVector.hh" 59 60 G4RDeIonisationParameters:: G4RDeIonisationParameters(G4int minZ, G4int maxZ) 61 : zMin(minZ), zMax(maxZ), 62 length(24) 63 { 64 LoadData(); 65 } 66 67 68 G4RDeIonisationParameters::~G4RDeIonisationParameters() 69 { 70 // Reset the map of data sets: remove the data sets from the map 71 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos; 72 73 for (pos = param.begin(); pos != param.end(); ++pos) 74 { 75 G4RDVEMDataSet* dataSet = (*pos).second; 76 delete dataSet; 77 } 78 79 for (pos = excit.begin(); pos != excit.end(); ++pos) 80 { 81 G4RDVEMDataSet* dataSet = (*pos).second; 82 delete dataSet; 83 } 84 85 activeZ.clear(); 86 } 87 88 89 G4double G4RDeIonisationParameters::Parameter(G4int Z, G4int shellIndex, 90 G4int parameterIndex, 91 G4double e) const 92 { 93 G4double value = 0.; 94 G4int id = Z*100 + parameterIndex; 95 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos; 96 97 pos = param.find(id); 98 if (pos!= param.end()) { 99 G4RDVEMDataSet* dataSet = (*pos).second; 100 G4int nShells = dataSet->NumberOfComponents(); 101 102 if(shellIndex < nShells) { 103 const G4RDVEMDataSet* component = dataSet->GetComponent(shellIndex); 104 const G4DataVector ener = component->GetEnergies(0); 105 G4double ee = std::max(ener.front(),std::min(ener.back(),e)); 106 value = component->FindValue(ee); 107 } else { 108 G4cout << "WARNING: G4IonisationParameters::FindParameter " 109 << "has no parameters for shell= " << shellIndex 110 << "; Z= " << Z 111 << G4endl; 112 } 113 } else { 114 G4cout << "WARNING: G4IonisationParameters::Parameter " 115 << "did not find ID = " 116 << shellIndex << G4endl; 117 } 118 119 return value; 120 } 121 122 G4double G4RDeIonisationParameters::Excitation(G4int Z, G4double e) const 123 { 124 G4double value = 0.; 125 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos; 126 127 pos = excit.find(Z); 128 if (pos!= excit.end()) { 129 G4RDVEMDataSet* dataSet = (*pos).second; 130 131 const G4DataVector ener = dataSet->GetEnergies(0); 132 G4double ee = std::max(ener.front(),std::min(ener.back(),e)); 133 value = dataSet->FindValue(ee); 134 } else { 135 G4cout << "WARNING: G4IonisationParameters::Excitation " 136 << "did not find ID = " 137 << Z << G4endl; 138 } 139 140 return value; 141 } 142 143 144 void G4RDeIonisationParameters::LoadData() 145 { 146 // --------------------------------------- 147 // Please document what are the parameters 148 // --------------------------------------- 149 150 // define active elements 151 152 const G4MaterialTable* materialTable = G4Material::GetMaterialTable(); 153 if (materialTable == 0) 154 G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup", 155 FatalException, "No MaterialTable found!"); 156 157 G4int nMaterials = G4Material::GetNumberOfMaterials(); 158 159 for (G4int m=0; m<nMaterials; m++) { 160 161 const G4Material* material= (*materialTable)[m]; 162 const G4ElementVector* elementVector = material->GetElementVector(); 163 const size_t nElements = material->GetNumberOfElements(); 164 165 for (size_t iEl=0; iEl<nElements; iEl++) { 166 G4Element* element = (*elementVector)[iEl]; 167 G4double Z = element->GetZ(); 168 if (!(activeZ.contains(Z))) { 169 activeZ.push_back(Z); 170 } 171 } 172 } 173 174 const char* path = G4FindDataDir("G4LEDATA"); 175 if (!path) 176 { 177 G4String excep("G4LEDATA environment variable not set!"); 178 G4Exception("G4RDeIonisationParameters::LoadData()", "InvalidSetup", 179 FatalException, excep); 180 } 181 182 G4String pathString(path); 183 G4String path2("/ioni/ion-sp-"); 184 pathString += path2; 185 186 G4double energy, sum; 187 188 size_t nZ = activeZ.size(); 189 190 for (size_t i=0; i<nZ; i++) { 191 192 G4int Z = (G4int)activeZ[i]; 193 std::ostringstream ost; 194 ost << pathString << Z << ".dat"; 195 G4String name(ost.str()); 196 197 std::ifstream file(name); 198 std::filebuf* lsdp = file.rdbuf(); 199 200 if (! (lsdp->is_open()) ) { 201 G4String excep = G4String("Data file: ") 202 + name + G4String(" not found. The version 1.# of G4LEDATA should be used"); 203 G4Exception("G4RDeIonisationParameters::LoadData()", "DataNotFound", 204 FatalException, excep); 205 } 206 207 // The file is organized into: 208 // For each shell there are two lines: 209 // 1st line: 210 // 1st column is the energy of incident e-, 211 // 2d column is the parameter of screan term; 212 // 2d line: 213 // 3 energy (MeV) subdividing different approximation area of the spectrum 214 // 20 point on the spectrum 215 // The shell terminates with the pattern: -1 -1 216 // The file terminates with the pattern: -2 -2 217 218 std::vector<G4RDVEMDataSet*> p; 219 for (size_t k=0; k<length; k++) 220 { 221 G4RDVDataSetAlgorithm* inter = new G4RDLinLogLogInterpolation(); 222 G4RDVEMDataSet* composite = new G4RDCompositeEMDataSet(inter,1.,1.); 223 p.push_back(composite); 224 } 225 226 G4int shell = 0; 227 std::vector<G4DataVector*> a; 228 for (size_t j=0; j<length; j++) 229 { 230 G4DataVector* aa = new G4DataVector(); 231 a.push_back(aa); 232 } 233 G4DataVector e; 234 e.clear(); 235 do { 236 file >> energy >> sum; 237 if (energy == -2) break; 238 239 if (energy > -1) { 240 e.push_back(energy); 241 a[0]->push_back(sum); 242 for (size_t j=0; j<length-1; j++) { 243 G4double qRead; 244 file >> qRead; 245 a[j + 1]->push_back(qRead); 246 } 247 248 } else { 249 250 // End of set for a shell, fill the map 251 for (size_t k=0; k<length; k++) { 252 253 G4RDVDataSetAlgorithm* interp; 254 if(0 == k) interp = new G4RDLinLogLogInterpolation(); 255 else interp = new G4RDLogLogInterpolation(); 256 257 G4DataVector* eVector = new G4DataVector; 258 size_t eSize = e.size(); 259 for (size_t s=0; s<eSize; s++) { 260 eVector->push_back(e[s]); 261 } 262 G4RDVEMDataSet* set = new G4RDEMDataSet(shell,eVector,a[k],interp,1.,1.); 263 264 p[k]->AddComponent(set); 265 } 266 267 // clear vectors 268 for (size_t j2=0; j2<length; j2++) { 269 a[j2] = new G4DataVector(); 270 } 271 shell++; 272 e.clear(); 273 } 274 } while (energy > -2); 275 276 file.close(); 277 278 for (size_t kk=0; kk<length; kk++) 279 { 280 G4int id = Z*100 + kk; 281 param[id] = p[kk]; 282 } 283 } 284 285 G4String pathString_a(path); 286 G4String name_a = pathString_a + G4String("/ioni/ion-ex-av.dat"); 287 std::ifstream file_a(name_a); 288 std::filebuf* lsdp_a = file_a.rdbuf(); 289 G4String pathString_b(path); 290 G4String name_b = pathString_b + G4String("/ioni/ion-ex-sig.dat"); 291 std::ifstream file_b(name_b); 292 std::filebuf* lsdp_b = file_b.rdbuf(); 293 294 if (! (lsdp_a->is_open()) ) { 295 G4String excep = G4String("Cannot open file ") 296 + name_a; 297 G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile", 298 FatalException, excep); 299 } 300 if (! (lsdp_b->is_open()) ) { 301 G4String excep = G4String("Cannot open file ") 302 + name_b; 303 G4Exception("G4RDeIonisationParameters::LoadData()", "CannotOpenFile", 304 FatalException, excep); 305 } 306 307 // The file is organized into two columns: 308 // 1st column is the energy 309 // 2nd column is the corresponding value 310 // The file terminates with the pattern: -1 -1 311 // -2 -2 312 313 G4double ener, ener1, sig, sig1; 314 G4int z = 0; 315 316 G4DataVector e; 317 e.clear(); 318 G4DataVector d; 319 d.clear(); 320 321 do { 322 file_a >> ener >> sig; 323 file_b >> ener1 >> sig1; 324 325 if(ener != ener1) { 326 G4cout << "G4RDeIonisationParameters: problem in excitation data " 327 << "ener= " << ener 328 << " ener1= " << ener1 329 << G4endl; 330 } 331 332 // End of file 333 if (ener == -2) { 334 break; 335 336 // End of next element 337 } else if (ener == -1) { 338 339 z++; 340 G4double Z = (G4double)z; 341 342 // fill map if Z is used 343 if (activeZ.contains(Z)) { 344 345 G4RDVDataSetAlgorithm* inter = new G4RDLinInterpolation(); 346 G4DataVector* eVector = new G4DataVector; 347 G4DataVector* dVector = new G4DataVector; 348 size_t eSize = e.size(); 349 for (size_t s=0; s<eSize; s++) { 350 eVector->push_back(e[s]); 351 dVector->push_back(d[s]); 352 } 353 G4RDVEMDataSet* set = new G4RDEMDataSet(z,eVector,dVector,inter,1.,1.); 354 excit[z] = set; 355 } 356 e.clear(); 357 d.clear(); 358 359 } else { 360 361 e.push_back(ener*MeV); 362 d.push_back(sig1*sig*barn*MeV); 363 } 364 } while (ener != -2); 365 366 file_a.close(); 367 368 } 369 370 371 void G4RDeIonisationParameters::PrintData() const 372 { 373 G4cout << G4endl; 374 G4cout << "===== G4RDeIonisationParameters =====" << G4endl; 375 G4cout << G4endl; 376 377 size_t nZ = activeZ.size(); 378 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos; 379 380 for (size_t i=0; i<nZ; i++) { 381 G4int Z = (G4int)activeZ[i]; 382 383 for (size_t j=0; j<length; j++) { 384 385 G4int index = Z*100 + j; 386 387 pos = param.find(index); 388 if (pos!= param.end()) { 389 G4RDVEMDataSet* dataSet = (*pos).second; 390 size_t nShells = dataSet->NumberOfComponents(); 391 392 for (size_t k=0; k<nShells; k++) { 393 394 G4cout << "===== Z= " << Z << " shell= " << k 395 << " parameter[" << j << "] =====" 396 << G4endl; 397 const G4RDVEMDataSet* comp = dataSet->GetComponent(k); 398 comp->PrintData(); 399 } 400 } 401 } 402 } 403 G4cout << "====================================" << G4endl; 404 } 405 406 407