Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 29 // 30 // History: 31 // ----------- 32 // 1 Aug 2001 MGP Created 33 // 09 Oct 2001 VI Add FindValue with 34 // + NumberOfComponen 35 // 19 Jul 2002 VI Create composite d 36 // 21 Jan 2003 VI Cut per region 37 // 38 // ------------------------------------------- 39 40 #include "G4RDVCrossSectionHandler.hh" 41 #include "G4RDVDataSetAlgorithm.hh" 42 #include "G4RDLogLogInterpolation.hh" 43 #include "G4RDVEMDataSet.hh" 44 #include "G4RDEMDataSet.hh" 45 #include "G4RDCompositeEMDataSet.hh" 46 #include "G4RDShellEMDataSet.hh" 47 #include "G4ProductionCutsTable.hh" 48 #include "G4Material.hh" 49 #include "G4Element.hh" 50 #include "Randomize.hh" 51 #include <map> 52 #include <vector> 53 #include <fstream> 54 #include <sstream> 55 56 57 G4RDVCrossSectionHandler::G4RDVCrossSectionHan 58 { 59 crossSections = 0; 60 interpolation = 0; 61 Initialise(); 62 ActiveElements(); 63 } 64 65 66 G4RDVCrossSectionHandler::G4RDVCrossSectionHan 67 G4double minE, 68 G4double maxE, 69 G4int bins, 70 G4double unitE, 71 G4double unitData, 72 G4int minZ, 73 G4int maxZ) 74 : interpolation(algorithm), eMin(minE), eMax 75 unit1(unitE), unit2(unitData), zMin(minZ), 76 { 77 crossSections = 0; 78 ActiveElements(); 79 } 80 81 G4RDVCrossSectionHandler::~G4RDVCrossSectionHa 82 { 83 delete interpolation; 84 interpolation = 0; 85 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 86 87 for (pos = dataMap.begin(); pos != dataMap.e 88 { 89 // The following is a workaround for STL 90 // which does not support the standard a 91 // the syntax pos->second 92 // G4RDVEMDataSet* dataSet = pos->second 93 G4RDVEMDataSet* dataSet = (*pos).second; 94 delete dataSet; 95 } 96 97 if (crossSections != 0) 98 { 99 size_t n = crossSections->size(); 100 for (size_t i=0; i<n; i++) 101 { 102 delete (*crossSections)[i]; 103 } 104 delete crossSections; 105 crossSections = 0; 106 } 107 } 108 109 void G4RDVCrossSectionHandler::Initialise(G4RD 110 G4double minE, G4double maxE, 111 G4int numberOfBins, 112 G4double unitE, G4double unitData, 113 G4int minZ, G4int maxZ) 114 { 115 if (algorithm != 0) 116 { 117 delete interpolation; 118 interpolation = algorithm; 119 } 120 else 121 { 122 interpolation = CreateInterpolation(); 123 } 124 125 eMin = minE; 126 eMax = maxE; 127 nBins = numberOfBins; 128 unit1 = unitE; 129 unit2 = unitData; 130 zMin = minZ; 131 zMax = maxZ; 132 } 133 134 void G4RDVCrossSectionHandler::PrintData() con 135 { 136 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 137 138 for (pos = dataMap.begin(); pos != dataMap.e 139 { 140 // The following is a workaround for STL 141 // which does not support the standard a 142 // the syntax pos->first or pos->second 143 // G4int z = pos->first; 144 // G4RDVEMDataSet* dataSet = pos->second 145 G4int z = (*pos).first; 146 G4RDVEMDataSet* dataSet = (*pos).second; 147 G4cout << "---- Data set for Z = " 148 << z 149 << G4endl; 150 dataSet->PrintData(); 151 G4cout << "----------------------------- 152 } 153 } 154 155 void G4RDVCrossSectionHandler::LoadData(const 156 { 157 size_t nZ = activeZ.size(); 158 for (size_t i=0; i<nZ; i++) 159 { 160 G4int Z = (G4int) activeZ[i]; 161 162 // Build the complete string identifying 163 164 const char* path = G4FindDataDir("G4LEDA 165 if (!path) 166 { 167 G4String excep = "G4LEDATA environment var 168 G4Exception("G4RDVCrossSectionHandler::Loa 169 "InvalidSetup", FatalExc 170 } 171 172 std::ostringstream ost; 173 ost << path << '/' << fileName << Z << " 174 std::ifstream file(ost.str().c_str()); 175 std::filebuf* lsdp = file.rdbuf(); 176 177 if (! (lsdp->is_open()) ) 178 { 179 G4String excep = "Data file: " + ost.str() 180 G4Exception("G4RDVCrossSectionHandler::Loa 181 "DataNotFound", FatalExc 182 } 183 G4double a = 0; 184 G4int k = 1; 185 G4DataVector* energies = new G4DataVecto 186 G4DataVector* data = new G4DataVector; 187 do 188 { 189 file >> a; 190 G4int nColumns = 2; 191 // The file is organized into two columns: 192 // 1st column is the energy 193 // 2nd column is the corresponding value 194 // The file terminates with the pattern: - 195 // - 196 if (a == -1 || a == -2) 197 { 198 } 199 else 200 { 201 if (k%nColumns != 0) 202 { 203 G4double e = a * unit1; 204 energies->push_back(e); 205 k++; 206 } 207 else if (k%nColumns == 0) 208 { 209 G4double value = a * unit2; 210 data->push_back(value); 211 k = 1; 212 } 213 } 214 } while (a != -2); // end of file 215 216 file.close(); 217 G4RDVDataSetAlgorithm* algo = interpolat 218 G4RDVEMDataSet* dataSet = new G4RDEMData 219 dataMap[Z] = dataSet; 220 } 221 } 222 223 void G4RDVCrossSectionHandler::LoadShellData(c 224 { 225 size_t nZ = activeZ.size(); 226 for (size_t i=0; i<nZ; i++) 227 { 228 G4int Z = (G4int) activeZ[i]; 229 230 // Riccardo Capra <capra@ge.infn.it>: PL 231 // "energies" AND "data" G4DataVector AR 232 // DELETED. WHATSMORE LOADING FILE OPERA 233 // EVEN BEFORE THE CHANGES I DID ON THIS 234 // OPINION SHOULD BE USELESS AND SHOULD 235 236 // Build the complete string identifying 237 238 const char* path = G4FindDataDir("G4LEDA 239 if (!path) 240 { 241 G4String excep = "G4LEDATA environment var 242 G4Exception("G4RDVCrossSectionHandler::Loa 243 "InvalidSetup", FatalExc 244 } 245 246 std::ostringstream ost; 247 248 ost << path << '/' << fileName << Z << " 249 250 std::ifstream file(ost.str().c_str()); 251 std::filebuf* lsdp = file.rdbuf(); 252 253 if (! (lsdp->is_open()) ) 254 { 255 G4String excep = "Data file: " + ost.str() 256 G4Exception("G4RDVCrossSectionHandler::Loa 257 "DataNotFound", FatalExc 258 } 259 G4double a = 0; 260 G4int k = 1; 261 G4DataVector* energies = new G4DataVecto 262 G4DataVector* data = new G4DataVector; 263 do 264 { 265 file >> a; 266 G4int nColumns = 2; 267 // The file is organized into two columns: 268 // 1st column is the energy 269 // 2nd column is the corresponding value 270 // The file terminates with the pattern: - 271 // - 272 if (a == -1 || a == -2) 273 { 274 } 275 else 276 { 277 if (k%nColumns != 0) 278 { 279 G4double e = a * unit1; 280 energies->push_back(e); 281 k++; 282 } 283 else if (k%nColumns == 0) 284 { 285 G4double value = a * unit2; 286 data->push_back(value); 287 k = 1; 288 } 289 } 290 } while (a != -2); // end of file 291 292 file.close(); 293 294 // Riccardo Capra <capra@ge.infn.it>: EN 295 // REMOVED. 296 297 G4RDVDataSetAlgorithm* algo = interpolat 298 G4RDVEMDataSet* dataSet = new G4RDShellE 299 dataSet->LoadData(fileName); 300 dataMap[Z] = dataSet; 301 } 302 } 303 304 void G4RDVCrossSectionHandler::Clear() 305 { 306 // Reset the map of data sets: remove the da 307 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 308 309 if(! dataMap.empty()) 310 { 311 for (pos = dataMap.begin(); pos != dat 312 { 313 // The following is a workaround for STL O 314 // which does not support the standard and 315 // the syntax pos->first or pos->second 316 // G4RDVEMDataSet* dataSet = pos->second; 317 G4RDVEMDataSet* dataSet = (*pos).second; 318 delete dataSet; 319 dataSet = 0; 320 G4int i = (*pos).first; 321 dataMap[i] = 0; 322 } 323 dataMap.clear(); 324 } 325 326 activeZ.clear(); 327 ActiveElements(); 328 } 329 330 G4double G4RDVCrossSectionHandler::FindValue(G 331 { 332 G4double value = 0.; 333 334 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 335 pos = dataMap.find(Z); 336 if (pos!= dataMap.end()) 337 { 338 // The following is a workaround for STL 339 // which does not support the standard a 340 // the syntax pos->first or pos->second 341 // G4RDVEMDataSet* dataSet = pos->second 342 G4RDVEMDataSet* dataSet = (*pos).second; 343 value = dataSet->FindValue(energy); 344 } 345 else 346 { 347 G4cout << "WARNING: G4RDVCrossSectionHan 348 << Z << G4endl; 349 } 350 return value; 351 } 352 353 G4double G4RDVCrossSectionHandler::FindValue(G 354 G4i 355 { 356 G4double value = 0.; 357 358 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 359 pos = dataMap.find(Z); 360 if (pos!= dataMap.end()) 361 { 362 // The following is a workaround for STL 363 // which does not support the standard a 364 // the syntax pos->first or pos->second 365 // G4RDVEMDataSet* dataSet = pos->second 366 G4RDVEMDataSet* dataSet = (*pos).second; 367 if (shellIndex >= 0) 368 { 369 G4int nComponents = dataSet->NumberOfCompo 370 if(shellIndex < nComponents) 371 // - MGP - Why doesn't it use G4RDVEMDat 372 value = dataSet->GetComponent(shellIndex 373 else 374 { 375 G4cout << "WARNING: G4RDVCrossSectionH 376 << " shellIndex= " << shellIndex 377 << " for Z= " 378 << Z << G4endl; 379 } 380 } else { 381 value = dataSet->FindValue(energy); 382 } 383 } 384 else 385 { 386 G4cout << "WARNING: G4RDVCrossSectionHan 387 << Z << G4endl; 388 } 389 return value; 390 } 391 392 393 G4double G4RDVCrossSectionHandler::ValueForMat 394 G4double energy) const 395 { 396 G4double value = 0.; 397 398 const G4ElementVector* elementVector = mater 399 const G4double* nAtomsPerVolume = material-> 400 G4int nElements = material->GetNumberOfEleme 401 402 for (G4int i=0 ; i<nElements ; i++) 403 { 404 G4int Z = (G4int) (*elementVector)[i]->G 405 G4double elementValue = FindValue(Z,ener 406 G4double nAtomsVol = nAtomsPerVolume[i]; 407 value += nAtomsVol * elementValue; 408 } 409 410 return value; 411 } 412 413 414 G4RDVEMDataSet* G4RDVCrossSectionHandler::Buil 415 { 416 // Builds a CompositeDataSet containing the 417 // in the material table 418 419 G4DataVector energyVector; 420 G4double dBin = std::log10(eMax/eMin) / nBin 421 422 for (G4int i=0; i<nBins+1; i++) 423 { 424 energyVector.push_back(std::pow(10., std 425 } 426 427 // Factory method to build cross sections in 428 // related to the type of physics process 429 430 if (crossSections != 0) 431 { // Reset the list of cross sections 432 std::vector<G4RDVEMDataSet*>::iterator m 433 if (! crossSections->empty()) 434 { 435 for (mat = crossSections->begin(); mat!= c 436 { 437 G4RDVEMDataSet* set = *mat; 438 delete set; 439 set = 0; 440 } 441 crossSections->clear(); 442 delete crossSections; 443 crossSections = 0; 444 } 445 } 446 447 crossSections = BuildCrossSectionsForMateria 448 449 if (crossSections == 0) 450 G4Exception("G4RDVCrossSectionHandler::Bui 451 "InvalidCondition", FatalExcep 452 453 G4RDVDataSetAlgorithm* algo = CreateInterpol 454 G4RDVEMDataSet* materialSet = new G4RDCompos 455 456 G4DataVector* energies; 457 G4DataVector* data; 458 459 const G4ProductionCutsTable* theCoupleTable= 460 G4ProductionCutsTable::GetProductionCu 461 size_t numOfCouples = theCoupleTable->GetTab 462 463 464 for (size_t m=0; m<numOfCouples; m++) 465 { 466 energies = new G4DataVector; 467 data = new G4DataVector; 468 for (G4int bin=0; bin<nBins; bin++) 469 { 470 G4double energy = energyVector[bin]; 471 energies->push_back(energy); 472 G4RDVEMDataSet* matCrossSet = (*crossSecti 473 G4double materialCrossSection = 0.0; 474 G4int nElm = matCrossSet->NumberOfCo 475 for(G4int j=0; j<nElm; j++) { 476 materialCrossSection += matCrossSe 477 } 478 479 if (materialCrossSection > 0.) 480 { 481 data->push_back(1./materialCrossSectio 482 } 483 else 484 { 485 data->push_back(DBL_MAX); 486 } 487 } 488 G4RDVDataSetAlgorithm* algo = CreateInte 489 G4RDVEMDataSet* dataSet = new G4RDEMData 490 materialSet->AddComponent(dataSet); 491 } 492 493 return materialSet; 494 } 495 496 G4int G4RDVCrossSectionHandler::SelectRandomAt 497 498 { 499 // Select randomly an element within the mat 500 // determined by the cross sections in the d 501 502 const G4Material* material = couple->GetMate 503 G4int nElements = material->GetNumberOfEleme 504 505 // Special case: the material consists of on 506 if (nElements == 1) 507 { 508 G4int Z = (G4int) material->GetZ(); 509 return Z; 510 } 511 512 // Composite material 513 514 const G4ElementVector* elementVector = mater 515 size_t materialIndex = couple->GetIndex(); 516 517 G4RDVEMDataSet* materialSet = (*crossSection 518 G4double materialCrossSection0 = 0.0; 519 G4DataVector cross; 520 cross.clear(); 521 for ( G4int i=0; i < nElements; i++ ) 522 { 523 G4double cr = materialSet->GetComponent( 524 materialCrossSection0 += cr; 525 cross.push_back(materialCrossSection0); 526 } 527 528 G4double random = G4UniformRand() * material 529 530 for (G4int k=0 ; k < nElements ; k++ ) 531 { 532 if (random <= cross[k]) return (G4int) ( 533 } 534 // It should never get here 535 return 0; 536 } 537 538 const G4Element* G4RDVCrossSectionHandler::Sel 539 G4double e) const 540 { 541 // Select randomly an element within the mat 542 // by the cross sections in the data set 543 544 const G4Material* material = couple->GetMate 545 G4Element* nullElement = 0; 546 G4int nElements = material->GetNumberOfEleme 547 const G4ElementVector* elementVector = mater 548 549 // Special case: the material consists of on 550 if (nElements == 1) 551 { 552 G4Element* element = (*elementVector)[0] 553 return element; 554 } 555 else 556 { 557 // Composite material 558 559 size_t materialIndex = couple->GetIndex( 560 561 G4RDVEMDataSet* materialSet = (*crossSec 562 G4double materialCrossSection0 = 0.0; 563 G4DataVector cross; 564 cross.clear(); 565 for (G4int i=0; i<nElements; i++) 566 { 567 G4double cr = materialSet->GetCompon 568 materialCrossSection0 += cr; 569 cross.push_back(materialCrossSection 570 } 571 572 G4double random = G4UniformRand() * mate 573 574 for (G4int k=0 ; k < nElements ; k++ ) 575 { 576 if (random <= cross[k]) return (*ele 577 } 578 // It should never end up here 579 G4cout << "G4RDVCrossSectionHandler::Sel 580 return nullElement; 581 } 582 } 583 584 G4int G4RDVCrossSectionHandler::SelectRandomSh 585 { 586 // Select randomly a shell, according to the 587 // in the data set 588 589 // Note for later improvement: it would be u 590 // used shells to improve performance 591 592 G4int shell = 0; 593 594 G4double totCrossSection = FindValue(Z,e); 595 G4double random = G4UniformRand() * totCross 596 G4double partialSum = 0.; 597 598 G4RDVEMDataSet* dataSet = 0; 599 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 600 pos = dataMap.find(Z); 601 // The following is a workaround for STL Obj 602 // which does not support the standard and d 603 // the syntax pos->first or pos->second 604 // if (pos != dataMap.end()) dataSet = pos-> 605 if (pos != dataMap.end()) dataSet = (*pos).s 606 607 size_t nShells = dataSet->NumberOfComponents 608 for (size_t i=0; i<nShells; i++) 609 { 610 const G4RDVEMDataSet* shellDataSet = dat 611 if (shellDataSet != 0) 612 { 613 G4double value = shellDataSet->FindValue(e 614 partialSum += value; 615 if (random <= partialSum) return i; 616 } 617 } 618 // It should never get here 619 return shell; 620 } 621 622 void G4RDVCrossSectionHandler::ActiveElements( 623 { 624 const G4MaterialTable* materialTable = G4Mat 625 if (materialTable == 0) 626 G4Exception("G4RDVCrossSectionHandler::Act 627 "InvalidSetup", FatalException 628 629 G4int nMaterials = G4Material::GetNumberOfMa 630 631 for (G4int m=0; m<nMaterials; m++) 632 { 633 const G4Material* material= (*materialTa 634 const G4ElementVector* elementVector = m 635 const G4int nElements = material->GetNum 636 637 for (G4int iEl=0; iEl<nElements; iEl++) 638 { 639 G4Element* element = (*elementVector)[iEl] 640 G4double Z = element->GetZ(); 641 if (!(activeZ.contains(Z)) && Z >= zMin && 642 { 643 activeZ.push_back(Z); 644 } 645 } 646 } 647 } 648 649 G4RDVDataSetAlgorithm* G4RDVCrossSectionHandle 650 { 651 G4RDVDataSetAlgorithm* algorithm = new G4RDL 652 return algorithm; 653 } 654 655 G4int G4RDVCrossSectionHandler::NumberOfCompon 656 { 657 G4int n = 0; 658 659 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 660 pos = dataMap.find(Z); 661 if (pos!= dataMap.end()) 662 { 663 G4RDVEMDataSet* dataSet = (*pos).second; 664 n = dataSet->NumberOfComponents(); 665 } 666 else 667 { 668 G4cout << "WARNING: G4RDVCrossSectionHan 669 << "find Z = " 670 << Z << G4endl; 671 } 672 return n; 673 } 674 675 676