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 // Author: Luciano Pandola and Gianfranco Pate 27 // 28 // ------------------------------------------- 29 // History: 30 // 03 Dec 2009 L. Pandola 1st implementati 31 // 25 May 2011 L. Pandola Renamed (make v2 32 // 27 Sep 2013 L. Pandola Migration to MT 33 // 20 Aug 2017 G. Paterno Molecular Interf 34 // 24 Mar 2019 G. Paterno Improved Molecul 35 // 20 Jun 2020 G. Paterno Read qext separa 36 // 27 Aug 2020 G. Paterno Further improvem 37 // 38 // ------------------------------------------- 39 // Class description: 40 // Low Energy Electromagnetic Physics, Rayleig 41 // with the model from Penelope, version 2008 42 // ------------------------------------------- 43 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 46 #include "G4PenelopeRayleighModelMI.hh" 47 48 #include "G4PenelopeSamplingData.hh" 49 #include "G4ParticleDefinition.hh" 50 #include "G4MaterialCutsCouple.hh" 51 #include "G4ProductionCutsTable.hh" 52 #include "G4DynamicParticle.hh" 53 #include "G4ElementTable.hh" 54 #include "G4Element.hh" 55 #include "G4PhysicsFreeVector.hh" 56 #include "G4AutoLock.hh" 57 #include "G4Exp.hh" 58 #include "G4ExtendedMaterial.hh" 59 #include "G4CrystalExtension.hh" 60 #include "G4EmParameters.hh" 61 62 #include "G4PhysicalConstants.hh" 63 #include "G4SystemOfUnits.hh" 64 //....oooOO0OOooo........oooOO0OOooo........oo 65 66 const G4int G4PenelopeRayleighModelMI::fMaxZ; 67 G4PhysicsFreeVector* G4PenelopeRayleighModelMI 68 G4PhysicsFreeVector* G4PenelopeRayleighModelMI 69 70 //....oooOO0OOooo........oooOO0OOooo........oo 71 72 G4PenelopeRayleighModelMI::G4PenelopeRayleighM 73 const G4String& nam) : 74 G4VEmModel(nam), 75 fParticleChange(nullptr),fParticle(nullptr), 76 fSamplingTable(nullptr),fMolInterferenceData 77 fIsInitialised(false),fLocalTable(false),fIs 78 { 79 fIntrinsicLowEnergyLimit = 100.0*eV; 80 fIntrinsicHighEnergyLimit = 100.0*GeV; 81 //SetLowEnergyLimit(fIntrinsicLowEnergyLimit 82 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 83 84 if (part) SetParticle(part); 85 86 fVerboseLevel = 0; 87 // Verbosity scale: 88 // 0 = nothing 89 // 1 = warning for energy non-conservation 90 // 2 = details of energy budget 91 // 3 = calculation of FF and CS, file openin 92 // 4 = entering in methods 93 94 //build the energy grid. It is the same for 95 G4double logenergy = G4Log(fIntrinsicLowEner 96 G4double logmaxenergy = G4Log(1.5*fIntrinsic 97 //finer grid below 160 keV 98 G4double logtransitionenergy = G4Log(160*keV 99 G4double logfactor1 = G4Log(10.)/250.; 100 G4double logfactor2 = logfactor1*10; 101 fLogEnergyGridPMax.push_back(logenergy); 102 do { 103 if (logenergy < logtransitionenergy) 104 logenergy += logfactor1; 105 else 106 logenergy += logfactor2; 107 fLogEnergyGridPMax.push_back(logenergy); 108 } while (logenergy < logmaxenergy); 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oo 112 113 G4PenelopeRayleighModelMI::~G4PenelopeRayleigh 114 { 115 if (IsMaster() || fLocalTable) { 116 117 for(G4int i=0; i<=fMaxZ; ++i) 118 { 119 if(fLogAtomicCrossSection[i]) 120 { 121 delete fLogAtomicCrossSection[i]; 122 fLogAtomicCrossSection[i] = nullptr; 123 } 124 if(fAtomicFormFactor[i]) 125 { 126 delete fAtomicFormFactor[i]; 127 fAtomicFormFactor[i] = nullptr; 128 } 129 } 130 if (fMolInterferenceData) { 131 for (auto& item : (*fMolInterferenceData 132 if (item.second) delete item.second; 133 delete fMolInterferenceData; 134 fMolInterferenceData = nullptr; 135 } 136 if (fKnownMaterials) 137 { 138 delete fKnownMaterials; 139 fKnownMaterials = nullptr; 140 } 141 if (fAngularFunction) 142 { 143 delete fAngularFunction; 144 fAngularFunction = nullptr; 145 } 146 ClearTables(); 147 } 148 } 149 150 //....oooOO0OOooo........oooOO0OOooo........oo 151 152 void G4PenelopeRayleighModelMI::ClearTables() 153 { 154 if (fLogFormFactorTable) { 155 for (auto& item : (*fLogFormFactorTable)) 156 if (item.second) delete item.second; 157 delete fLogFormFactorTable; 158 fLogFormFactorTable = nullptr; //zero expl 159 } 160 161 if (fPMaxTable) { 162 for (auto& item : (*fPMaxTable)) 163 if (item.second) delete item.second; 164 delete fPMaxTable; 165 fPMaxTable = nullptr; //zero explicitly 166 } 167 168 if (fSamplingTable) { 169 for (auto& item : (*fSamplingTable)) 170 if (item.second) delete item.second; 171 delete fSamplingTable; 172 fSamplingTable = nullptr; //zero explicitl 173 } 174 175 return; 176 } 177 178 //....oooOO0OOooo........oooOO0OOooo........oo 179 180 void G4PenelopeRayleighModelMI::Initialise(con 181 const G4DataVector& ) 182 { 183 if (fVerboseLevel > 3) 184 G4cout << "Calling G4PenelopeRayleighModel 185 186 SetParticle(part); 187 188 if (fVerboseLevel) 189 G4cout << "# Molecular Interference is " < 190 191 //Only the master model creates/fills/destro 192 if (IsMaster() && part == fParticle) { 193 //clear tables depending on materials, not 194 ClearTables(); 195 196 //Use here the highest verbosity, from G4E 197 G4int globVerb = G4EmParameters::Instance( 198 if (globVerb > fVerboseLevel) 199 { 200 fVerboseLevel = globVerb; 201 if (fVerboseLevel) 202 G4cout << "Verbosity level of G4PenelopeRa 203 " from G4EmParameters()" << G4endl; 204 } 205 if (fVerboseLevel > 3) 206 G4cout << "Calling G4PenelopeRayleighMod 207 208 //Load the list of known materials and the 209 if (fIsMIActive) 210 { 211 if (!fKnownMaterials) 212 fKnownMaterials = new std::map<G4String,G4 213 if (!fKnownMaterials->size()) 214 LoadKnownMIFFMaterials(); 215 if (!fAngularFunction) 216 { 217 //Create and fill once 218 fAngularFunction = new G4PhysicsFreeVect 219 CalculateThetaAndAngFun(); //angular fun 220 } 221 } 222 if (fIsMIActive && !fMolInterferenceData) 223 fMolInterferenceData = new std::map<G4St 224 if (!fLogFormFactorTable) 225 fLogFormFactorTable = new std::map<const 226 if (!fPMaxTable) 227 fPMaxTable = new std::map<const G4Materi 228 if (!fSamplingTable) 229 fSamplingTable = new std::map<const G4Ma 230 231 //loop on the used materials 232 G4ProductionCutsTable* theCoupleTable = G4 233 234 for (G4int i=0;i<(G4int)theCoupleTable->Ge 235 const G4Material* material = 236 theCoupleTable->GetMaterialCutsCouple(i)->Ge 237 const G4ElementVector* theElementVector 238 239 for (std::size_t j=0;j<material->GetNumb 240 G4int iZ = theElementVector->at(j)->GetZasIn 241 //read data files only in the master 242 if (!fLogAtomicCrossSection[iZ]) 243 ReadDataFile(iZ); 244 } 245 246 //1) Read MI form factors 247 if (fIsMIActive && !fMolInterferenceData 248 ReadMolInterferenceData(material->GetName()) 249 250 //2) If the table has not been built for 251 if (!fLogFormFactorTable->count(material 252 BuildFormFactorTable(material); 253 254 //3) retrieve or build the sampling tabl 255 if (!(fSamplingTable->count(material))) 256 InitializeSamplingAlgorithm(material); 257 258 //4) retrieve or build the pMax data 259 if (!fPMaxTable->count(material)) 260 GetPMaxTable(material); 261 } 262 263 if (fVerboseLevel > 1) { 264 G4cout << G4endl << "Penelope Rayleigh m 265 << "Energy range: " 266 << LowEnergyLimit() / keV << " keV - " 267 << HighEnergyLimit() / GeV << " GeV" 268 << G4endl; 269 } 270 } 271 272 if (fIsInitialised) 273 return; 274 fParticleChange = GetParticleChangeForGamma( 275 fIsInitialised = true; 276 } 277 278 //....oooOO0OOooo........oooOO0OOooo........oo 279 280 void G4PenelopeRayleighModelMI::InitialiseLoca 281 G4VEmModel *masterModel) 282 { 283 if (fVerboseLevel > 3) 284 G4cout << "Calling G4PenelopeRayleighMode 285 286 //Check that particle matches: one might hav 287 //(e.g. for e+ and e-) 288 if (part == fParticle) { 289 290 //Get the const table pointers from the ma 291 const G4PenelopeRayleighModelMI* theModel 292 static_cast<G4PenelopeRayleighModelMI*> 293 294 //Copy pointers to the data tables 295 for(G4int i=0; i<=fMaxZ; ++i) 296 { 297 fLogAtomicCrossSection[i] = theModel->fLogAt 298 fAtomicFormFactor[i] = theModel->fAtomicForm 299 } 300 fMolInterferenceData = theModel->fMolInter 301 fLogFormFactorTable = theModel->fLogFormFa 302 fPMaxTable = theModel->fPMaxTable; 303 fSamplingTable = theModel->fSamplingTable; 304 fKnownMaterials = theModel->fKnownMaterial 305 fAngularFunction = theModel->fAngularFunct 306 307 //Copy the G4DataVector with the grid 308 fLogQSquareGrid = theModel->fLogQSquareGri 309 310 //Same verbosity for all workers, as the m 311 fVerboseLevel = theModel->fVerboseLevel; 312 } 313 return; 314 } 315 316 //....oooOO0OOooo........oooOO0OOooo........oo 317 318 namespace {G4Mutex PenelopeRayleighModelMutex 319 G4double G4PenelopeRayleighModelMI::ComputeCro 320 G4double energy, 321 G4double Z, 322 G4double, 323 G4double, 324 G4double) 325 { 326 //Cross section of Rayleigh scattering in Pe 327 //tabulation, Cuellen et al. (1997), with no 328 //et al. J. Phys. Chem. Ref. Data 4 (1975) 4 329 330 if (fVerboseLevel > 3) 331 G4cout << "Calling CrossSectionPerAtom() o 332 333 G4int iZ = G4int(Z); 334 if (!fLogAtomicCrossSection[iZ]) { 335 //If we are here, it means that Initialize 336 //not filled up. This can happen in a Unit 337 if (fVerboseLevel > 0) { 338 //Issue a G4Exception (warning) only in 339 G4ExceptionDescription ed; 340 ed << "Unable to retrieve the cross sect 341 ed << "This can happen only in Unit Test 342 G4Exception("G4PenelopeRayleighModelMI:: 343 "em2040",JustWarning,ed); 344 } 345 346 //protect file reading via autolock 347 G4AutoLock lock(&PenelopeRayleighModelMute 348 ReadDataFile(iZ); 349 lock.unlock(); 350 } 351 352 G4double cross = 0; 353 G4PhysicsFreeVector* atom = fLogAtomicCrossS 354 if (!atom) { 355 G4ExceptionDescription ed; 356 ed << "Unable to find Z=" << iZ << " in th 357 G4Exception("G4PenelopeRayleighModelMI::Co 358 "em2041",FatalException,ed); 359 return 0; 360 } 361 362 G4double logene = G4Log(energy); 363 G4double logXS = atom->Value(logene); 364 cross = G4Exp(logXS); 365 366 if (fVerboseLevel > 2) { 367 G4cout << "Rayleigh cross section at " << 368 << Z << " = " << cross/barn << " barn" << 369 } 370 return cross; 371 } 372 373 //....oooOO0OOooo........oooOO0OOooo........oo 374 375 void G4PenelopeRayleighModelMI::CalculateTheta 376 { 377 G4double theta = 0; 378 for(G4int k=0; k<fNtheta; k++) { 379 theta += fDTheta; 380 G4double value = (1+std::cos(theta)*std::c 381 fAngularFunction->PutValue(k,theta,value); 382 if (fVerboseLevel > 3) 383 G4cout << "theta[" << k << "]: " << fAn 384 << ", angFun[" << k << "]: " << (*fAngu 385 } 386 } 387 388 //....oooOO0OOooo........oooOO0OOooo........oo 389 390 G4double G4PenelopeRayleighModelMI::CalculateQ 391 { 392 G4double lambda,x,q,q2 = 0; 393 394 lambda = hbarc*twopi/energy; 395 x = 1./lambda*std::sin(angle/2.); 396 q = 2.*h_Planck*x/(electron_mass_c2/c_light) 397 398 if (fVerboseLevel > 3) { 399 G4cout << "E: " << energy/keV << " keV, la 400 << ", x: " << x*nm << ", q: " << q << G4e 401 } 402 q2 = std::pow(q,2); 403 return q2; 404 } 405 406 //....oooOO0OOooo........oooOO0OOooo........oo 407 408 //Overriding of parent's (G4VEmModel) method 409 G4double G4PenelopeRayleighModelMI::CrossSecti 410 const G4ParticleDefinition* p, 411 G4double energy, 412 G4double, 413 G4double) 414 { 415 //check if we are in a Unit Test (only for t 416 static G4bool amInAUnitTest = false; 417 if (G4ProductionCutsTable::GetProductionCuts 418 { 419 amInAUnitTest = true; 420 G4ExceptionDescription ed; 421 ed << "The ProductionCuts table is empty 422 ed << "This should happen only in Unit T 423 G4Exception("G4PenelopeRayleighModelMI:: 424 "em2019",JustWarning,ed); 425 } 426 //If the material does not have a MIFF, cont 427 G4String matname = material->GetName(); 428 if (amInAUnitTest) 429 { 430 //No need to lock, as this is always cal 431 const G4ElementVector* theElementVector 432 //protect file reading via autolock 433 for (std::size_t j=0;j<material->GetNumb 434 G4int iZ = theElementVector->at(j)->GetZasIn 435 if (!fLogAtomicCrossSection[iZ]) { 436 ReadDataFile(iZ); 437 } 438 } 439 if (fIsMIActive) 440 ReadMolInterferenceData(matname); 441 if (!fLogFormFactorTable->count(material 442 BuildFormFactorTable(material); 443 if (!(fSamplingTable->count(material))) 444 InitializeSamplingAlgorithm(material); 445 if (!fPMaxTable->count(material)) 446 GetPMaxTable(material); 447 } 448 G4bool useMIFF = fIsMIActive && (fMolInterfe 449 if (!useMIFF) 450 { 451 if (fVerboseLevel > 2) 452 G4cout << "Rayleigh CS of: " << matname << " 453 return G4VEmModel::CrossSectionPerVolume 454 } 455 456 // If we are here, it means that we have to 457 if (fVerboseLevel > 2) 458 G4cout << "Rayleigh CS of: " << matname 459 << " calculated through integration of th 460 461 G4double cs = 0; 462 463 //force null cross-section if below the low- 464 if (energy < LowEnergyLimit()) 465 return cs; 466 467 //if the material is a CRYSTAL, forbid this 468 if (material->IsExtended() && material->GetN 469 G4ExtendedMaterial* extendedmaterial = (G4 470 G4CrystalExtension* crystalExtension = (G4 471 if (crystalExtension != 0) { 472 G4cout << "The material has a crystallin 473 return 0; 474 } 475 } 476 477 //GET MATERIAL INFORMATION 478 G4double atomDensity = material->GetTotNbOfA 479 std::size_t nElements = material->GetNumberO 480 const G4ElementVector* elementVector = mater 481 const G4double* fractionVector = material->G 482 483 //Stoichiometric factors 484 std::vector<G4double> *StoichiometricFactors 485 for (std::size_t i=0;i<nElements;++i) { 486 G4double fraction = fractionVector[i]; 487 G4double atomicWeigth = (*elementVector)[i 488 StoichiometricFactors->push_back(fraction/ 489 } 490 G4double MaxStoichiometricFactor = 0.; 491 for (std::size_t i=0;i<nElements;++i) { 492 if ((*StoichiometricFactors)[i] > MaxStoic 493 MaxStoichiometricFactor = (*Stoichiometr 494 } 495 for (std::size_t i=0;i<nElements;++i) { 496 (*StoichiometricFactors)[i] /= MaxStoichi 497 } 498 499 //Equivalent atoms per molecule 500 G4double atPerMol = 0.; 501 for (std::size_t i=0;i<nElements;++i) 502 atPerMol += (*StoichiometricFactors)[i]; 503 G4double moleculeDensity = 0.; 504 if (atPerMol != 0.) moleculeDensity = atomDe 505 506 if (fVerboseLevel > 2) 507 G4cout << "Material " << material->GetName 508 << "per molecule and " << moleculeDensity 509 510 //Equivalent molecular weight (dimensionless 511 G4double MolWeight = 0.; 512 for (std::size_t i=0;i<nElements;++i) 513 MolWeight += (*StoichiometricFactors)[i]*( 514 515 if (fVerboseLevel > 2) 516 G4cout << "Molecular weight of " << matnam 517 518 G4double IntegrandFun[fNtheta]; 519 for (G4int k=0; k<fNtheta; k++) { 520 G4double theta = fAngularFunction->Energy( 521 G4double F2 = GetFSquared(material,Calcula 522 IntegrandFun[k] = (*fAngularFunction)[k]*F 523 } 524 525 G4double constant = pi*classic_electr_radius 526 cs = constant*IntegrateFun(IntegrandFun,fNth 527 528 //Now cs is the cross section per molecule, 529 G4double csvolume = cs*moleculeDensity; 530 531 //print CS and mfp 532 if (fVerboseLevel > 2) 533 G4cout << "Rayleigh CS of " << matname << 534 << " keV: " << cs/barn << " barn" 535 << ", mean free path: " << 1./csvolume/mm 536 537 delete StoichiometricFactors; 538 //return CS 539 return csvolume; 540 } 541 542 //....oooOO0OOooo........oooOO0OOooo........oo 543 544 void G4PenelopeRayleighModelMI::BuildFormFacto 545 { 546 if (fVerboseLevel > 3) 547 G4cout << "Calling BuildFormFactorTable() 548 549 //GET MATERIAL INFORMATION 550 std::size_t nElements = material->GetNumberO 551 const G4ElementVector* elementVector = mater 552 const G4double* fractionVector = material->G 553 554 //Stoichiometric factors 555 std::vector<G4double> *StoichiometricFactors 556 for (std::size_t i=0;i<nElements;++i) { 557 G4double fraction = fractionVector[i]; 558 G4double atomicWeigth = (*elementVector)[i 559 StoichiometricFactors->push_back(fraction/ 560 } 561 G4double MaxStoichiometricFactor = 0.; 562 for (std::size_t i=0;i<nElements;++i) { 563 if ((*StoichiometricFactors)[i] > MaxStoic 564 MaxStoichiometricFactor = (*Stoichiometr 565 } 566 if (MaxStoichiometricFactor<1e-16) { 567 G4ExceptionDescription ed; 568 ed << "Inconsistent data of atomic composi 569 G4Exception("G4PenelopeRayleighModelMI::Bu 570 "em2042",FatalException,ed); 571 } 572 for (std::size_t i=0;i<nElements;++i) 573 (*StoichiometricFactors)[i] /= MaxStoichi 574 575 //Equivalent molecular weight (dimensionless 576 G4double MolWeight = 0.; 577 for (std::size_t i=0;i<nElements;++i) 578 MolWeight += (*StoichiometricFactors)[i]*( 579 580 //CREATE THE FORM FACTOR TABLE 581 //First, the form factors are retrieved [F/s 582 //Then, they are squared and multiplied for 583 //This makes difference for CS calculation, 584 G4PhysicsFreeVector* theFFVec = new G4Physic 585 /*spline=*/true); 586 587 G4String matname = material->GetName(); 588 G4String aFileNameFF = ""; 589 590 //retrieve MIdata (fFileNameFF) 591 G4MIData* dataMI = GetMIData(material); 592 if (dataMI) 593 aFileNameFF = dataMI->GetFilenameFF(); 594 595 //read the MIFF from a file passed by the us 596 if (fIsMIActive && aFileNameFF != "") { 597 if (fVerboseLevel) 598 G4cout << "Read MIFF for " << matname << 599 600 ReadMolInterferenceData(matname,aFileNameF 601 G4PhysicsFreeVector* Fvector = fMolInterf 602 603 for (std::size_t k=0;k<fLogQSquareGrid.siz 604 G4double q = std::pow(G4Exp(fLogQSquareG 605 G4double f = Fvector->Value(q); 606 G4double ff2 = f*f*MolWeight; 607 if (ff2) 608 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo 609 } 610 } 611 //retrieve the MIFF from the database or use 612 else { 613 //medical material: composition of fat, wa 614 if (fIsMIActive && matname.find("MedMat") 615 if (fVerboseLevel) 616 G4cout << "Build MIFF from components for " 617 618 //get the material composition from its 619 G4int ki, kf=6, ktot=19; 620 G4double comp[4]; 621 G4String compstring = matname.substr(kf+ 622 for (std::size_t j=0; j<4; ++j) { 623 ki = kf+1; 624 kf = ki+4; 625 compstring = matname.substr(ki, 4); 626 comp[j] = atof(compstring.c_str()); 627 if (fVerboseLevel > 2) 628 G4cout << " -- MedMat comp[" << j+1 << "]: 629 } 630 631 const char* path = G4FindDataDir("G4LEDA 632 if (!path) { 633 G4String excep = "G4LEDATA environment varia 634 G4Exception("G4PenelopeRayleighModelMI::Buil 635 "em0006",FatalException,excep); 636 } 637 638 if (!fMolInterferenceData->count("Fat_MI 639 ReadMolInterferenceData("Fat_MI"); 640 G4PhysicsFreeVector* fatFF = fMolInterfe 641 642 if (!fMolInterferenceData->count("Water_ 643 ReadMolInterferenceData("Water_MI"); 644 G4PhysicsFreeVector* waterFF = fMolInter 645 646 if (!fMolInterferenceData->count("BoneMa 647 ReadMolInterferenceData("BoneMatrix_MI"); 648 G4PhysicsFreeVector* bonematrixFF = fMol 649 650 if (!fMolInterferenceData->count("Minera 651 ReadMolInterferenceData("Mineral_MI"); 652 G4PhysicsFreeVector* mineralFF = fMolInt 653 654 //get and combine the molecular form fac 655 for (std::size_t k=0;k<fLogQSquareGrid.s 656 G4double ff2 = 0; 657 G4double q = std::pow(G4Exp(fLogQSquareGrid[ 658 G4double ffat = fatFF->Value(q); 659 G4double fwater = waterFF->Value(q); 660 G4double fbonematrix = bonematrixFF->Value(q 661 G4double fmineral = mineralFF->Value(q); 662 ff2 = comp[0]*ffat*ffat+comp[1]*fwater*fwate 663 comp[2]*fbonematrix*fbonematrix+comp[3]*fm 664 ff2 *= MolWeight; 665 if (ff2) theFFVec->PutValue(k,fLogQSquareGri 666 } 667 } 668 //other materials with MIFF (from the data 669 else if (fIsMIActive && fMolInterferenceDa 670 if (fVerboseLevel) 671 G4cout << "Read MIFF from database " << matn 672 G4PhysicsFreeVector* FF = fMolInterferen 673 for (std::size_t k=0;k<fLogQSquareGrid.s 674 G4double ff2 = 0; 675 G4double q = std::pow(G4Exp(fLogQSquareGrid[ 676 G4double f = FF->Value(q); 677 ff2 = f*f*MolWeight; 678 if (ff2) theFFVec->PutValue(k,fLogQSquareGri 679 } 680 } 681 //IAM 682 else { 683 if (fVerboseLevel) 684 G4cout << "FF of " << matname << " calculate 685 for (std::size_t k=0;k<fLogQSquareGrid.s 686 G4double ff2 = 0; 687 for (std::size_t i=0;i<nElements;++i) { 688 G4int iZ = (*elementVector)[i]->GetZasInt( 689 G4PhysicsFreeVector* theAtomVec = fAtomicF 690 G4double q = std::pow(G4Exp(fLogQSquareGri 691 G4double f = theAtomVec->Value(q); 692 ff2 += f*f*(*StoichiometricFactors)[i]; 693 } 694 if (ff2) theFFVec->PutValue(k,fLogQSquareGri 695 } 696 } 697 } 698 theFFVec->FillSecondDerivatives(); 699 fLogFormFactorTable->insert(std::make_pair(m 700 701 if (fVerboseLevel > 3) 702 DumpFormFactorTable(material); 703 delete StoichiometricFactors; 704 705 return; 706 } 707 708 //....oooOO0OOooo........oooOO0OOooo........oo 709 710 void G4PenelopeRayleighModelMI::SampleSecondar 711 const G4MaterialCutsCouple* coup 712 const G4DynamicParticle* aDynami 713 G4double, 714 G4double) 715 { 716 // Sampling of the Rayleigh final state (nam 717 // from the Penelope2008 model. The scatteri 718 // cross section dOmega/d(cosTheta) from Bor 719 // anomalous scattering effects. The Form Fa 720 // analytical cross section is retrieved via 721 // are tabulated for F(Q). Form factor for c 722 // the additivity rule. The sampling from th 723 // Transform with Aliasing (RITA) algorithm; 724 // for each material and managed by G4Penelo 725 // The sampling algorithm (rejection method) 726 // increases with energy. For E=100 keV the 727 // hydrogen and uranium, respectively. 728 729 if (fVerboseLevel > 3) 730 G4cout << "Calling SamplingSecondaries() o 731 732 G4double photonEnergy0 = aDynamicGamma->GetK 733 734 if (photonEnergy0 <= fIntrinsicLowEnergyLimi 735 fParticleChange->ProposeTrackStatus(fStopA 736 fParticleChange->SetProposedKineticEnergy( 737 fParticleChange->ProposeLocalEnergyDeposit 738 return; 739 } 740 741 G4ParticleMomentum photonDirection0 = aDynam 742 const G4Material* theMat = couple->GetMateri 743 744 //1) Verify if tables are ready 745 //Either Initialize() was not called, or we 746 //not invoked 747 if (!fPMaxTable || !fSamplingTable || !fLogF 748 //create a **thread-local** version of the 749 //Unit Tests 750 fLocalTable = true; 751 if (!fLogFormFactorTable) 752 fLogFormFactorTable = new std::map<const 753 if (!fPMaxTable) 754 fPMaxTable = new std::map<const G4Materi 755 if (!fSamplingTable) 756 fSamplingTable = new std::map<const G4Ma 757 if (fIsMIActive && !fMolInterferenceData) 758 fMolInterferenceData = new std::map<G4St 759 } 760 761 if (!fSamplingTable->count(theMat)) { 762 //If we are here, it means that Initialize 763 //not filled up. This can happen in a Unit 764 if (fVerboseLevel > 0) { 765 //Issue a G4Exception (warning) only in 766 G4ExceptionDescription ed; 767 ed << "Unable to find the fSamplingTable 768 theMat->GetName() << G4endl; 769 ed << "This can happen only in Unit Test 770 G4Exception("G4PenelopeRayleighModelMI:: 771 "em2019",JustWarning,ed); 772 } 773 const G4ElementVector* theElementVector = 774 //protect file reading via autolock 775 G4AutoLock lock(&PenelopeRayleighModelMute 776 for (std::size_t j=0;j<theMat->GetNumberOf 777 G4int iZ = theElementVector->at(j)->GetZ 778 if (!fLogAtomicCrossSection[iZ]) { 779 lock.lock(); 780 ReadDataFile(iZ); 781 lock.unlock(); 782 } 783 } 784 lock.lock(); 785 786 //1) If the table has not been built for t 787 if (!fLogFormFactorTable->count(theMat)) 788 BuildFormFactorTable(theMat); 789 790 //2) retrieve or build the sampling table 791 if (!(fSamplingTable->count(theMat))) 792 InitializeSamplingAlgorithm(theMat); 793 794 //3) retrieve or build the pMax data 795 if (!fPMaxTable->count(theMat)) 796 GetPMaxTable(theMat); 797 798 lock.unlock(); 799 } 800 801 //Ok, restart the job 802 G4PenelopeSamplingData* theDataTable = fSamp 803 G4PhysicsFreeVector* thePMax = fPMaxTable->f 804 G4double cosTheta = 1.0; 805 806 //OK, ready to go! 807 G4double qmax = 2.0*photonEnergy0/electron_m 808 809 if (qmax < 1e-10) { //very low momentum tran 810 G4bool loopAgain=false; 811 do { 812 loopAgain = false; 813 cosTheta = 1.0-2.0*G4UniformRand(); 814 G4double G = 0.5*(1+cosTheta*cosTheta); 815 if (G4UniformRand()>G) 816 loopAgain = true; 817 } while(loopAgain); 818 } 819 else { //larger momentum transfer 820 std::size_t nData = theDataTable->GetNumbe 821 G4double LastQ2inTheTable = theDataTable-> 822 G4double q2max = std::min(qmax*qmax,LastQ2 823 824 G4bool loopAgain = false; 825 G4double MaxPValue = thePMax->Value(photon 826 G4double xx=0; 827 828 //Sampling by rejection method. The reject 829 //G = 0.5*(1+cos^2(theta)) 830 do { 831 loopAgain = false; 832 G4double RandomMax = G4UniformRand()*Max 833 xx = theDataTable->SampleValue(RandomMax 834 835 //xx is a random value of q^2 in (0,q2ma 836 //F(Q^2) via the RITA algorithm 837 if (xx > q2max) 838 loopAgain = true; 839 cosTheta = 1.0-2.0*xx/q2max; 840 G4double G = 0.5*(1+cosTheta*cosTheta); 841 if (G4UniformRand()>G) 842 loopAgain = true; 843 } while(loopAgain); 844 } 845 846 G4double sinTheta = std::sqrt(1-cosTheta*cos 847 848 //Scattered photon angles. ( Z - axis along 849 G4double phi = twopi * G4UniformRand() ; 850 G4double dirX = sinTheta*std::cos(phi); 851 G4double dirY = sinTheta*std::sin(phi); 852 G4double dirZ = cosTheta; 853 854 //Update G4VParticleChange for the scattered 855 G4ThreeVector photonDirection1(dirX, dirY, d 856 857 photonDirection1.rotateUz(photonDirection0); 858 fParticleChange->ProposeMomentumDirection(ph 859 fParticleChange->SetProposedKineticEnergy(ph 860 861 return; 862 } 863 864 //....oooOO0OOooo........oooOO0OOooo........oo 865 866 void G4PenelopeRayleighModelMI::ReadDataFile(c 867 { 868 if (fVerboseLevel > 2) { 869 G4cout << "G4PenelopeRayleighModelMI::Read 870 G4cout << "Going to read Rayleigh data fil 871 } 872 873 const char* path = G4FindDataDir("G4LEDATA") 874 if (!path) { 875 G4String excep = "G4LEDATA environment var 876 G4Exception("G4PenelopeRayleighModelMI::Re 877 "em0006",FatalException,excep); 878 return; 879 } 880 881 //Read first the cross section file (all the 882 std::ostringstream ost; 883 if (Z>9) 884 ost << path << "/penelope/rayleigh/pdgra" 885 else 886 ost << path << "/penelope/rayleigh/pdgra0" 887 std::ifstream file(ost.str().c_str()); 888 889 if (!file.is_open()) { 890 G4String excep = "Data file " + G4String(o 891 G4Exception("G4PenelopeRayleighModelMI::Re 892 "em0003",FatalException,excep); 893 } 894 895 G4int readZ = 0; 896 std::size_t nPoints = 0; 897 file >> readZ >> nPoints; 898 899 if (readZ != Z || nPoints <= 0 || nPoints >= 900 G4ExceptionDescription ed; 901 ed << "Corrupted data file for Z=" << Z << 902 G4Exception("G4PenelopeRayleighModelMI::Re 903 "em0005",FatalException,ed); 904 return; 905 } 906 907 fLogAtomicCrossSection[Z] = new G4PhysicsFre 908 G4double ene=0,f1=0,f2=0,xs=0; 909 for (std::size_t i=0;i<nPoints;++i) { 910 file >> ene >> f1 >> f2 >> xs; 911 //dimensional quantities 912 ene *= eV; 913 xs *= cm2; 914 fLogAtomicCrossSection[Z]->PutValue(i,G4Lo 915 if (file.eof() && i != (nPoints-1)) { //fi 916 G4ExceptionDescription ed ; 917 ed << "Corrupted data file for Z=" << Z 918 ed << "Found less than " << nPoints << " 919 G4Exception("G4PenelopeRayleighModelMI:: 920 "em0005",FatalException,ed); 921 } 922 } 923 file.close(); 924 925 //Then, read the extended momentum transfer 926 std::ostringstream ost2; 927 ost2 << path << "/penelope/rayleigh/MIFF/qex 928 file.open(ost2.str().c_str()); 929 930 if (!file.is_open()) { 931 G4String excep = "Data file " + G4String(o 932 G4Exception("G4PenelopeRayleighModelMI::Re 933 "em0003",FatalException,excep); 934 } 935 G4bool fillQGrid = false; 936 if (!fLogQSquareGrid.size()) { 937 fillQGrid = true; 938 nPoints = 1142; 939 } 940 G4double qext = 0; 941 if (fillQGrid) { //fLogQSquareGrid filled o 942 for (std::size_t i=0;i<nPoints;++i) { 943 file >> qext; 944 fLogQSquareGrid.push_back(2.0*G4Log(qext 945 } 946 } 947 file.close(); 948 949 //Finally, read the form factor file 950 std::ostringstream ost3; 951 if (Z>9) 952 ost3 << path << "/penelope/rayleigh/pdaff" 953 else 954 ost3 << path << "/penelope/rayleigh/pdaff0 955 file.open(ost3.str().c_str()); 956 957 if (!file.is_open()) { 958 G4String excep = "Data file " + G4String(o 959 G4Exception("G4PenelopeRayleighModelMI::Re 960 "em0003",FatalException,excep); 961 } 962 963 file >> readZ >> nPoints; 964 965 if (readZ != Z || nPoints <= 0 || nPoints >= 966 G4ExceptionDescription ed; 967 ed << "Corrupted data file for Z=" << Z << 968 G4Exception("G4PenelopeRayleighModelMI::Re 969 "em0005",FatalException,ed); 970 return; 971 } 972 973 fAtomicFormFactor[Z] = new G4PhysicsFreeVect 974 G4double q=0,ff=0,incoh=0; 975 for (std::size_t i=0;i<nPoints;++i) { 976 file >> q >> ff >> incoh; //q and ff are d 977 fAtomicFormFactor[Z]->PutValue(i,q,ff); 978 if (file.eof() && i != (nPoints-1)) { //fi 979 G4ExceptionDescription ed; 980 ed << "Corrupted data file for Z=" << Z 981 ed << "Found less than " << nPoints << " 982 G4Exception("G4PenelopeRayleighModelMI:: 983 "em0005",FatalException,ed); 984 } 985 } 986 file.close(); 987 return; 988 } 989 990 //....oooOO0OOooo........oooOO0OOooo........oo 991 992 void G4PenelopeRayleighModelMI::ReadMolInterfe 993 const G4String& FFfilename) 994 { 995 if (fVerboseLevel > 2) { 996 G4cout << "G4PenelopeRayleighModelMI::Read 997 matname << G4endl; 998 } 999 G4bool isLocalFile = (FFfilename != "NULL"); 1000 1001 const char* path = G4FindDataDir("G4LEDATA" 1002 if (!path) { 1003 G4String excep = "G4LEDATA environment va 1004 G4Exception("G4PenelopeRayleighModelMI::R 1005 "em0006",FatalException,excep); 1006 return; 1007 } 1008 1009 if (!(fKnownMaterials->count(matname)) && ! 1010 return; 1011 1012 G4String aFileName = (isLocalFile) ? FFfile 1013 1014 //if the material has a MIFF, read it from 1015 if (aFileName != "") { 1016 if (fVerboseLevel > 2) 1017 G4cout << "ReadMolInterferenceData(). R 1018 aFileName << " " << 1019 (isLocalFile ? "(local)" : "(database)") << 1020 std::ifstream file; 1021 std::ostringstream ostIMFF; 1022 if (isLocalFile) 1023 ostIMFF << aFileName; 1024 else 1025 ostIMFF << path << "/penelope/rayleigh/ 1026 file.open(ostIMFF.str().c_str()); 1027 1028 if (!file.is_open()) { 1029 G4String excep = "Data file " + G4Strin 1030 G4Exception("G4PenelopeRayleighModelMI: 1031 "em1031",FatalException,excep); 1032 return; 1033 } 1034 1035 //check the number of points in the file 1036 std::size_t nPoints = 0; 1037 G4double x=0,y=0; 1038 while (file.good()) { 1039 file >> x >> y; 1040 nPoints++; 1041 } 1042 file.close(); 1043 nPoints--; 1044 if (fVerboseLevel > 3) 1045 G4cout << "Number of nPoints: " << nPoi 1046 1047 //read the file 1048 file.open(ostIMFF.str().c_str()); 1049 G4PhysicsFreeVector* theFFVec = new G4Phy 1050 G4double q=0,ff=0; 1051 for (std::size_t i=0; i<nPoints; ++i) { 1052 file >> q >> ff; 1053 1054 //q and ff are dimensionless (q is in u 1055 theFFVec->PutValue(i,q,ff); 1056 1057 //file ended too early 1058 if (file.eof() && i != (nPoints-1)) { 1059 G4ExceptionDescription ed; 1060 ed << "Corrupted data file" << G4endl; 1061 ed << "Found less than " << nPoints << " en 1062 G4Exception("G4PenelopeRayleighModelMI::Rea 1063 "em1005",FatalException,ed); 1064 } 1065 } 1066 if (!fMolInterferenceData) { 1067 G4Exception("G4PenelopeRayleighModelMI: 1068 "em2145",FatalException, 1069 "Unable to allocate the Molecular Inter 1070 delete theFFVec; 1071 return; 1072 } 1073 file.close(); 1074 fMolInterferenceData->insert(std::make_pa 1075 } 1076 return; 1077 } 1078 1079 //....oooOO0OOooo........oooOO0OOooo........o 1080 1081 G4double G4PenelopeRayleighModelMI::GetFSquar 1082 { 1083 G4double f2 = 0; 1084 //Input value QSquared could be zero: prote 1085 //the FPE exception 1086 1087 //If Q<1e-10, set Q to 1e-10 1088 G4double logQSquared = (QSquared>1e-10) ? G 1089 //last value of the table 1090 G4double maxlogQ2 = fLogQSquareGrid[fLogQSq 1091 1092 //now it should be all right 1093 G4PhysicsFreeVector* theVec = fLogFormFacto 1094 1095 if (!theVec) { 1096 G4ExceptionDescription ed; 1097 ed << "Unable to retrieve F squared table 1098 G4Exception("G4PenelopeRayleighModelMI::G 1099 "em2046",FatalException,ed); 1100 return 0; 1101 } 1102 1103 if (logQSquared < -20) { //Q < 1e-9 1104 G4double logf2 = (*theVec)[0]; //first va 1105 f2 = G4Exp(logf2); 1106 } 1107 else if (logQSquared > maxlogQ2) 1108 f2 =0; 1109 else { 1110 //log(Q^2) vs. log(F^2) 1111 G4double logf2 = theVec->Value(logQSquare 1112 f2 = G4Exp(logf2); 1113 } 1114 1115 if (fVerboseLevel > 3) { 1116 G4cout << "G4PenelopeRayleighModelMI::Get 1117 G4cout << "Q^2 = " << QSquared << " (uni 1118 } 1119 return f2; 1120 } 1121 1122 //....oooOO0OOooo........oooOO0OOooo........o 1123 1124 void G4PenelopeRayleighModelMI::InitializeSam 1125 { 1126 G4double q2min = 0; 1127 G4double q2max = 0; 1128 const std::size_t np = 150; //hard-coded in 1129 for (std::size_t i=1;i<fLogQSquareGrid.size 1130 { 1131 G4double Q2 = G4Exp(fLogQSquareGrid[i]) 1132 if (GetFSquared(mat,Q2) > 1e-35) 1133 { 1134 q2max = G4Exp(fLogQSquareGrid[i-1]); 1135 } 1136 } 1137 std::size_t nReducedPoints = np/4; 1138 1139 //check for errors 1140 if (np < 16) 1141 { 1142 G4Exception("G4PenelopeRayleighModelMI: 1143 "em2047",FatalException, 1144 "Too few points to initialize the sampl 1145 } 1146 if (q2min > (q2max-1e-10)) 1147 { 1148 G4cout << "q2min= " << q2min << " q2max 1149 G4Exception("G4PenelopeRayleighModelMI: 1150 "em2048",FatalException, 1151 "Too narrow grid to initialize the samp 1152 } 1153 1154 //This is subroutine RITAI0 of Penelope 1155 //Create an object of type G4PenelopeRaylei 1156 1157 //temporary vectors --> Then everything is 1158 G4DataVector* x = new G4DataVector(); 1159 1160 /****************************************** 1161 Start with a grid of NUNIF points uniform 1162 ******************************************* 1163 std::size_t NUNIF = std::min(std::max(((std 1164 const G4int nip = 51; //hard-coded in Penel 1165 1166 G4double dx = (q2max-q2min)/((G4double) NUN 1167 x->push_back(q2min); 1168 for (std::size_t i=1;i<NUNIF-1;++i) 1169 { 1170 G4double app = q2min + i*dx; 1171 x->push_back(app); //increase 1172 } 1173 x->push_back(q2max); 1174 1175 if (fVerboseLevel> 3) 1176 G4cout << "Vector x has " << x->size() << 1177 1178 //Allocate temporary storage vectors 1179 G4DataVector* area = new G4DataVector(); 1180 G4DataVector* a = new G4DataVector(); 1181 G4DataVector* b = new G4DataVector(); 1182 G4DataVector* c = new G4DataVector(); 1183 G4DataVector* err = new G4DataVector(); 1184 1185 for (std::size_t i=0;i<NUNIF-1;++i) //build 1186 { 1187 //Temporary vectors for this loop 1188 G4DataVector* pdfi = new G4DataVector() 1189 G4DataVector* pdfih = new G4DataVector( 1190 G4DataVector* sumi = new G4DataVector() 1191 1192 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4d 1193 G4double pdfmax = 0; 1194 for (G4int k=0;k<nip;k++) 1195 { 1196 G4double xik = (*x)[i]+k*dxi; 1197 G4double pdfk = std::max(GetFSquared(mat, 1198 pdfi->push_back(pdfk); 1199 pdfmax = std::max(pdfmax,pdfk); 1200 if (k < (nip-1)) 1201 { 1202 G4double xih = xik + 0.5*dxi; 1203 G4double pdfIK = std::max(GetFSquared 1204 pdfih->push_back(pdfIK); 1205 pdfmax = std::max(pdfmax,pdfIK); 1206 } 1207 } 1208 1209 //Simpson's integration 1210 G4double cons = dxi*0.5*(1./3.); 1211 sumi->push_back(0.); 1212 for (G4int k=1;k<nip;k++) 1213 { 1214 G4double previous = (*sumi)[k-1]; 1215 G4double next = previous + cons*((*pdfi)[ 1216 sumi->push_back(next); 1217 } 1218 1219 G4double lastIntegral = (*sumi)[sumi->s 1220 area->push_back(lastIntegral); 1221 //Normalize cumulative function 1222 G4double factor = 1.0/lastIntegral; 1223 for (std::size_t k=0;k<sumi->size();++k 1224 (*sumi)[k] *= factor; 1225 1226 //When the PDF vanishes at one of the i 1227 if ((*pdfi)[0] < 1e-35) 1228 (*pdfi)[0] = 1e-5*pdfmax; 1229 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1230 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1231 1232 G4double pli = (*pdfi)[0]*factor; 1233 G4double pui = (*pdfi)[pdfi->size()-1]* 1234 G4double B_temp = 1.0-1.0/(pli*pui*dx*d 1235 G4double A_temp = (1.0/(pli*dx))-1.0-B_ 1236 G4double C_temp = 1.0+A_temp+B_temp; 1237 if (C_temp < 1e-35) 1238 { 1239 a->push_back(0.); 1240 b->push_back(0.); 1241 c->push_back(1.); 1242 } 1243 else 1244 { 1245 a->push_back(A_temp); 1246 b->push_back(B_temp); 1247 c->push_back(C_temp); 1248 } 1249 1250 //OK, now get ERR(I), the integral of t 1251 //and the true pdf, extended over the i 1252 G4int icase = 1; //loop code 1253 G4bool reLoop = false; 1254 err->push_back(0.); 1255 do 1256 { 1257 reLoop = false; 1258 (*err)[i] = 0.; //zero variable 1259 for (G4int k=0;k<nip;k++) 1260 { 1261 G4double rr = (*sumi)[k]; 1262 G4double pap = (*area)[i]*(1.0+((*a)[ 1263 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-( 1264 if (k == 0 || k == nip-1) 1265 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k] 1266 else 1267 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 1268 } 1269 (*err)[i] *= dxi; 1270 1271 //If err(I) is too large, the pdf is appr 1272 if ((*err)[i] > 0.1*(*area)[i] && icase = 1273 { 1274 (*b)[i] = 0; 1275 (*a)[i] = 0; 1276 (*c)[i] = 1.; 1277 icase = 2; 1278 reLoop = true; 1279 } 1280 }while(reLoop); 1281 1282 delete pdfi; 1283 delete pdfih; 1284 delete sumi; 1285 } //end of first loop over i 1286 1287 //Now assign last point 1288 (*x)[x->size()-1] = q2max; 1289 a->push_back(0.); 1290 b->push_back(0.); 1291 c->push_back(0.); 1292 err->push_back(0.); 1293 area->push_back(0.); 1294 1295 if (x->size() != NUNIF || a->size() != NUNI 1296 err->size() != NUNIF || area->size() != 1297 { 1298 G4ExceptionDescription ed; 1299 ed << "Problem in building the Table fo 1300 G4Exception("G4PenelopeRayleighModelMI: 1301 "em2049",FatalException,ed); 1302 } 1303 1304 /****************************************** 1305 New grid points are added by halving the s 1306 This is done up to np=150 points in the gri 1307 ******************************************* 1308 do 1309 { 1310 G4double maxError = 0.0; 1311 std::size_t iErrMax = 0; 1312 for (std::size_t i=0;i<err->size()-2;++ 1313 { 1314 //maxError is the lagest of the interval 1315 if ((*err)[i] > maxError) 1316 { 1317 maxError = (*err)[i]; 1318 iErrMax = i; 1319 } 1320 } 1321 1322 //OK, now I have to insert one new poin 1323 G4double newx = 0.5*((*x)[iErrMax]+(*x) 1324 1325 x->insert(x->begin()+iErrMax+1,newx); 1326 //Add place-holders in the other vector 1327 area->insert(area->begin()+iErrMax+1,0. 1328 a->insert(a->begin()+iErrMax+1,0.); 1329 b->insert(b->begin()+iErrMax+1,0.); 1330 c->insert(c->begin()+iErrMax+1,0.); 1331 err->insert(err->begin()+iErrMax+1,0.); 1332 1333 //Now calculate the other parameters 1334 for (std::size_t i=iErrMax;i<=iErrMax+1 1335 { 1336 //Temporary vectors for this loop 1337 G4DataVector* pdfi = new G4DataVector(); 1338 G4DataVector* pdfih = new G4DataVector(); 1339 G4DataVector* sumi = new G4DataVector(); 1340 1341 G4double dxLocal = (*x)[i+1]-(*x)[i]; 1342 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4dou 1343 G4double pdfmax = 0; 1344 for (G4int k=0;k<nip;k++) 1345 { 1346 G4double xik = (*x)[i]+k*dxi; 1347 G4double pdfk = std::max(GetFSquared( 1348 pdfi->push_back(pdfk); 1349 pdfmax = std::max(pdfmax,pdfk); 1350 if (k < (nip-1)) 1351 { 1352 G4double xih = xik + 0.5*dxi; 1353 G4double pdfIK = std::max(GetFSquared(m 1354 pdfih->push_back(pdfIK); 1355 pdfmax = std::max(pdfmax,pdfIK); 1356 } 1357 } 1358 1359 //Simpson's integration 1360 G4double cons = dxi*0.5*(1./3.); 1361 sumi->push_back(0.); 1362 for (G4int k=1;k<nip;k++) 1363 { 1364 G4double previous = (*sumi)[k-1]; 1365 G4double next = previous + cons*((*pd 1366 sumi->push_back(next); 1367 } 1368 G4double lastIntegral = (*sumi)[sumi->siz 1369 (*area)[i] = lastIntegral; 1370 1371 //Normalize cumulative function 1372 G4double factor = 1.0/lastIntegral; 1373 for (std::size_t k=0;k<sumi->size();++k) 1374 (*sumi)[k] *= factor; 1375 1376 //When the PDF vanishes at one of the int 1377 if ((*pdfi)[0] < 1e-35) 1378 (*pdfi)[0] = 1e-5*pdfmax; 1379 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1380 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1381 1382 G4double pli = (*pdfi)[0]*factor; 1383 G4double pui = (*pdfi)[pdfi->size()-1]*fa 1384 G4double B_temp = 1.0-1.0/(pli*pui*dxLoca 1385 G4double A_temp = (1.0/(pli*dxLocal))-1.0 1386 G4double C_temp = 1.0+A_temp+B_temp; 1387 if (C_temp < 1e-35) 1388 { 1389 (*a)[i]= 0.; 1390 (*b)[i] = 0.; 1391 (*c)[i] = 1; 1392 } 1393 else 1394 { 1395 (*a)[i]= A_temp; 1396 (*b)[i] = B_temp; 1397 (*c)[i] = C_temp; 1398 } 1399 //OK, now get ERR(I), the integral of the 1400 //and the true pdf, extended over the int 1401 G4int icase = 1; //loop code 1402 G4bool reLoop = false; 1403 do 1404 { 1405 reLoop = false; 1406 (*err)[i] = 0.; //zero variable 1407 for (G4int k=0;k<nip;k++) 1408 { 1409 G4double rr = (*sumi)[k]; 1410 G4double pap = (*area)[i]*(1.0+((*a)[i] 1411 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+ 1412 if (k == 0 || k == nip-1) 1413 (*err)[i] += 0.5*std::fabs(pap-(*pdfi 1414 else 1415 (*err)[i] += std::fabs(pap-(*pdfi)[k] 1416 } 1417 (*err)[i] *= dxi; 1418 1419 //If err(I) is too large, the pdf is 1420 if ((*err)[i] > 0.1*(*area)[i] && ica 1421 { 1422 (*b)[i] = 0; 1423 (*a)[i] = 0; 1424 (*c)[i] = 1.; 1425 icase = 2; 1426 reLoop = true; 1427 } 1428 }while(reLoop); 1429 delete pdfi; 1430 delete pdfih; 1431 delete sumi; 1432 } 1433 }while(x->size() < np); 1434 1435 if (x->size() != np || a->size() != np || 1436 err->size() != np || area->size() != np 1437 { 1438 G4Exception("G4PenelopeRayleighModelMI: 1439 "em2050",FatalException, 1440 "Problem in building the extended Table 1441 } 1442 1443 /****************************************** 1444 Renormalization 1445 ******************************************* 1446 G4double ws = 0; 1447 for (std::size_t i=0;i<np-1;++i) 1448 ws += (*area)[i]; 1449 ws = 1.0/ws; 1450 G4double errMax = 0; 1451 for (std::size_t i=0;i<np-1;++i) 1452 { 1453 (*area)[i] *= ws; 1454 (*err)[i] *= ws; 1455 errMax = std::max(errMax,(*err)[i]); 1456 } 1457 1458 //Vector with the normalized cumulative dis 1459 G4DataVector* PAC = new G4DataVector(); 1460 PAC->push_back(0.); 1461 for (std::size_t i=0;i<np-1;++i) 1462 { 1463 G4double previous = (*PAC)[i]; 1464 PAC->push_back(previous+(*area)[i]); 1465 } 1466 (*PAC)[PAC->size()-1] = 1.; 1467 1468 /****************************************** 1469 Pre-calculated limits for the initial binar 1470 ******************************************* 1471 std::vector<std::size_t> *ITTL = new std::v 1472 std::vector<std::size_t> *ITTU = new std::v 1473 1474 //Just create place-holders 1475 for (std::size_t i=0;i<np;++i) 1476 { 1477 ITTL->push_back(0); 1478 ITTU->push_back(0); 1479 } 1480 1481 G4double bin = 1.0/(np-1); 1482 (*ITTL)[0]=0; 1483 for (std::size_t i=1;i<(np-1);++i) 1484 { 1485 G4double ptst = i*bin; 1486 G4bool found = false; 1487 for (std::size_t j=(*ITTL)[i-1];j<np && 1488 { 1489 if ((*PAC)[j] > ptst) 1490 { 1491 (*ITTL)[i] = j-1; 1492 (*ITTU)[i-1] = j; 1493 found = true; 1494 } 1495 } 1496 } 1497 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; 1498 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; 1499 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; 1500 1501 if (ITTU->size() != np || ITTU->size() != n 1502 { 1503 G4Exception("G4PenelopeRayleighModelMI: 1504 "em2051",FatalException, 1505 "Problem in building the Limit Tables f 1506 } 1507 1508 /****************************************** 1509 Copy tables 1510 ******************************************* 1511 G4PenelopeSamplingData* theTable = new G4Pe 1512 for (std::size_t i=0;i<np;++i) 1513 { 1514 theTable->AddPoint((*x)[i],(*PAC)[i],(* 1515 } 1516 if (fVerboseLevel > 2) 1517 { 1518 G4cout << "**************************** 1519 G4endl; 1520 G4cout << "Sampling table for Penelope 1521 theTable->DumpTable(); 1522 } 1523 fSamplingTable->insert(std::make_pair(mat,t 1524 1525 //Clean up temporary vectors 1526 delete x; 1527 delete a; 1528 delete b; 1529 delete c; 1530 delete err; 1531 delete area; 1532 delete PAC; 1533 delete ITTL; 1534 delete ITTU; 1535 1536 //DONE! 1537 return; 1538 } 1539 1540 //....oooOO0OOooo........oooOO0OOooo........o 1541 1542 void G4PenelopeRayleighModelMI::GetPMaxTable( 1543 { 1544 if (!fPMaxTable) 1545 { 1546 G4cout << "G4PenelopeRayleighModelMI::B 1547 G4cout << "Going to instanziate the fPM 1548 G4cout << "That should _not_ be here! " 1549 fPMaxTable = new std::map<const G4Mater 1550 } 1551 //check if the table is already there 1552 if (fPMaxTable->count(mat)) 1553 return; 1554 1555 //otherwise build it 1556 if (!fSamplingTable) 1557 { 1558 G4Exception("G4PenelopeRayleighModelMI: 1559 "em2052",FatalException, 1560 "SamplingTable is not properly instanti 1561 return; 1562 } 1563 1564 //This should not be: the sampling table is 1565 if (!fSamplingTable->count(mat)) 1566 { 1567 G4ExceptionDescription ed; 1568 ed << "Sampling table for material " << 1569 G4Exception("G4PenelopeRayleighModelMI: 1570 "em2052",FatalException, 1571 ed); 1572 return; 1573 } 1574 1575 G4PenelopeSamplingData *theTable = fSamplin 1576 std::size_t tablePoints = theTable->GetNumb 1577 std::size_t nOfEnergyPoints = fLogEnergyGri 1578 G4PhysicsFreeVector* theVec = new G4Physics 1579 1580 const std::size_t nip = 51; //hard-coded in 1581 1582 for (std::size_t ie=0;ie<fLogEnergyGridPMax 1583 { 1584 G4double energy = G4Exp(fLogEnergyGridP 1585 G4double Qm = 2.0*energy/electron_mass_ 1586 G4double Qm2 = Qm*Qm; 1587 G4double firstQ2 = theTable->GetX(0); 1588 G4double lastQ2 = theTable->GetX(tableP 1589 G4double thePMax = 0; 1590 1591 if (Qm2 > firstQ2) 1592 { 1593 if (Qm2 < lastQ2) 1594 { 1595 //bisection to look for the index of 1596 std::size_t lowerBound = 0; 1597 std::size_t upperBound = tablePoints- 1598 while (lowerBound <= upperBound) 1599 { 1600 std::size_t midBin = (lowerBound + uppe 1601 if( Qm2 < theTable->GetX(midBin)) 1602 { upperBound = midBin-1; } 1603 else 1604 { lowerBound = midBin+1; } 1605 } 1606 //upperBound is the output (but also 1607 G4double Q1 = theTable->GetX(upperBou 1608 G4double Q2 = Qm2; 1609 G4double DQ = (Q2-Q1)/((G4double)(nip 1610 G4double theA = theTable->GetA(upperB 1611 G4double theB = theTable->GetB(upperB 1612 G4double thePAC = theTable->GetPAC(up 1613 G4DataVector* fun = new G4DataVector( 1614 for (std::size_t k=0;k<nip;++k) 1615 { 1616 G4double qi = Q1 + k*DQ; 1617 G4double tau = (qi-Q1)/ 1618 (theTable->GetX(upperBound+1)-Q1); 1619 G4double con1 = 2.0*theB*tau; 1620 G4double ci = 1.0+theA+theB; 1621 G4double con2 = ci-theA*tau; 1622 G4double etap = 0; 1623 if (std::fabs(con1) > 1.0e-16*std::fabs 1624 etap = con2*(1.0-std::sqrt(1.0-2.0*ta 1625 else 1626 etap = tau/con2; 1627 G4double theFun = (theTable->GetPAC(upp 1628 (1.0+(theA+theB*etap)*etap)*(1.0+(the 1629 ((1.0-theB*etap*etap)*ci*(theTable->G 1630 fun->push_back(theFun); 1631 } 1632 //Now intergrate numerically the fun 1633 G4DataVector* sum = new G4DataVector; 1634 G4double CONS = DQ*(1./12.); 1635 G4double HCONS = 0.5*CONS; 1636 sum->push_back(0.); 1637 G4double secondPoint = (*sum)[0] + 1638 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C 1639 sum->push_back(secondPoint); 1640 for (std::size_t hh=2;hh<nip-1;++hh) 1641 { 1642 G4double previous = (*sum)[hh-1]; 1643 G4double next = previous+(13.0*((*fun)[ 1644 (*fun)[hh+1]-(*fun)[hh-2])*HCON 1645 sum->push_back(next); 1646 } 1647 G4double last = (*sum)[nip-2]+(5.0*(* 1648 (*fun)[nip-3])*CONS; 1649 sum->push_back(last); 1650 thePMax = thePAC + (*sum)[sum->size() 1651 delete fun; 1652 delete sum; 1653 } 1654 else 1655 { 1656 thePMax = 1.0; 1657 } 1658 } 1659 else 1660 { 1661 thePMax = theTable->GetPAC(0); 1662 } 1663 1664 //Write number in the table 1665 theVec->PutValue(ie,energy,thePMax); 1666 } 1667 1668 fPMaxTable->insert(std::make_pair(mat,theVe 1669 return; 1670 } 1671 1672 //....oooOO0OOooo........oooOO0OOooo........o 1673 1674 void G4PenelopeRayleighModelMI::DumpFormFacto 1675 { 1676 G4cout << "******************************** 1677 G4cout << "G4PenelopeRayleighModelMI: Form 1678 //try to use the same format as Penelope-Fo 1679 G4cout << "Q/(m_e*c) F(Q) 1680 G4cout << "******************************** 1681 if (!fLogFormFactorTable->count(mat)) 1682 BuildFormFactorTable(mat); 1683 1684 G4PhysicsFreeVector* theVec = fLogFormFacto 1685 for (std::size_t i=0;i<theVec->GetVectorLen 1686 { 1687 G4double logQ2 = theVec->GetLowEdgeEner 1688 G4double Q = G4Exp(0.5*logQ2); 1689 G4double logF2 = (*theVec)[i]; 1690 G4double F = G4Exp(0.5*logF2); 1691 G4cout << Q << " " << F << 1692 } 1693 //DONE 1694 return; 1695 } 1696 1697 //....oooOO0OOooo........oooOO0OOooo........o 1698 1699 void G4PenelopeRayleighModelMI::SetParticle(c 1700 { 1701 if(!fParticle) { 1702 fParticle = p; 1703 } 1704 } 1705 1706 //....oooOO0OOooo........oooOO0OOooo........o 1707 void G4PenelopeRayleighModelMI::LoadKnownMIFF 1708 { 1709 fKnownMaterials->insert(std::pair<G4String, 1710 fKnownMaterials->insert(std::pair<G4String, 1711 fKnownMaterials->insert(std::pair<G4String, 1712 fKnownMaterials->insert(std::pair<G4String, 1713 fKnownMaterials->insert(std::pair<G4String, 1714 fKnownMaterials->insert(std::pair<G4String, 1715 fKnownMaterials->insert(std::pair<G4String, 1716 fKnownMaterials->insert(std::pair<G4String, 1717 fKnownMaterials->insert(std::pair<G4String, 1718 fKnownMaterials->insert(std::pair<G4String, 1719 fKnownMaterials->insert(std::pair<G4String, 1720 fKnownMaterials->insert(std::pair<G4String, 1721 fKnownMaterials->insert(std::pair<G4String, 1722 fKnownMaterials->insert(std::pair<G4String, 1723 fKnownMaterials->insert(std::pair<G4String, 1724 fKnownMaterials->insert(std::pair<G4String, 1725 fKnownMaterials->insert(std::pair<G4String, 1726 fKnownMaterials->insert(std::pair<G4String, 1727 fKnownMaterials->insert(std::pair<G4String, 1728 fKnownMaterials->insert(std::pair<G4String, 1729 fKnownMaterials->insert(std::pair<G4String, 1730 fKnownMaterials->insert(std::pair<G4String, 1731 fKnownMaterials->insert(std::pair<G4String, 1732 fKnownMaterials->insert(std::pair<G4String, 1733 fKnownMaterials->insert(std::pair<G4String, 1734 fKnownMaterials->insert(std::pair<G4String, 1735 fKnownMaterials->insert(std::pair<G4String, 1736 fKnownMaterials->insert(std::pair<G4String, 1737 fKnownMaterials->insert(std::pair<G4String, 1738 fKnownMaterials->insert(std::pair<G4String, 1739 fKnownMaterials->insert(std::pair<G4String, 1740 fKnownMaterials->insert(std::pair<G4String, 1741 fKnownMaterials->insert(std::pair<G4String, 1742 } 1743 1744 //....oooOO0OOooo........oooOO0OOooo........o 1745 1746 G4double G4PenelopeRayleighModelMI::Integrate 1747 { 1748 G4double integral = 0.; 1749 for (G4int k=0; k<n-1; k++) { 1750 integral += (y[k]+y[k+1]); 1751 } 1752 integral *= dTheta*0.5; 1753 return integral; 1754 } 1755 1756 //....oooOO0OOooo........oooOO0OOooo........o 1757 G4MIData* G4PenelopeRayleighModelMI::GetMIDat 1758 { 1759 if (material->IsExtended()) { 1760 G4ExtendedMaterial* aEM = (G4ExtendedMate 1761 G4MIData* dataMI = (G4MIData*)aEM->Retrie 1762 return dataMI; 1763 } else { 1764 return nullptr; 1765 } 1766 } 1767