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 // 16 Jun 2008 MGP Created; Cross section ma 33 // Documented in: 34 // M.G. Pia et al., PIXE Sim 35 // IEEE Trans. Nucl. Sci., v 36 // 37 // ------------------------------------------- 38 39 #include "G4PixeCrossSectionHandler.hh" 40 #include "G4PhysicalConstants.hh" 41 #include "G4IInterpolator.hh" 42 #include "G4LogLogInterpolator.hh" 43 #include "G4IDataSet.hh" 44 #include "G4DataSet.hh" 45 #include "G4CompositeDataSet.hh" 46 #include "G4PixeShellDataSet.hh" 47 #include "G4ProductionCutsTable.hh" 48 #include "G4Material.hh" 49 #include "G4Element.hh" 50 #include "Randomize.hh" 51 #include "G4SystemOfUnits.hh" 52 #include "G4ParticleDefinition.hh" 53 54 #include <map> 55 #include <vector> 56 #include <fstream> 57 #include <sstream> 58 #include <memory> 59 60 61 G4PixeCrossSectionHandler::G4PixeCrossSectionH 62 { 63 crossSections = 0; 64 interpolation = 0; 65 // Initialise with default values 66 Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV 67 ActiveElements(); 68 } 69 70 71 G4PixeCrossSectionHandler::G4PixeCrossSectionH 72 const G4String& modelK, 73 const G4String& modelL, 74 const G4String& modelM, 75 G4double minE, 76 G4double maxE, 77 G4int bins, 78 G4double unitE, 79 G4double unitData, 80 G4int minZ, 81 G4int maxZ) 82 : interpolation(algorithm), eMin(minE), eMax 83 unit1(unitE), unit2(unitData), zMin(minZ), 84 { 85 crossSections = 0; 86 87 crossModel.push_back(modelK); 88 crossModel.push_back(modelL); 89 crossModel.push_back(modelM); 90 91 //std::cout << "PixeCrossSectionHandler cons 92 // << crossModel[0] 93 // << std::endl; 94 95 ActiveElements(); 96 } 97 98 G4PixeCrossSectionHandler::~G4PixeCrossSection 99 { 100 delete interpolation; 101 interpolation = 0; 102 std::map<G4int,G4IDataSet*,std::less<G4int> 103 104 for (pos = dataMap.begin(); pos != dataMap.e 105 { 106 // The following is a workaround for STL 107 // which does not support the standard a 108 // the syntax pos->second 109 // G4IDataSet* dataSet = pos->second; 110 G4IDataSet* dataSet = (*pos).second; 111 delete dataSet; 112 } 113 114 if (crossSections != 0) 115 { 116 std::size_t n = crossSections->size(); 117 for (std::size_t i=0; i<n; ++i) 118 { 119 delete (*crossSections)[i]; 120 } 121 delete crossSections; 122 crossSections = 0; 123 } 124 } 125 126 void G4PixeCrossSectionHandler::Initialise(G4I 127 const G4String& modelK, 128 const G4String& modelL, 129 const G4String& modelM, 130 G4double minE, G4double maxE, 131 G4int numberOfBins, 132 G4double unitE, G4double unitData 133 G4int minZ, G4int maxZ) 134 { 135 if (algorithm != 0) 136 { 137 delete interpolation; 138 interpolation = algorithm; 139 } 140 else 141 { 142 interpolation = CreateInterpolation(); 143 } 144 145 eMin = minE; 146 eMax = maxE; 147 nBins = numberOfBins; 148 unit1 = unitE; 149 unit2 = unitData; 150 zMin = minZ; 151 zMax = maxZ; 152 153 crossModel.push_back(modelK); 154 crossModel.push_back(modelL); 155 crossModel.push_back(modelM); 156 157 } 158 159 void G4PixeCrossSectionHandler::PrintData() co 160 { 161 std::map<G4int,G4IDataSet*,std::less<G4int> 162 163 for (pos = dataMap.begin(); pos != dataMap.e 164 { 165 // The following is a workaround for STL 166 // which does not support the standard a 167 // the syntax pos->first or pos->second 168 // G4int z = pos->first; 169 // G4IDataSet* dataSet = pos->second; 170 G4int z = (*pos).first; 171 G4IDataSet* dataSet = (*pos).second; 172 G4cout << "---- Data set for Z = " 173 << z 174 << G4endl; 175 dataSet->PrintData(); 176 G4cout << "----------------------------- 177 } 178 } 179 180 void G4PixeCrossSectionHandler::LoadShellData( 181 { 182 std::size_t nZ = activeZ.size(); 183 for (std::size_t i=0; i<nZ; ++i) 184 { 185 G4int Z = (G4int) activeZ[i]; 186 G4IInterpolator* algo = interpolation->C 187 G4IDataSet* dataSet = new G4PixeShellDat 188 189 // Degug printing 190 //std::cout << "PixeCrossSectionHandler: 191 // << Z 192 // << ", modelK = " 193 // << crossModel[0] 194 // << " fileName = " 195 // << fileName 196 // << std::endl; 197 198 dataSet->LoadData(fileName); 199 dataMap[Z] = dataSet; 200 } 201 202 // Build cross sections for materials if not 203 if (! crossSections) 204 { 205 BuildForMaterials(); 206 } 207 208 } 209 210 void G4PixeCrossSectionHandler::Clear() 211 { 212 // Reset the map of data sets: remove the da 213 std::map<G4int,G4IDataSet*,std::less<G4int> 214 215 if(! dataMap.empty()) 216 { 217 for (pos = dataMap.begin(); pos != dataM 218 { 219 // The following is a workaround for STL O 220 // which does not support the standard and 221 // the syntax pos->first or pos->second 222 // G4IDataSet* dataSet = pos->second; 223 G4IDataSet* dataSet = (*pos).second; 224 delete dataSet; 225 dataSet = 0; 226 G4int i = (*pos).first; 227 dataMap[i] = 0; 228 } 229 dataMap.clear(); 230 } 231 232 activeZ.clear(); 233 ActiveElements(); 234 } 235 236 G4double G4PixeCrossSectionHandler::FindValue( 237 { 238 G4double value = 0.; 239 240 std::map<G4int,G4IDataSet*,std::less<G4int> 241 pos = dataMap.find(Z); 242 if (pos!= dataMap.end()) 243 { 244 // The following is a workaround for STL 245 // which does not support the standard a 246 // the syntax pos->first or pos->second 247 // G4IDataSet* dataSet = pos->second; 248 G4IDataSet* dataSet = (*pos).second; 249 value = dataSet->FindValue(energy); 250 } 251 else 252 { 253 G4cout << "WARNING: G4PixeCrossSectionHa 254 << Z << G4endl; 255 } 256 return value; 257 } 258 259 G4double G4PixeCrossSectionHandler::FindValue( 260 G4int shellIndex) const 261 { 262 G4double value = 0.; 263 264 std::map<G4int,G4IDataSet*,std::less<G4int> 265 pos = dataMap.find(Z); 266 if (pos!= dataMap.end()) 267 { 268 // The following is a workaround for STL 269 // which does not support the standard a 270 // the syntax pos->first or pos->second 271 // G4IDataSet* dataSet = pos->second; 272 G4IDataSet* dataSet = (*pos).second; 273 if (shellIndex >= 0) 274 { 275 G4int nComponents = (G4int)dataSet->Number 276 if(shellIndex < nComponents) 277 // The value is the cross section for sh 278 value = dataSet->GetComponent(shellIndex 279 else 280 { 281 G4cout << "WARNING: G4PixeCrossSection 282 << " shellIndex= " << shellIndex 283 << " for Z= " 284 << Z << G4endl; 285 } 286 } else { 287 value = dataSet->FindValue(energy); 288 } 289 } 290 else 291 { 292 G4cout << "WARNING: G4PixeCrossSectionHa 293 << Z << G4endl; 294 } 295 return value; 296 } 297 298 299 G4double G4PixeCrossSectionHandler::ValueForMa 300 G4double energy) const 301 { 302 G4double value = 0.; 303 304 const G4ElementVector* elementVector = mater 305 const G4double* nAtomsPerVolume = material-> 306 std::size_t nElements = material->GetNumberO 307 308 for (std::size_t i=0 ; i<nElements ; ++i) 309 { 310 G4int Z = (G4int) (*elementVector)[i]->G 311 G4double elementValue = FindValue(Z,ener 312 G4double nAtomsVol = nAtomsPerVolume[i]; 313 value += nAtomsVol * elementValue; 314 } 315 316 return value; 317 } 318 319 /* 320 G4IDataSet* G4PixeCrossSectionHandler::Build 321 { 322 // Builds a CompositeDataSet containing the 323 // in the material table 324 325 G4DataVector energyVector; 326 G4double dBin = std::log10(eMax/eMin) / nBin 327 328 for (G4int i=0; i<nBins+1; i++) 329 { 330 energyVector.push_back(std::pow(10., std::lo 331 } 332 333 // Factory method to build cross sections in 334 // related to the type of physics process 335 336 if (crossSections != 0) 337 { // Reset the list of cross sections 338 std::vector<G4IDataSet*>::iterator mat; 339 if (! crossSections->empty()) 340 { 341 for (mat = crossSections->begin(); mat!= cro 342 { 343 G4IDataSet* set = *mat; 344 delete set; 345 set = 0; 346 } 347 crossSections->clear(); 348 delete crossSections; 349 crossSections = 0; 350 } 351 } 352 353 crossSections = BuildCrossSectionsForMateria 354 355 if (crossSections == 0) 356 G4Exception("G4PixeCrossSectionHandler::Buil 357 "pii00000201", 358 FatalException, 359 "crossSections = 0"); 360 361 G4IInterpolator* algo = CreateInterpolation( 362 G4IDataSet* materialSet = new G4CompositeDat 363 364 G4DataVector* energies; 365 G4DataVector* data; 366 367 const G4ProductionCutsTable* theCoupleTable= 368 G4ProductionCutsTable::GetProductionCutsTabl 369 std::size_t numOfCouples = theCoupleTable->G 370 371 372 for (std::size_t m=0; m<numOfCouples; m++) 373 { 374 energies = new G4DataVector; 375 data = new G4DataVector; 376 for (G4int bin=0; bin<nBins; bin++) 377 { 378 G4double energy = energyVector[bin]; 379 energies->push_back(energy); 380 G4IDataSet* matCrossSet = (*crossSections)[m 381 G4double materialCrossSection = 0.0; 382 G4int nElm = matCrossSet->NumberOfComponents 383 for(G4int j=0; j<nElm; j++) { 384 materialCrossSection += matCrossSet->GetComp 385 } 386 387 if (materialCrossSection > 0.) 388 { 389 data->push_back(1./materialCrossSection); 390 } 391 else 392 { 393 data->push_back(DBL_MAX); 394 } 395 } 396 G4IInterpolator* algo = CreateInterpolation( 397 G4IDataSet* dataSet = new G4DataSet(m,energi 398 materialSet->AddComponent(dataSet); 399 } 400 401 return materialSet; 402 } 403 404 */ 405 406 void G4PixeCrossSectionHandler::BuildForMateri 407 { 408 // Builds a CompositeDataSet containing the 409 // in the material table 410 411 G4DataVector energyVector; 412 G4double dBin = std::log10(eMax/eMin) / nBin 413 414 for (G4int i=0; i<nBins+1; i++) 415 { 416 energyVector.push_back(std::pow(10., std 417 } 418 419 if (crossSections != 0) 420 { // Reset the list of cross sections 421 std::vector<G4IDataSet*>::iterator mat; 422 if (! crossSections->empty()) 423 { 424 for (mat = crossSections->begin(); mat!= c 425 { 426 G4IDataSet* set = *mat; 427 delete set; 428 set = 0; 429 } 430 crossSections->clear(); 431 delete crossSections; 432 crossSections = 0; 433 } 434 } 435 436 crossSections = BuildCrossSectionsForMateria 437 438 if (crossSections == 0) 439 G4Exception("G4PixeCrossSectionHandler::Bu 440 "pii00000210", 441 FatalException, 442 ", crossSections = 0"); 443 444 return; 445 } 446 447 448 G4int G4PixeCrossSectionHandler::SelectRandomA 449 G4double e) const 450 { 451 // Select randomly an element within the mat 452 // determined by the cross sections in the d 453 454 G4int nElements = (G4int)material->GetNumber 455 456 // Special case: the material consists of on 457 if (nElements == 1) 458 { 459 G4int Z = (G4int) material->GetZ(); 460 return Z; 461 } 462 463 // Composite material 464 465 const G4ElementVector* elementVector = mater 466 std::size_t materialIndex = material->GetInd 467 468 G4IDataSet* materialSet = (*crossSections)[m 469 G4double materialCrossSection0 = 0.0; 470 G4DataVector cross; 471 cross.clear(); 472 for ( G4int i=0; i < nElements; ++i ) 473 { 474 G4double cr = materialSet->GetComponent( 475 materialCrossSection0 += cr; 476 cross.push_back(materialCrossSection0); 477 } 478 479 G4double random = G4UniformRand() * material 480 481 for (G4int k=0 ; k < nElements ; ++k ) 482 { 483 if (random <= cross[k]) return (G4int) ( 484 } 485 // It should never get here 486 return 0; 487 } 488 489 /* 490 const G4Element* G4PixeCrossSectionHandler:: 491 G4double e) const 492 { 493 // Select randomly an element within the mat 494 // by the cross sections in the data set 495 496 const G4Material* material = couple->GetMate 497 G4Element* nullElement = 0; 498 G4int nElements = material->GetNumberOfEleme 499 const G4ElementVector* elementVector = mater 500 501 // Special case: the material consists of on 502 if (nElements == 1) 503 { 504 G4Element* element = (*elementVector)[0]; 505 return element; 506 } 507 else 508 { 509 // Composite material 510 511 std::size_t materialIndex = couple->GetIndex 512 513 G4IDataSet* materialSet = (*crossSections)[m 514 G4double materialCrossSection0 = 0.0; 515 G4DataVector cross; 516 cross.clear(); 517 for (G4int i=0; i<nElements; i++) 518 { 519 G4double cr = materialSet->GetComponent(i)-> 520 materialCrossSection0 += cr; 521 cross.push_back(materialCrossSection0); 522 } 523 524 G4double random = G4UniformRand() * material 525 526 for (G4int k=0 ; k < nElements ; k++ ) 527 { 528 if (random <= cross[k]) return (*elementVect 529 } 530 // It should never end up here 531 G4cout << "G4PixeCrossSectionHandler::Select 532 return nullElement; 533 } 534 } 535 */ 536 537 538 G4int G4PixeCrossSectionHandler::SelectRandomS 539 { 540 // Select randomly a shell, according to the 541 // in the data set 542 543 // Note for later improvement: it would be u 544 // used shells to improve performance 545 546 G4int shell = 0; 547 548 G4double totCrossSection = FindValue(Z,e); 549 G4double random = G4UniformRand() * totCross 550 G4double partialSum = 0.; 551 552 G4IDataSet* dataSet = nullptr; 553 auto pos = dataMap.find(Z); 554 // The following is a workaround for STL Obj 555 // which does not support the standard and d 556 // the syntax pos->first or pos->second 557 // if (pos != dataMap.end()) dataSet = pos-> 558 if (pos != dataMap.end()) dataSet = (*pos).s 559 560 if (dataSet != nullptr) { 561 G4int nShells = (G4int)dataSet->NumberOfCo 562 for (G4int i=0; i<nShells; ++i) 563 { 564 const G4IDataSet* shellDataSet = dataSet 565 if (shellDataSet != 0) 566 { 567 G4double value = shellDataSet->FindVal 568 partialSum += value; 569 if (random <= partialSum) return i; 570 } 571 } 572 } 573 // It should never get here 574 return shell; 575 } 576 577 void G4PixeCrossSectionHandler::ActiveElements 578 { 579 const G4MaterialTable* materialTable = G4Mat 580 if (materialTable == 0) 581 G4Exception("G4PixeCrossSectionHandler::Ac 582 "pii00000220", 583 FatalException, 584 "no MaterialTable found"); 585 586 std::size_t nMaterials = G4Material::GetNumb 587 588 for (std::size_t mat=0; mat<nMaterials; ++ma 589 { 590 const G4Material* material= (*materialTa 591 const G4ElementVector* elementVector = m 592 const std::size_t nElements = material-> 593 594 for (std::size_t iEl=0; iEl<nElements; + 595 { 596 G4double Z = (*elementVector)[iEl]->GetZ() 597 if (!(activeZ.contains(Z)) && Z >= zMin && 598 { 599 activeZ.push_back(Z); 600 } 601 } 602 } 603 } 604 605 G4IInterpolator* G4PixeCrossSectionHandler::Cr 606 { 607 G4IInterpolator* algorithm = new G4LogLogInt 608 return algorithm; 609 } 610 611 G4int G4PixeCrossSectionHandler::NumberOfCompo 612 { 613 G4int n = 0; 614 615 std::map<G4int,G4IDataSet*,std::less<G4int> 616 pos = dataMap.find(Z); 617 if (pos!= dataMap.end()) 618 { 619 G4IDataSet* dataSet = (*pos).second; 620 n = (G4int)dataSet->NumberOfComponents() 621 } 622 else 623 { 624 G4cout << "WARNING: G4PixeCrossSectionHa 625 << "find Z = " 626 << Z << G4endl; 627 } 628 return n; 629 } 630 631 632 std::vector<G4IDataSet*>* 633 G4PixeCrossSectionHandler::BuildCrossSectionsF 634 { 635 G4DataVector* energies; 636 G4DataVector* data; 637 638 std::vector<G4IDataSet*>* matCrossSections = 639 640 //const G4ProductionCutsTable* theCoupleTabl 641 //std::size_t numOfCouples = theCoupleTable- 642 643 std::size_t nOfBins = energyVector.size(); 644 const auto interpolationAlgo = std::unique_p 645 646 const G4MaterialTable* materialTable = G4Mat 647 if (materialTable == 0) 648 G4Exception("G4PixeCrossSectionHandler::Bu 649 "pii00000230", 650 FatalException, 651 "no MaterialTable found"); 652 653 std::size_t nMaterials = G4Material::GetNumb 654 655 for (std::size_t mat=0; mat<nMaterials; ++ma 656 { 657 const G4Material* material = (*materialT 658 G4int nElements = (G4int)material->GetNu 659 const G4ElementVector* elementVector = m 660 const G4double* nAtomsPerVolume = materi 661 662 G4IInterpolator* algo = interpolationAlg 663 664 G4IDataSet* setForMat = new G4CompositeD 665 666 for (G4int i=0; i<nElements; ++i) { 667 668 G4int Z = (G4int) (*elementVector)[i]- 669 G4double density = nAtomsPerVolume[i]; 670 671 energies = new G4DataVector; 672 data = new G4DataVector; 673 674 675 for (std::size_t bin=0; bin<nOfBins; + 676 { 677 G4double e = energyVector[bin]; 678 energies->push_back(e); 679 G4double cross = 0.; 680 if (Z >= zMin && Z <= zMax) cross = dens 681 data->push_back(cross); 682 } 683 684 G4IInterpolator* algo1 = interpolation 685 G4IDataSet* elSet = new G4DataSet(i,en 686 setForMat->AddComponent(elSet); 687 } 688 689 matCrossSections->push_back(setForMat); 690 } 691 return matCrossSections; 692 } 693 694 695 G4double G4PixeCrossSectionHandler::Microscopi 696 G4double kineticEnergy, 697 G4double Z, 698 G4double deltaCut) const 699 { 700 // Cross section formula is OK for spin=0, 1 701 // Calculates the microscopic cross section 702 // Formula documented in Geant4 Phys. Ref. M 703 // ( it is called for elements, AtomicNumber 704 705 G4double cross = 0.; 706 707 // Particle mass and energy 708 G4double particleMass = particleDef->GetPD 709 G4double energy = kineticEnergy + particle 710 711 // Some kinematics 712 G4double gamma = energy / particleMass; 713 G4double beta2 = 1. - 1. / (gamma * gamma); 714 G4double var = electron_mass_c2 / particleMa 715 G4double tMax = 2. * electron_mass_c2 * (gam 716 717 // Calculate the total cross section 718 719 if ( tMax > deltaCut ) 720 { 721 var = deltaCut / tMax; 722 cross = (1. - var * (1. - beta2 * std::l 723 724 G4double spin = particleDef->GetPDGSpin( 725 726 // +term for spin=1/2 particle 727 if (spin == 0.5) 728 { 729 cross += 0.5 * (tMax - deltaCut) / (energ 730 } 731 // +term for spin=1 particle 732 else if (spin > 0.9 ) 733 { 734 cross += -std::log(var) / (3.*deltaCut) + 735 ((5.+1./var)*0.25 /(energy*energy) - bet 736 } 737 cross *= twopi_mc2_rcl2 * Z / beta2 ; 738 } 739 740 //std::cout << "Microscopic = " << cross/bar 741 // << ", e = " << kineticEnergy/MeV <<std 742 743 return cross; 744 } 745 746