Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4PenelopeRayleighModel.cc 97613 2016-06-06 12:24:51Z gcosmo $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 03 Dec 2009 L Pandola First implementa 32 // 03 Dec 2009 L Pandola First implementation 32 // 25 May 2011 L.Pandola Renamed (make v2 33 // 25 May 2011 L.Pandola Renamed (make v2008 as default Penelope) 33 // 19 Sep 2013 L.Pandola Migration to MT 34 // 19 Sep 2013 L.Pandola Migration to MT 34 // 35 // 35 36 36 #include "G4PenelopeRayleighModel.hh" 37 #include "G4PenelopeRayleighModel.hh" 37 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh" 39 #include "G4PenelopeSamplingData.hh" 40 #include "G4PenelopeSamplingData.hh" 40 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh" 41 #include "G4MaterialCutsCouple.hh" 42 #include "G4MaterialCutsCouple.hh" 42 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh" 43 #include "G4DynamicParticle.hh" 44 #include "G4DynamicParticle.hh" 44 #include "G4PhysicsTable.hh" 45 #include "G4PhysicsTable.hh" 45 #include "G4ElementTable.hh" 46 #include "G4ElementTable.hh" 46 #include "G4Element.hh" 47 #include "G4Element.hh" 47 #include "G4PhysicsFreeVector.hh" 48 #include "G4PhysicsFreeVector.hh" 48 #include "G4AutoLock.hh" 49 #include "G4AutoLock.hh" 49 #include "G4Exp.hh" 50 #include "G4Exp.hh" 50 51 51 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 52 53 53 const G4int G4PenelopeRayleighModel::fMaxZ; << 54 G4PhysicsFreeVector* G4PenelopeRayleighModel:: << 55 G4PhysicsFreeVector* G4PenelopeRayleighModel:: << 56 << 57 //....oooOO0OOooo........oooOO0OOooo........oo << 58 << 59 G4PenelopeRayleighModel::G4PenelopeRayleighMod 54 G4PenelopeRayleighModel::G4PenelopeRayleighModel(const G4ParticleDefinition* part, 60 const G4String& nam) 55 const G4String& nam) 61 :G4VEmModel(nam),fParticleChange(nullptr),fP << 56 :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false), 62 fLogFormFactorTable(nullptr),fPMaxTable(nul << 57 logAtomicCrossSection(0),atomicFormFactor(0),logFormFactorTable(0), 63 fIsInitialised(false),fLocalTable(false) << 58 pMaxTable(0),samplingTable(0),fLocalTable(false) 64 { 59 { 65 fIntrinsicLowEnergyLimit = 100.0*eV; 60 fIntrinsicLowEnergyLimit = 100.0*eV; 66 fIntrinsicHighEnergyLimit = 100.0*GeV; 61 fIntrinsicHighEnergyLimit = 100.0*GeV; 67 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 62 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 68 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 63 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 69 64 70 if (part) 65 if (part) 71 SetParticle(part); 66 SetParticle(part); 72 67 73 // 68 // 74 fVerboseLevel= 0; << 69 verboseLevel= 0; 75 // Verbosity scale: 70 // Verbosity scale: 76 // 0 = nothing 71 // 0 = nothing 77 // 1 = warning for energy non-conservation 72 // 1 = warning for energy non-conservation 78 // 2 = details of energy budget 73 // 2 = details of energy budget 79 // 3 = calculation of cross sections, file o 74 // 3 = calculation of cross sections, file openings, sampling of atoms 80 // 4 = entering in methods 75 // 4 = entering in methods 81 76 82 //build the energy grid. It is the same for 77 //build the energy grid. It is the same for all materials 83 G4double logenergy = G4Log(fIntrinsicLowEner << 78 G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.); 84 G4double logmaxenergy = G4Log(1.5*fIntrinsic << 79 G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit); 85 //finer grid below 160 keV 80 //finer grid below 160 keV 86 G4double logtransitionenergy = G4Log(160*keV << 81 G4double logtransitionenergy = std::log(160*keV); 87 G4double logfactor1 = G4Log(10.)/250.; << 82 G4double logfactor1 = std::log(10.)/250.; 88 G4double logfactor2 = logfactor1*10; 83 G4double logfactor2 = logfactor1*10; 89 fLogEnergyGridPMax.push_back(logenergy); << 84 logEnergyGridPMax.push_back(logenergy); 90 do{ 85 do{ 91 if (logenergy < logtransitionenergy) 86 if (logenergy < logtransitionenergy) 92 logenergy += logfactor1; 87 logenergy += logfactor1; 93 else 88 else 94 logenergy += logfactor2; 89 logenergy += logfactor2; 95 fLogEnergyGridPMax.push_back(logenergy); << 90 logEnergyGridPMax.push_back(logenergy); 96 }while (logenergy < logmaxenergy); 91 }while (logenergy < logmaxenergy); 97 } 92 } 98 93 99 //....oooOO0OOooo........oooOO0OOooo........oo 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 100 95 101 G4PenelopeRayleighModel::~G4PenelopeRayleighMo 96 G4PenelopeRayleighModel::~G4PenelopeRayleighModel() 102 { 97 { 103 if (IsMaster() || fLocalTable) 98 if (IsMaster() || fLocalTable) 104 { 99 { 105 << 100 std::map <G4int,G4PhysicsFreeVector*>::iterator i; 106 for(G4int i=0; i<=fMaxZ; ++i) << 101 if (logAtomicCrossSection) 107 { 102 { 108 if(fLogAtomicCrossSection[i]) << 103 for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++) 109 { << 104 if (i->second) delete i->second; 110 delete fLogAtomicCrossSection[i]; << 105 111 fLogAtomicCrossSection[i] = nullptr; << 106 delete logAtomicCrossSection; 112 } << 107 logAtomicCrossSection = 0; 113 if(fAtomicFormFactor[i]) << 114 { << 115 delete fAtomicFormFactor[i]; << 116 fAtomicFormFactor[i] = nullptr; << 117 } << 118 } 108 } >> 109 if (atomicFormFactor) >> 110 { >> 111 for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++) >> 112 if (i->second) delete i->second; >> 113 delete atomicFormFactor; >> 114 atomicFormFactor = 0; >> 115 } >> 116 119 ClearTables(); 117 ClearTables(); 120 } 118 } 121 } 119 } 122 120 123 //....oooOO0OOooo........oooOO0OOooo........oo 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 124 void G4PenelopeRayleighModel::ClearTables() 122 void G4PenelopeRayleighModel::ClearTables() 125 { 123 { 126 if (fLogFormFactorTable) << 124 /* >> 125 if (!IsMaster()) >> 126 //Should not be here! >> 127 G4Exception("G4PenelopeRayleighModel::ClearTables()", >> 128 "em0100",FatalException,"Worker thread in this method"); >> 129 */ >> 130 >> 131 std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i; >> 132 >> 133 if (logFormFactorTable) 127 { 134 { 128 for (auto& item : (*fLogFormFactorTable << 135 for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++) 129 if (item.second) delete item.second; << 136 if (i->second) delete i->second; 130 delete fLogFormFactorTable; << 137 delete logFormFactorTable; 131 fLogFormFactorTable = nullptr; //zero e << 138 logFormFactorTable = 0; //zero explicitely 132 } 139 } 133 if (fPMaxTable) << 140 >> 141 if (pMaxTable) 134 { 142 { 135 for (auto& item : (*fPMaxTable)) << 143 for (i=pMaxTable->begin(); i != pMaxTable->end(); i++) 136 if (item.second) delete item.second; << 144 if (i->second) delete i->second; 137 delete fPMaxTable; << 145 delete pMaxTable; 138 fPMaxTable = nullptr; //zero explicitly << 146 pMaxTable = 0; //zero explicitely 139 } 147 } 140 if (fSamplingTable) << 148 >> 149 std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii; >> 150 if (samplingTable) 141 { 151 { 142 for (auto& item : (*fSamplingTable)) << 152 for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++) 143 if (item.second) delete item.second; << 153 if (ii->second) delete ii->second; 144 delete fSamplingTable; << 154 delete samplingTable; 145 fSamplingTable = nullptr; //zero explic << 155 samplingTable = 0; //zero explicitely 146 } 156 } >> 157 147 return; 158 return; 148 } 159 } 149 160 150 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 151 162 152 void G4PenelopeRayleighModel::Initialise(const 163 void G4PenelopeRayleighModel::Initialise(const G4ParticleDefinition* part, 153 const G4DataVector& ) 164 const G4DataVector& ) 154 { 165 { 155 if (fVerboseLevel > 3) << 166 if (verboseLevel > 3) 156 G4cout << "Calling G4PenelopeRayleighModel 167 G4cout << "Calling G4PenelopeRayleighModel::Initialise()" << G4endl; 157 168 158 SetParticle(part); 169 SetParticle(part); 159 170 160 //Only the master model creates/fills/destro 171 //Only the master model creates/fills/destroys the tables 161 if (IsMaster() && part == fParticle) 172 if (IsMaster() && part == fParticle) 162 { 173 { 163 //clear tables depending on materials, n 174 //clear tables depending on materials, not the atomic ones 164 ClearTables(); 175 ClearTables(); 165 176 166 if (fVerboseLevel > 3) << 177 if (verboseLevel > 3) 167 G4cout << "Calling G4PenelopeRayleighModel:: 178 G4cout << "Calling G4PenelopeRayleighModel::Initialise() [master]" << G4endl; 168 179 169 //create new tables << 180 //create new tables 170 if (!fLogFormFactorTable) << 181 // 171 fLogFormFactorTable = new std::map<const G4M << 182 // logAtomicCrossSection and atomicFormFactor are created only once, 172 if (!fPMaxTable) << 183 // since they are never cleared 173 fPMaxTable = new std::map<const G4Material*, << 184 if (!logAtomicCrossSection) 174 if (!fSamplingTable) << 185 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>; 175 fSamplingTable = new std::map<const G4Materi << 186 if (!atomicFormFactor) >> 187 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>; >> 188 >> 189 if (!logFormFactorTable) >> 190 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; >> 191 if (!pMaxTable) >> 192 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; >> 193 if (!samplingTable) >> 194 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>; >> 195 176 196 177 G4ProductionCutsTable* theCoupleTable = 197 G4ProductionCutsTable* theCoupleTable = 178 G4ProductionCutsTable::GetProductionCutsTabl 198 G4ProductionCutsTable::GetProductionCutsTable(); 179 199 180 for (G4int i=0;i<(G4int)theCoupleTable-> << 200 for (size_t i=0;i<theCoupleTable->GetTableSize();i++) 181 { 201 { 182 const G4Material* material = 202 const G4Material* material = 183 theCoupleTable->GetMaterialCutsCouple(i) 203 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 184 const G4ElementVector* theElementVector = 204 const G4ElementVector* theElementVector = material->GetElementVector(); 185 205 186 for (std::size_t j=0;j<material->GetNumber << 206 for (size_t j=0;j<material->GetNumberOfElements();j++) 187 { 207 { 188 G4int iZ = theElementVector->at(j)->Ge << 208 G4int iZ = (G4int) theElementVector->at(j)->GetZ(); 189 //read data files only in the master 209 //read data files only in the master 190 if (!fLogAtomicCrossSection[iZ]) << 210 if (!logAtomicCrossSection->count(iZ)) 191 ReadDataFile(iZ); 211 ReadDataFile(iZ); 192 } 212 } 193 213 194 //1) If the table has not been built for t 214 //1) If the table has not been built for the material, do it! 195 if (!fLogFormFactorTable->count(material)) << 215 if (!logFormFactorTable->count(material)) 196 BuildFormFactorTable(material); 216 BuildFormFactorTable(material); 197 217 198 //2) retrieve or build the sampling table 218 //2) retrieve or build the sampling table 199 if (!(fSamplingTable->count(material))) << 219 if (!(samplingTable->count(material))) 200 InitializeSamplingAlgorithm(material); 220 InitializeSamplingAlgorithm(material); 201 221 202 //3) retrieve or build the pMax data 222 //3) retrieve or build the pMax data 203 if (!fPMaxTable->count(material)) << 223 if (!pMaxTable->count(material)) 204 GetPMaxTable(material); 224 GetPMaxTable(material); >> 225 205 } 226 } 206 227 207 if (fVerboseLevel > 1) { << 228 if (verboseLevel > 1) { 208 G4cout << "Penelope Rayleigh model v2008 is 229 G4cout << "Penelope Rayleigh model v2008 is initialized " << G4endl 209 << "Energy range: " 230 << "Energy range: " 210 << LowEnergyLimit() / keV << " keV - 231 << LowEnergyLimit() / keV << " keV - " 211 << HighEnergyLimit() / GeV << " GeV" 232 << HighEnergyLimit() / GeV << " GeV" 212 << G4endl; 233 << G4endl; 213 } 234 } 214 } 235 } 215 236 216 if(fIsInitialised) return; << 237 if(isInitialised) return; 217 fParticleChange = GetParticleChangeForGamma( 238 fParticleChange = GetParticleChangeForGamma(); 218 fIsInitialised = true; << 239 isInitialised = true; 219 } 240 } 220 241 221 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 222 243 223 void G4PenelopeRayleighModel::InitialiseLocal( 244 void G4PenelopeRayleighModel::InitialiseLocal(const G4ParticleDefinition* part, 224 G4VEmModel *masterModel) 245 G4VEmModel *masterModel) 225 { 246 { 226 if (fVerboseLevel > 3) << 247 if (verboseLevel > 3) 227 G4cout << "Calling G4PenelopeRayleighMode 248 G4cout << "Calling G4PenelopeRayleighModel::InitialiseLocal()" << G4endl; >> 249 228 // 250 // 229 //Check that particle matches: one might hav 251 //Check that particle matches: one might have multiple master models (e.g. 230 //for e+ and e-). 252 //for e+ and e-). 231 // 253 // 232 if (part == fParticle) 254 if (part == fParticle) 233 { 255 { 234 //Get the const table pointers from the 256 //Get the const table pointers from the master to the workers 235 const G4PenelopeRayleighModel* theModel 257 const G4PenelopeRayleighModel* theModel = 236 static_cast<G4PenelopeRayleighModel*> (maste 258 static_cast<G4PenelopeRayleighModel*> (masterModel); 237 259 238 //Copy pointers to the data tables 260 //Copy pointers to the data tables 239 for(G4int i=0; i<=fMaxZ; ++i) << 261 logAtomicCrossSection = theModel->logAtomicCrossSection; 240 { << 262 atomicFormFactor = theModel->atomicFormFactor; 241 fLogAtomicCrossSection[i] = theModel->fLog << 263 logFormFactorTable = theModel->logFormFactorTable; 242 fAtomicFormFactor[i] = theModel->fAtomicFo << 264 pMaxTable = theModel->pMaxTable; 243 } << 265 samplingTable = theModel->samplingTable; 244 fLogFormFactorTable = theModel->fLogForm << 245 fPMaxTable = theModel->fPMaxTable; << 246 fSamplingTable = theModel->fSamplingTabl << 247 266 248 //copy the G4DataVector with the grid 267 //copy the G4DataVector with the grid 249 fLogQSquareGrid = theModel->fLogQSquareG << 268 logQSquareGrid = theModel->logQSquareGrid; 250 269 251 //Same verbosity for all workers, as the 270 //Same verbosity for all workers, as the master 252 fVerboseLevel = theModel->fVerboseLevel; << 271 verboseLevel = theModel->verboseLevel; 253 } 272 } 254 273 255 return; 274 return; 256 } 275 } 257 276 258 277 259 //....oooOO0OOooo........oooOO0OOooo........oo 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 260 namespace { G4Mutex PenelopeRayleighModelMute 279 namespace { G4Mutex PenelopeRayleighModelMutex = G4MUTEX_INITIALIZER; } 261 G4double G4PenelopeRayleighModel::ComputeCross 280 G4double G4PenelopeRayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 262 G4double energy, 281 G4double energy, 263 G4double Z, 282 G4double Z, 264 G4double, 283 G4double, 265 G4double, 284 G4double, 266 G4double) 285 G4double) 267 { 286 { 268 // Cross section of Rayleigh scattering in P 287 // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97 269 // tabulation, Cuellen et al. (1997), with n 288 // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel 270 // et al. J. Phys. Chem. Ref. Data 4 (1975) 289 // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615. 271 290 272 if (fVerboseLevel > 3) << 291 if (verboseLevel > 3) 273 G4cout << "Calling CrossSectionPerAtom() o 292 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl; 274 293 275 G4int iZ = G4int(Z); << 294 G4int iZ = (G4int) Z; 276 295 277 if (!fLogAtomicCrossSection[iZ]) << 296 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was >> 297 //not invoked >> 298 if (!logAtomicCrossSection) >> 299 { >> 300 //create a **thread-local** version of the table. Used only for G4EmCalculator and >> 301 //Unit Tests >> 302 fLocalTable = true; >> 303 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>; >> 304 } >> 305 //now it should be ok >> 306 if (!logAtomicCrossSection->count(iZ)) 278 { 307 { 279 //If we are here, it means that Initial 308 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 280 //not filled up. This can happen in a U 309 //not filled up. This can happen in a UnitTest or via G4EmCalculator 281 if (fVerboseLevel > 0) << 310 if (verboseLevel > 0) 282 { 311 { 283 //Issue a G4Exception (warning) only in ve 312 //Issue a G4Exception (warning) only in verbose mode 284 G4ExceptionDescription ed; 313 G4ExceptionDescription ed; 285 ed << "Unable to retrieve the cross sectio 314 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl; 286 ed << "This can happen only in Unit Tests 315 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl; 287 G4Exception("G4PenelopeRayleighModel::Comp 316 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()", 288 "em2040",JustWarning,ed); 317 "em2040",JustWarning,ed); 289 } 318 } 290 //protect file reading via autolock 319 //protect file reading via autolock 291 G4AutoLock lock(&PenelopeRayleighModelM 320 G4AutoLock lock(&PenelopeRayleighModelMutex); 292 ReadDataFile(iZ); 321 ReadDataFile(iZ); 293 lock.unlock(); 322 lock.unlock(); 294 } 323 } 295 324 296 G4double cross = 0; 325 G4double cross = 0; 297 G4PhysicsFreeVector* atom = fLogAtomicCross << 326 >> 327 G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second; 298 if (!atom) 328 if (!atom) 299 { 329 { 300 G4ExceptionDescription ed; 330 G4ExceptionDescription ed; 301 ed << "Unable to find Z=" << iZ << " in 331 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl; 302 G4Exception("G4PenelopeRayleighModel::C 332 G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()", 303 "em2041",FatalException,ed); 333 "em2041",FatalException,ed); 304 return 0; 334 return 0; 305 } 335 } 306 G4double logene = G4Log(energy); << 336 G4double logene = std::log(energy); 307 G4double logXS = atom->Value(logene); 337 G4double logXS = atom->Value(logene); 308 cross = G4Exp(logXS); 338 cross = G4Exp(logXS); 309 339 310 if (fVerboseLevel > 2) << 340 if (verboseLevel > 2) 311 { 341 { 312 G4cout << "Rayleigh cross section at " 342 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z << 313 " = " << cross/barn << " barn" << G4endl; 343 " = " << cross/barn << " barn" << G4endl; 314 } 344 } >> 345 315 return cross; 346 return cross; 316 } 347 } 317 348 318 349 319 //....oooOO0OOooo........oooOO0OOooo........oo 350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 320 void G4PenelopeRayleighModel::BuildFormFactorT 351 void G4PenelopeRayleighModel::BuildFormFactorTable(const G4Material* material) 321 { 352 { 322 /* << 353 323 1) get composition and equivalent molecula << 354 /* 324 */ << 355 1) get composition and equivalent molecular density 325 std::size_t nElements = material->GetNumberO << 356 */ >> 357 >> 358 G4int nElements = material->GetNumberOfElements(); 326 const G4ElementVector* elementVector = mater 359 const G4ElementVector* elementVector = material->GetElementVector(); 327 const G4double* fractionVector = material->G 360 const G4double* fractionVector = material->GetFractionVector(); 328 361 329 std::vector<G4double> *StechiometricFactors 362 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>; 330 for (std::size_t i=0;i<nElements;++i) << 363 for (G4int i=0;i<nElements;i++) 331 { 364 { 332 G4double fraction = fractionVector[i]; 365 G4double fraction = fractionVector[i]; 333 G4double atomicWeigth = (*elementVector) 366 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole); 334 StechiometricFactors->push_back(fraction 367 StechiometricFactors->push_back(fraction/atomicWeigth); 335 } 368 } 336 //Find max 369 //Find max 337 G4double MaxStechiometricFactor = 0.; 370 G4double MaxStechiometricFactor = 0.; 338 for (std::size_t i=0;i<nElements;++i) << 371 for (G4int i=0;i<nElements;i++) 339 { 372 { 340 if ((*StechiometricFactors)[i] > MaxStec 373 if ((*StechiometricFactors)[i] > MaxStechiometricFactor) 341 MaxStechiometricFactor = (*Stechiometr 374 MaxStechiometricFactor = (*StechiometricFactors)[i]; 342 } 375 } 343 if (MaxStechiometricFactor<1e-16) 376 if (MaxStechiometricFactor<1e-16) 344 { 377 { 345 G4ExceptionDescription ed; 378 G4ExceptionDescription ed; 346 ed << "Inconsistent data of atomic compo 379 ed << "Inconsistent data of atomic composition for " << 347 material->GetName() << G4endl; 380 material->GetName() << G4endl; 348 G4Exception("G4PenelopeRayleighModel::Bu 381 G4Exception("G4PenelopeRayleighModel::BuildFormFactorTable()", 349 "em2042",FatalException,ed); 382 "em2042",FatalException,ed); 350 } 383 } 351 //Normalize 384 //Normalize 352 for (std::size_t i=0;i<nElements;++i) << 385 for (G4int i=0;i<nElements;i++) 353 (*StechiometricFactors)[i] /= MaxStechiom 386 (*StechiometricFactors)[i] /= MaxStechiometricFactor; 354 387 >> 388 // Equivalent atoms per molecule >> 389 G4double atomsPerMolecule = 0; >> 390 for (G4int i=0;i<nElements;i++) >> 391 atomsPerMolecule += (*StechiometricFactors)[i]; >> 392 355 /* 393 /* 356 CREATE THE FORM FACTOR TABLE 394 CREATE THE FORM FACTOR TABLE 357 */ 395 */ 358 G4PhysicsFreeVector* theFFVec = new G4Physic << 396 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size()); >> 397 theFFVec->SetSpline(true); 359 398 360 for (std::size_t k=0;k<fLogQSquareGrid.size( << 399 for (size_t k=0;k<logQSquareGrid.size();k++) 361 { 400 { 362 G4double ff2 = 0; //squared form factor 401 G4double ff2 = 0; //squared form factor 363 for (std::size_t i=0;i<nElements;++i) << 402 for (G4int i=0;i<nElements;i++) 364 { 403 { 365 G4int iZ = (*elementVector)[i]->GetZasInt( << 404 G4int iZ = (G4int) (*elementVector)[i]->GetZ(); 366 G4PhysicsFreeVector* theAtomVec = fAtomicF << 405 G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second; 367 G4double f = (*theAtomVec)[k]; //the q-gri 406 G4double f = (*theAtomVec)[k]; //the q-grid is always the same 368 ff2 += f*f*(*StechiometricFactors)[i]; 407 ff2 += f*f*(*StechiometricFactors)[i]; 369 } 408 } 370 if (ff2) 409 if (ff2) 371 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo << 410 theFFVec->PutValue(k,logQSquareGrid[k],std::log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2) 372 } 411 } 373 theFFVec->FillSecondDerivatives(); //vector << 412 logFormFactorTable->insert(std::make_pair(material,theFFVec)); 374 fLogFormFactorTable->insert(std::make_pair(m << 375 413 376 delete StechiometricFactors; 414 delete StechiometricFactors; >> 415 377 return; 416 return; 378 } 417 } 379 418 >> 419 380 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 381 421 382 void G4PenelopeRayleighModel::SampleSecondarie 422 void G4PenelopeRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* , 383 const G4MaterialCutsCouple* couple 423 const G4MaterialCutsCouple* couple, 384 const G4DynamicParticle* aDynamicG 424 const G4DynamicParticle* aDynamicGamma, 385 G4double, 425 G4double, 386 G4double) 426 G4double) 387 { 427 { 388 // Sampling of the Rayleigh final state (nam 428 // Sampling of the Rayleigh final state (namely, scattering angle of the photon) 389 // from the Penelope2008 model. The scatteri 429 // from the Penelope2008 model. The scattering angle is sampled from the atomic 390 // cross section dOmega/d(cosTheta) from Bor 430 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding 391 // anomalous scattering effects. The Form Fa 431 // anomalous scattering effects. The Form Factor F(Q) function which appears in the 392 // analytical cross section is retrieved via 432 // analytical cross section is retrieved via the method GetFSquared(); atomic data 393 // are tabulated for F(Q). Form factor for c 433 // are tabulated for F(Q). Form factor for compounds is calculated according to 394 // the additivity rule. The sampling from th 434 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse 395 // Transform with Aliasing (RITA) algorithm; 435 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once 396 // for each material and managed by G4Penelo 436 // for each material and managed by G4PenelopeSamplingData objects. 397 // The sampling algorithm (rejection method) 437 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and 398 // increases with energy. For E=100 keV the 438 // increases with energy. For E=100 keV the efficiency is 100% and 86% for 399 // hydrogen and uranium, respectively. 439 // hydrogen and uranium, respectively. 400 440 401 if (fVerboseLevel > 3) << 441 if (verboseLevel > 3) 402 G4cout << "Calling SamplingSecondaries() o 442 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModel" << G4endl; 403 443 404 G4double photonEnergy0 = aDynamicGamma->GetK 444 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 405 445 406 if (photonEnergy0 <= fIntrinsicLowEnergyLimi 446 if (photonEnergy0 <= fIntrinsicLowEnergyLimit) 407 { 447 { 408 fParticleChange->ProposeTrackStatus(fSto 448 fParticleChange->ProposeTrackStatus(fStopAndKill); 409 fParticleChange->SetProposedKineticEnerg 449 fParticleChange->SetProposedKineticEnergy(0.); 410 fParticleChange->ProposeLocalEnergyDepos 450 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); 411 return ; 451 return ; 412 } 452 } 413 453 414 G4ParticleMomentum photonDirection0 = aDynam 454 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 415 455 416 const G4Material* theMat = couple->GetMateri 456 const G4Material* theMat = couple->GetMaterial(); 417 457 >> 458 418 //1) Verify if tables are ready 459 //1) Verify if tables are ready 419 //Either Initialize() was not called, or we 460 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 420 //not invoked 461 //not invoked 421 if (!fPMaxTable || !fSamplingTable || !fLogF << 462 if (!pMaxTable || !samplingTable || !logAtomicCrossSection || !atomicFormFactor || >> 463 !logFormFactorTable) 422 { 464 { 423 //create a **thread-local** version of t 465 //create a **thread-local** version of the table. Used only for G4EmCalculator and 424 //Unit Tests 466 //Unit Tests 425 fLocalTable = true; 467 fLocalTable = true; 426 if (!fLogFormFactorTable) << 468 if (!logAtomicCrossSection) 427 fLogFormFactorTable = new std::map<const G4M << 469 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>; 428 if (!fPMaxTable) << 470 if (!atomicFormFactor) 429 fPMaxTable = new std::map<const G4Material*, << 471 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>; 430 if (!fSamplingTable) << 472 if (!logFormFactorTable) 431 fSamplingTable = new std::map<const G4Materi << 473 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; >> 474 if (!pMaxTable) >> 475 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; >> 476 if (!samplingTable) >> 477 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>; 432 } 478 } 433 479 434 if (!fSamplingTable->count(theMat)) << 480 if (!samplingTable->count(theMat)) 435 { 481 { 436 //If we are here, it means that Initiali 482 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 437 //not filled up. This can happen in a Un 483 //not filled up. This can happen in a UnitTest 438 if (fVerboseLevel > 0) << 484 if (verboseLevel > 0) 439 { 485 { 440 //Issue a G4Exception (warning) only in ve 486 //Issue a G4Exception (warning) only in verbose mode 441 G4ExceptionDescription ed; 487 G4ExceptionDescription ed; 442 ed << "Unable to find the fSamplingTable d << 488 ed << "Unable to find the samplingTable data for " << 443 theMat->GetName() << G4endl; 489 theMat->GetName() << G4endl; 444 ed << "This can happen only in Unit Tests" 490 ed << "This can happen only in Unit Tests" << G4endl; 445 G4Exception("G4PenelopeRayleighModel::Samp 491 G4Exception("G4PenelopeRayleighModel::SampleSecondaries()", 446 "em2019",JustWarning,ed); 492 "em2019",JustWarning,ed); 447 } 493 } 448 const G4ElementVector* theElementVector 494 const G4ElementVector* theElementVector = theMat->GetElementVector(); 449 //protect file reading via autolock 495 //protect file reading via autolock 450 G4AutoLock lock(&PenelopeRayleighModelMu 496 G4AutoLock lock(&PenelopeRayleighModelMutex); 451 for (std::size_t j=0;j<theMat->GetNumber << 497 for (size_t j=0;j<theMat->GetNumberOfElements();j++) 452 { 498 { 453 G4int iZ = theElementVector->at(j)->GetZas << 499 G4int iZ = (G4int) theElementVector->at(j)->GetZ(); 454 if (!fLogAtomicCrossSection[iZ]) << 500 if (!logAtomicCrossSection->count(iZ)) 455 { 501 { 456 lock.lock(); 502 lock.lock(); 457 ReadDataFile(iZ); 503 ReadDataFile(iZ); 458 lock.unlock(); 504 lock.unlock(); 459 } 505 } 460 } 506 } 461 lock.lock(); 507 lock.lock(); 462 //1) If the table has not been built for 508 //1) If the table has not been built for the material, do it! 463 if (!fLogFormFactorTable->count(theMat)) << 509 if (!logFormFactorTable->count(theMat)) 464 BuildFormFactorTable(theMat); 510 BuildFormFactorTable(theMat); 465 511 466 //2) retrieve or build the sampling tabl 512 //2) retrieve or build the sampling table 467 if (!(fSamplingTable->count(theMat))) << 513 if (!(samplingTable->count(theMat))) 468 InitializeSamplingAlgorithm(theMat); 514 InitializeSamplingAlgorithm(theMat); 469 515 470 //3) retrieve or build the pMax data 516 //3) retrieve or build the pMax data 471 if (!fPMaxTable->count(theMat)) << 517 if (!pMaxTable->count(theMat)) 472 GetPMaxTable(theMat); 518 GetPMaxTable(theMat); 473 lock.unlock(); 519 lock.unlock(); 474 } 520 } 475 521 476 //Ok, restart the job 522 //Ok, restart the job 477 G4PenelopeSamplingData* theDataTable = fSamp << 523 478 G4PhysicsFreeVector* thePMax = fPMaxTable->f << 524 G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second; >> 525 G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second; 479 526 480 G4double cosTheta = 1.0; 527 G4double cosTheta = 1.0; 481 528 482 //OK, ready to go! 529 //OK, ready to go! 483 G4double qmax = 2.0*photonEnergy0/electron_m 530 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now 484 531 485 if (qmax < 1e-10) //very low momentum transf 532 if (qmax < 1e-10) //very low momentum transfer 486 { 533 { 487 G4bool loopAgain=false; 534 G4bool loopAgain=false; 488 do 535 do 489 { 536 { 490 loopAgain = false; 537 loopAgain = false; 491 cosTheta = 1.0-2.0*G4UniformRand(); 538 cosTheta = 1.0-2.0*G4UniformRand(); 492 G4double G = 0.5*(1+cosTheta*cosTheta); 539 G4double G = 0.5*(1+cosTheta*cosTheta); 493 if (G4UniformRand()>G) 540 if (G4UniformRand()>G) 494 loopAgain = true; 541 loopAgain = true; 495 }while(loopAgain); 542 }while(loopAgain); 496 } 543 } 497 else //larger momentum transfer 544 else //larger momentum transfer 498 { 545 { 499 std::size_t nData = theDataTable->GetNum << 546 size_t nData = theDataTable->GetNumberOfStoredPoints(); 500 G4double LastQ2inTheTable = theDataTable 547 G4double LastQ2inTheTable = theDataTable->GetX(nData-1); 501 G4double q2max = std::min(qmax*qmax,Last 548 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable); 502 549 503 G4bool loopAgain = false; 550 G4bool loopAgain = false; 504 G4double MaxPValue = thePMax->Value(phot 551 G4double MaxPValue = thePMax->Value(photonEnergy0); 505 G4double xx=0; 552 G4double xx=0; 506 553 507 //Sampling by rejection method. The reje 554 //Sampling by rejection method. The rejection function is 508 //G = 0.5*(1+cos^2(theta)) 555 //G = 0.5*(1+cos^2(theta)) 509 // 556 // 510 do{ 557 do{ 511 loopAgain = false; 558 loopAgain = false; 512 G4double RandomMax = G4UniformRand()*MaxPVal 559 G4double RandomMax = G4UniformRand()*MaxPValue; 513 xx = theDataTable->SampleValue(RandomMax); 560 xx = theDataTable->SampleValue(RandomMax); 514 //xx is a random value of q^2 in (0,q2max),s 561 //xx is a random value of q^2 in (0,q2max),sampled according to 515 //F(Q^2) via the RITA algorithm 562 //F(Q^2) via the RITA algorithm 516 if (xx > q2max) 563 if (xx > q2max) 517 loopAgain = true; 564 loopAgain = true; 518 cosTheta = 1.0-2.0*xx/q2max; 565 cosTheta = 1.0-2.0*xx/q2max; 519 G4double G = 0.5*(1+cosTheta*cosTheta); 566 G4double G = 0.5*(1+cosTheta*cosTheta); 520 if (G4UniformRand()>G) 567 if (G4UniformRand()>G) 521 loopAgain = true; 568 loopAgain = true; 522 }while(loopAgain); 569 }while(loopAgain); 523 } 570 } 524 571 525 G4double sinTheta = std::sqrt(1-cosTheta*cos 572 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); 526 573 527 // Scattered photon angles. ( Z - axis along 574 // Scattered photon angles. ( Z - axis along the parent photon) 528 G4double phi = twopi * G4UniformRand() ; 575 G4double phi = twopi * G4UniformRand() ; 529 G4double dirX = sinTheta*std::cos(phi); 576 G4double dirX = sinTheta*std::cos(phi); 530 G4double dirY = sinTheta*std::sin(phi); 577 G4double dirY = sinTheta*std::sin(phi); 531 G4double dirZ = cosTheta; 578 G4double dirZ = cosTheta; 532 579 533 // Update G4VParticleChange for the scattere 580 // Update G4VParticleChange for the scattered photon 534 G4ThreeVector photonDirection1(dirX, dirY, d 581 G4ThreeVector photonDirection1(dirX, dirY, dirZ); 535 582 536 photonDirection1.rotateUz(photonDirection0); 583 photonDirection1.rotateUz(photonDirection0); 537 fParticleChange->ProposeMomentumDirection(ph 584 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 538 fParticleChange->SetProposedKineticEnergy(ph 585 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ; 539 586 540 return; 587 return; 541 } 588 } 542 589 543 590 544 //....oooOO0OOooo........oooOO0OOooo........oo 591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 545 592 546 void G4PenelopeRayleighModel::ReadDataFile(con 593 void G4PenelopeRayleighModel::ReadDataFile(const G4int Z) 547 { 594 { 548 if (fVerboseLevel > 2) << 595 >> 596 if (verboseLevel > 2) 549 { 597 { 550 G4cout << "G4PenelopeRayleighModel::Read 598 G4cout << "G4PenelopeRayleighModel::ReadDataFile()" << G4endl; 551 G4cout << "Going to read Rayleigh data f 599 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl; 552 } 600 } 553 const char* path = G4FindDataDir("G4LEDATA << 601 554 if(!path) << 602 char* path = getenv("G4LEDATA"); >> 603 if (!path) 555 { 604 { 556 G4String excep = "G4LEDATA environment v 605 G4String excep = "G4LEDATA environment variable not set!"; 557 G4Exception("G4PenelopeRayleighModel::Re 606 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 558 "em0006",FatalException,excep); 607 "em0006",FatalException,excep); 559 return; 608 return; 560 } 609 } 561 610 562 /* 611 /* 563 Read first the cross section file 612 Read first the cross section file 564 */ 613 */ 565 std::ostringstream ost; 614 std::ostringstream ost; 566 if (Z>9) 615 if (Z>9) 567 ost << path << "/penelope/rayleigh/pdgra" 616 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08"; 568 else 617 else 569 ost << path << "/penelope/rayleigh/pdgra0" 618 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08"; 570 std::ifstream file(ost.str().c_str()); 619 std::ifstream file(ost.str().c_str()); 571 if (!file.is_open()) 620 if (!file.is_open()) 572 { 621 { 573 G4String excep = "Data file " + G4String 622 G4String excep = "Data file " + G4String(ost.str()) + " not found!"; 574 G4Exception("G4PenelopeRayleighModel::Re 623 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 575 "em0003",FatalException,excep); 624 "em0003",FatalException,excep); 576 } 625 } 577 G4int readZ =0; 626 G4int readZ =0; 578 std::size_t nPoints= 0; << 627 size_t nPoints= 0; 579 file >> readZ >> nPoints; 628 file >> readZ >> nPoints; 580 //check the right file is opened. 629 //check the right file is opened. 581 if (readZ != Z || nPoints <= 0 || nPoints >= 630 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) 582 { 631 { 583 G4ExceptionDescription ed; 632 G4ExceptionDescription ed; 584 ed << "Corrupted data file for Z=" << Z 633 ed << "Corrupted data file for Z=" << Z << G4endl; 585 G4Exception("G4PenelopeRayleighModel::Re 634 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 586 "em0005",FatalException,ed); 635 "em0005",FatalException,ed); 587 return; 636 return; 588 } 637 } 589 << 638 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints); 590 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 591 G4double ene=0,f1=0,f2=0,xs=0; 639 G4double ene=0,f1=0,f2=0,xs=0; 592 for (std::size_t i=0;i<nPoints;++i) << 640 for (size_t i=0;i<nPoints;i++) 593 { 641 { 594 file >> ene >> f1 >> f2 >> xs; 642 file >> ene >> f1 >> f2 >> xs; 595 //dimensional quantities 643 //dimensional quantities 596 ene *= eV; 644 ene *= eV; 597 xs *= cm2; 645 xs *= cm2; 598 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 646 theVec->PutValue(i,std::log(ene),std::log(xs)); 599 if (file.eof() && i != (nPoints-1)) //fi 647 if (file.eof() && i != (nPoints-1)) //file ended too early 600 { 648 { 601 G4ExceptionDescription ed ; 649 G4ExceptionDescription ed ; 602 ed << "Corrupted data file for Z=" << Z << 650 ed << "Corrupted data file for Z=" << Z << G4endl; 603 ed << "Found less than " << nPoints << "en 651 ed << "Found less than " << nPoints << "entries " <<G4endl; 604 G4Exception("G4PenelopeRayleighModel::Read 652 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 605 "em0005",FatalException,ed); 653 "em0005",FatalException,ed); 606 } 654 } 607 } 655 } >> 656 if (!logAtomicCrossSection) >> 657 { >> 658 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", >> 659 "em2044",FatalException,"Unable to allocate the atomic cross section table"); >> 660 delete theVec; >> 661 return; >> 662 } >> 663 logAtomicCrossSection->insert(std::make_pair(Z,theVec)); 608 file.close(); 664 file.close(); 609 665 610 /* 666 /* 611 Then read the form factor file 667 Then read the form factor file 612 */ 668 */ 613 std::ostringstream ost2; 669 std::ostringstream ost2; 614 if (Z>9) 670 if (Z>9) 615 ost2 << path << "/penelope/rayleigh/pdaff" 671 ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08"; 616 else 672 else 617 ost2 << path << "/penelope/rayleigh/pdaff0 673 ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08"; 618 file.open(ost2.str().c_str()); 674 file.open(ost2.str().c_str()); 619 if (!file.is_open()) 675 if (!file.is_open()) 620 { 676 { 621 G4String excep = "Data file " + G4String 677 G4String excep = "Data file " + G4String(ost2.str()) + " not found!"; 622 G4Exception("G4PenelopeRayleighModel::Re 678 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 623 "em0003",FatalException,excep); 679 "em0003",FatalException,excep); 624 } 680 } 625 file >> readZ >> nPoints; 681 file >> readZ >> nPoints; 626 //check the right file is opened. 682 //check the right file is opened. 627 if (readZ != Z || nPoints <= 0 || nPoints >= 683 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) 628 { 684 { 629 G4ExceptionDescription ed; 685 G4ExceptionDescription ed; 630 ed << "Corrupted data file for Z=" << Z 686 ed << "Corrupted data file for Z=" << Z << G4endl; 631 G4Exception("G4PenelopeRayleighModel::Re 687 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 632 "em0005",FatalException,ed); 688 "em0005",FatalException,ed); 633 return; 689 return; 634 } 690 } 635 fAtomicFormFactor[Z] = new G4PhysicsFreeVect << 691 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints); 636 G4double q=0,ff=0,incoh=0; 692 G4double q=0,ff=0,incoh=0; 637 G4bool fillQGrid = false; 693 G4bool fillQGrid = false; 638 //fill this vector only the first time. 694 //fill this vector only the first time. 639 if (!fLogQSquareGrid.size()) << 695 if (!logQSquareGrid.size()) 640 fillQGrid = true; 696 fillQGrid = true; 641 for (std::size_t i=0;i<nPoints;++i) << 697 for (size_t i=0;i<nPoints;i++) 642 { 698 { 643 file >> q >> ff >> incoh; 699 file >> q >> ff >> incoh; 644 //q and ff are dimensionless (q is in un 700 //q and ff are dimensionless (q is in units of (m_e*c) 645 fAtomicFormFactor[Z]->PutValue(i,q,ff); << 701 theFFVec->PutValue(i,q,ff); 646 if (fillQGrid) 702 if (fillQGrid) 647 { 703 { 648 fLogQSquareGrid.push_back(2.0*G4Log(q)); << 704 logQSquareGrid.push_back(2.0*std::log(q)); 649 } 705 } 650 if (file.eof() && i != (nPoints-1)) //fi 706 if (file.eof() && i != (nPoints-1)) //file ended too early 651 { 707 { 652 G4ExceptionDescription ed; 708 G4ExceptionDescription ed; 653 ed << "Corrupted data file for Z=" << Z << 709 ed << "Corrupted data file for Z=" << Z << G4endl; 654 ed << "Found less than " << nPoints << "en 710 ed << "Found less than " << nPoints << "entries " <<G4endl; 655 G4Exception("G4PenelopeRayleighModel::Read 711 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", 656 "em0005",FatalException,ed); 712 "em0005",FatalException,ed); 657 } 713 } 658 } 714 } >> 715 if (!atomicFormFactor) >> 716 { >> 717 G4Exception("G4PenelopeRayleighModel::ReadDataFile()", >> 718 "em2045",FatalException, >> 719 "Unable to allocate the atomicFormFactor data table"); >> 720 delete theFFVec; >> 721 return; >> 722 } >> 723 atomicFormFactor->insert(std::make_pair(Z,theFFVec)); 659 file.close(); 724 file.close(); 660 return; 725 return; 661 } 726 } 662 727 663 //....oooOO0OOooo........oooOO0OOooo........oo 728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 664 729 665 G4double G4PenelopeRayleighModel::GetFSquared( 730 G4double G4PenelopeRayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared) 666 { 731 { 667 G4double f2 = 0; 732 G4double f2 = 0; 668 //Input value QSquared could be zero: protec 733 //Input value QSquared could be zero: protect the log() below against 669 //the FPE exception 734 //the FPE exception 670 //If Q<1e-10, set Q to 1e-10 735 //If Q<1e-10, set Q to 1e-10 671 G4double logQSquared = (QSquared>1e-10) ? G4 << 736 G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.; 672 //last value of the table 737 //last value of the table 673 G4double maxlogQ2 = fLogQSquareGrid[fLogQSqu << 738 G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1]; >> 739 674 740 675 //now it should be all right 741 //now it should be all right 676 G4PhysicsFreeVector* theVec = fLogFormFactor << 742 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; 677 743 678 if (!theVec) 744 if (!theVec) 679 { 745 { 680 G4ExceptionDescription ed; 746 G4ExceptionDescription ed; 681 ed << "Unable to retrieve F squared tabl 747 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl; 682 G4Exception("G4PenelopeRayleighModel::Ge 748 G4Exception("G4PenelopeRayleighModel::GetFSquared()", 683 "em2046",FatalException,ed); 749 "em2046",FatalException,ed); 684 return 0; 750 return 0; 685 } 751 } 686 if (logQSquared < -20) // Q < 1e-9 752 if (logQSquared < -20) // Q < 1e-9 687 { 753 { 688 G4double logf2 = (*theVec)[0]; //first v 754 G4double logf2 = (*theVec)[0]; //first value of the table 689 f2 = G4Exp(logf2); 755 f2 = G4Exp(logf2); 690 } 756 } 691 else if (logQSquared > maxlogQ2) 757 else if (logQSquared > maxlogQ2) 692 f2 =0; 758 f2 =0; 693 else 759 else 694 { 760 { 695 //log(Q^2) vs. log(F^2) 761 //log(Q^2) vs. log(F^2) 696 G4double logf2 = theVec->Value(logQSquar 762 G4double logf2 = theVec->Value(logQSquared); 697 f2 = G4Exp(logf2); 763 f2 = G4Exp(logf2); 698 764 699 } 765 } 700 if (fVerboseLevel > 3) << 766 if (verboseLevel > 3) 701 { 767 { 702 G4cout << "G4PenelopeRayleighModel::GetF 768 G4cout << "G4PenelopeRayleighModel::GetFSquared() in " << mat->GetName() << G4endl; 703 G4cout << "Q^2 = " << QSquared << " (un 769 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl; 704 } 770 } 705 return f2; 771 return f2; 706 } 772 } 707 773 708 //....oooOO0OOooo........oooOO0OOooo........oo 774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 709 775 710 void G4PenelopeRayleighModel::InitializeSampli 776 void G4PenelopeRayleighModel::InitializeSamplingAlgorithm(const G4Material* mat) 711 { 777 { >> 778 712 G4double q2min = 0; 779 G4double q2min = 0; 713 G4double q2max = 0; 780 G4double q2max = 0; 714 const std::size_t np = 150; //hard-coded in << 781 const size_t np = 150; //hard-coded in Penelope 715 //G4cout << "Init N= " << fLogQSquareGrid.si << 782 //G4cout << "Init N= " << logQSquareGrid.size() << G4endl; 716 for (std::size_t i=1;i<fLogQSquareGrid.size( << 783 for (size_t i=1;i<logQSquareGrid.size();i++) 717 { 784 { 718 G4double Q2 = G4Exp(fLogQSquareGrid[i]); << 785 G4double Q2 = G4Exp(logQSquareGrid[i]); 719 if (GetFSquared(mat,Q2) > 1e-35) 786 if (GetFSquared(mat,Q2) > 1e-35) 720 { 787 { 721 q2max = G4Exp(fLogQSquareGrid[i-1]); << 788 q2max = G4Exp(logQSquareGrid[i-1]); 722 } 789 } 723 //G4cout << "Q2= " << Q2 << " q2max= " < 790 //G4cout << "Q2= " << Q2 << " q2max= " << q2max << G4endl; 724 } 791 } 725 792 726 std::size_t nReducedPoints = np/4; << 793 size_t nReducedPoints = np/4; 727 794 728 //check for errors 795 //check for errors 729 if (np < 16) 796 if (np < 16) 730 { 797 { 731 G4Exception("G4PenelopeRayleighModel::In 798 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()", 732 "em2047",FatalException, 799 "em2047",FatalException, 733 "Too few points to initialize the sampli 800 "Too few points to initialize the sampling algorithm"); 734 } 801 } 735 if (q2min > (q2max-1e-10)) 802 if (q2min > (q2max-1e-10)) 736 { 803 { 737 G4cout << "q2min= " << q2min << " q2max= 804 G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl; 738 G4Exception("G4PenelopeRayleighModel::In 805 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()", 739 "em2048",FatalException, 806 "em2048",FatalException, 740 "Too narrow grid to initialize the sampl 807 "Too narrow grid to initialize the sampling algorithm"); 741 } 808 } 742 809 743 //This is subroutine RITAI0 of Penelope 810 //This is subroutine RITAI0 of Penelope 744 //Create an object of type G4PenelopeRayleig 811 //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material* 745 812 746 //temporary vectors --> Then everything is s 813 //temporary vectors --> Then everything is stored in G4PenelopeSamplingData 747 G4DataVector* x = new G4DataVector(); 814 G4DataVector* x = new G4DataVector(); 748 815 749 /******************************************* 816 /******************************************************************************* 750 Start with a grid of NUNIF points uniforml 817 Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max 751 ******************************************** 818 ********************************************************************************/ 752 std::size_t NUNIF = std::min(std::max(((std: << 819 size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2); 753 const G4int nip = 51; //hard-coded in Penelo 820 const G4int nip = 51; //hard-coded in Penelope 754 821 755 G4double dx = (q2max-q2min)/((G4double) NUNI 822 G4double dx = (q2max-q2min)/((G4double) NUNIF-1); 756 x->push_back(q2min); 823 x->push_back(q2min); 757 for (std::size_t i=1;i<NUNIF-1;++i) << 824 for (size_t i=1;i<NUNIF-1;i++) 758 { 825 { 759 G4double app = q2min + i*dx; 826 G4double app = q2min + i*dx; 760 x->push_back(app); //increase 827 x->push_back(app); //increase 761 } 828 } 762 x->push_back(q2max); 829 x->push_back(q2max); 763 830 764 if (fVerboseLevel> 3) << 831 if (verboseLevel> 3) 765 G4cout << "Vector x has " << x->size() << 832 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl; 766 833 767 //Allocate temporary storage vectors 834 //Allocate temporary storage vectors 768 G4DataVector* area = new G4DataVector(); 835 G4DataVector* area = new G4DataVector(); 769 G4DataVector* a = new G4DataVector(); 836 G4DataVector* a = new G4DataVector(); 770 G4DataVector* b = new G4DataVector(); 837 G4DataVector* b = new G4DataVector(); 771 G4DataVector* c = new G4DataVector(); 838 G4DataVector* c = new G4DataVector(); 772 G4DataVector* err = new G4DataVector(); 839 G4DataVector* err = new G4DataVector(); 773 840 774 for (std::size_t i=0;i<NUNIF-1;++i) //build << 841 for (size_t i=0;i<NUNIF-1;i++) //build all points but the last 775 { 842 { 776 //Temporary vectors for this loop 843 //Temporary vectors for this loop 777 G4DataVector* pdfi = new G4DataVector(); 844 G4DataVector* pdfi = new G4DataVector(); 778 G4DataVector* pdfih = new G4DataVector() 845 G4DataVector* pdfih = new G4DataVector(); 779 G4DataVector* sumi = new G4DataVector(); 846 G4DataVector* sumi = new G4DataVector(); 780 847 781 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4do 848 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); 782 G4double pdfmax = 0; 849 G4double pdfmax = 0; 783 for (G4int k=0;k<nip;k++) 850 for (G4int k=0;k<nip;k++) 784 { 851 { 785 G4double xik = (*x)[i]+k*dxi; 852 G4double xik = (*x)[i]+k*dxi; 786 G4double pdfk = std::max(GetFSquared(mat,x 853 G4double pdfk = std::max(GetFSquared(mat,xik),0.); 787 pdfi->push_back(pdfk); 854 pdfi->push_back(pdfk); 788 pdfmax = std::max(pdfmax,pdfk); 855 pdfmax = std::max(pdfmax,pdfk); 789 if (k < (nip-1)) 856 if (k < (nip-1)) 790 { 857 { 791 G4double xih = xik + 0.5*dxi; 858 G4double xih = xik + 0.5*dxi; 792 G4double pdfIK = std::max(GetFSquared( 859 G4double pdfIK = std::max(GetFSquared(mat,xih),0.); 793 pdfih->push_back(pdfIK); 860 pdfih->push_back(pdfIK); 794 pdfmax = std::max(pdfmax,pdfIK); 861 pdfmax = std::max(pdfmax,pdfIK); 795 } 862 } 796 } 863 } 797 864 798 //Simpson's integration 865 //Simpson's integration 799 G4double cons = dxi*0.5*(1./3.); 866 G4double cons = dxi*0.5*(1./3.); 800 sumi->push_back(0.); 867 sumi->push_back(0.); 801 for (G4int k=1;k<nip;k++) 868 for (G4int k=1;k<nip;k++) 802 { 869 { 803 G4double previous = (*sumi)[k-1]; 870 G4double previous = (*sumi)[k-1]; 804 G4double next = previous + cons*((*pdfi)[k 871 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); 805 sumi->push_back(next); 872 sumi->push_back(next); 806 } 873 } 807 874 808 G4double lastIntegral = (*sumi)[sumi->si 875 G4double lastIntegral = (*sumi)[sumi->size()-1]; 809 area->push_back(lastIntegral); 876 area->push_back(lastIntegral); 810 //Normalize cumulative function 877 //Normalize cumulative function 811 G4double factor = 1.0/lastIntegral; 878 G4double factor = 1.0/lastIntegral; 812 for (std::size_t k=0;k<sumi->size();++k) << 879 for (size_t k=0;k<sumi->size();k++) 813 (*sumi)[k] *= factor; 880 (*sumi)[k] *= factor; 814 881 815 //When the PDF vanishes at one of the in 882 //When the PDF vanishes at one of the interval end points, its value is modified 816 if ((*pdfi)[0] < 1e-35) 883 if ((*pdfi)[0] < 1e-35) 817 (*pdfi)[0] = 1e-5*pdfmax; 884 (*pdfi)[0] = 1e-5*pdfmax; 818 if ((*pdfi)[pdfi->size()-1] < 1e-35) 885 if ((*pdfi)[pdfi->size()-1] < 1e-35) 819 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 886 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 820 887 821 G4double pli = (*pdfi)[0]*factor; 888 G4double pli = (*pdfi)[0]*factor; 822 G4double pui = (*pdfi)[pdfi->size()-1]*f 889 G4double pui = (*pdfi)[pdfi->size()-1]*factor; 823 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx 890 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx); 824 G4double A_temp = (1.0/(pli*dx))-1.0-B_t 891 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp; 825 G4double C_temp = 1.0+A_temp+B_temp; 892 G4double C_temp = 1.0+A_temp+B_temp; 826 if (C_temp < 1e-35) 893 if (C_temp < 1e-35) 827 { 894 { 828 a->push_back(0.); 895 a->push_back(0.); 829 b->push_back(0.); 896 b->push_back(0.); 830 c->push_back(1.); 897 c->push_back(1.); 831 } 898 } 832 else 899 else 833 { 900 { 834 a->push_back(A_temp); 901 a->push_back(A_temp); 835 b->push_back(B_temp); 902 b->push_back(B_temp); 836 c->push_back(C_temp); 903 c->push_back(C_temp); 837 } 904 } 838 905 839 //OK, now get ERR(I), the integral of th 906 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 840 //and the true pdf, extended over the in 907 //and the true pdf, extended over the interval (X(I),X(I+1)) 841 G4int icase = 1; //loop code 908 G4int icase = 1; //loop code 842 G4bool reLoop = false; 909 G4bool reLoop = false; 843 err->push_back(0.); 910 err->push_back(0.); 844 do 911 do 845 { 912 { 846 reLoop = false; 913 reLoop = false; 847 (*err)[i] = 0.; //zero variable 914 (*err)[i] = 0.; //zero variable 848 for (G4int k=0;k<nip;k++) 915 for (G4int k=0;k<nip;k++) 849 { 916 { 850 G4double rr = (*sumi)[k]; 917 G4double rr = (*sumi)[k]; 851 G4double pap = (*area)[i]*(1.0+((*a)[i 918 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ 852 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(* 919 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); 853 if (k == 0 || k == nip-1) 920 if (k == 0 || k == nip-1) 854 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]) 921 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); 855 else 922 else 856 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 923 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 857 } 924 } 858 (*err)[i] *= dxi; 925 (*err)[i] *= dxi; 859 926 860 //If err(I) is too large, the pdf is appro 927 //If err(I) is too large, the pdf is approximated by a uniform distribution 861 if ((*err)[i] > 0.1*(*area)[i] && icase == 928 if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 862 { 929 { 863 (*b)[i] = 0; 930 (*b)[i] = 0; 864 (*a)[i] = 0; 931 (*a)[i] = 0; 865 (*c)[i] = 1.; 932 (*c)[i] = 1.; 866 icase = 2; 933 icase = 2; 867 reLoop = true; 934 reLoop = true; 868 } 935 } 869 }while(reLoop); 936 }while(reLoop); >> 937 870 delete pdfi; 938 delete pdfi; 871 delete pdfih; 939 delete pdfih; 872 delete sumi; 940 delete sumi; 873 } //end of first loop over i 941 } //end of first loop over i 874 942 875 //Now assign last point 943 //Now assign last point 876 (*x)[x->size()-1] = q2max; 944 (*x)[x->size()-1] = q2max; 877 a->push_back(0.); 945 a->push_back(0.); 878 b->push_back(0.); 946 b->push_back(0.); 879 c->push_back(0.); 947 c->push_back(0.); 880 err->push_back(0.); 948 err->push_back(0.); 881 area->push_back(0.); 949 area->push_back(0.); 882 950 883 if (x->size() != NUNIF || a->size() != NUNIF 951 if (x->size() != NUNIF || a->size() != NUNIF || 884 err->size() != NUNIF || area->size() != 952 err->size() != NUNIF || area->size() != NUNIF) 885 { 953 { 886 G4ExceptionDescription ed; 954 G4ExceptionDescription ed; 887 ed << "Problem in building the Table for 955 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl; 888 G4Exception("G4PenelopeRayleighModel::In 956 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()", 889 "em2049",FatalException,ed); 957 "em2049",FatalException,ed); 890 } 958 } 891 959 892 /******************************************* 960 /******************************************************************************* 893 New grid points are added by halving the su 961 New grid points are added by halving the sub-intervals with the largest absolute error 894 This is done up to np=150 points in the grid 962 This is done up to np=150 points in the grid 895 ******************************************** 963 ********************************************************************************/ 896 do 964 do 897 { 965 { 898 G4double maxError = 0.0; 966 G4double maxError = 0.0; 899 std::size_t iErrMax = 0; << 967 size_t iErrMax = 0; 900 for (std::size_t i=0;i<err->size()-2;++i << 968 for (size_t i=0;i<err->size()-2;i++) 901 { 969 { 902 //maxError is the lagest of the interval e 970 //maxError is the lagest of the interval errors err[i] 903 if ((*err)[i] > maxError) 971 if ((*err)[i] > maxError) 904 { 972 { 905 maxError = (*err)[i]; 973 maxError = (*err)[i]; 906 iErrMax = i; 974 iErrMax = i; 907 } 975 } 908 } 976 } 909 977 910 //OK, now I have to insert one new point 978 //OK, now I have to insert one new point in the position iErrMax 911 G4double newx = 0.5*((*x)[iErrMax]+(*x)[ 979 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]); 912 980 913 x->insert(x->begin()+iErrMax+1,newx); 981 x->insert(x->begin()+iErrMax+1,newx); 914 //Add place-holders in the other vectors 982 //Add place-holders in the other vectors 915 area->insert(area->begin()+iErrMax+1,0.) 983 area->insert(area->begin()+iErrMax+1,0.); 916 a->insert(a->begin()+iErrMax+1,0.); 984 a->insert(a->begin()+iErrMax+1,0.); 917 b->insert(b->begin()+iErrMax+1,0.); 985 b->insert(b->begin()+iErrMax+1,0.); 918 c->insert(c->begin()+iErrMax+1,0.); 986 c->insert(c->begin()+iErrMax+1,0.); 919 err->insert(err->begin()+iErrMax+1,0.); 987 err->insert(err->begin()+iErrMax+1,0.); 920 988 921 //Now calculate the other parameters 989 //Now calculate the other parameters 922 for (std::size_t i=iErrMax;i<=iErrMax+1; << 990 for (size_t i=iErrMax;i<=iErrMax+1;i++) 923 { 991 { 924 //Temporary vectors for this loop 992 //Temporary vectors for this loop 925 G4DataVector* pdfi = new G4DataVector(); 993 G4DataVector* pdfi = new G4DataVector(); 926 G4DataVector* pdfih = new G4DataVector(); 994 G4DataVector* pdfih = new G4DataVector(); 927 G4DataVector* sumi = new G4DataVector(); 995 G4DataVector* sumi = new G4DataVector(); 928 996 929 G4double dxLocal = (*x)[i+1]-(*x)[i]; 997 G4double dxLocal = (*x)[i+1]-(*x)[i]; 930 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4doub 998 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); 931 G4double pdfmax = 0; 999 G4double pdfmax = 0; 932 for (G4int k=0;k<nip;k++) 1000 for (G4int k=0;k<nip;k++) 933 { 1001 { 934 G4double xik = (*x)[i]+k*dxi; 1002 G4double xik = (*x)[i]+k*dxi; 935 G4double pdfk = std::max(GetFSquared(m 1003 G4double pdfk = std::max(GetFSquared(mat,xik),0.); 936 pdfi->push_back(pdfk); 1004 pdfi->push_back(pdfk); 937 pdfmax = std::max(pdfmax,pdfk); 1005 pdfmax = std::max(pdfmax,pdfk); 938 if (k < (nip-1)) 1006 if (k < (nip-1)) 939 { 1007 { 940 G4double xih = xik + 0.5*dxi; 1008 G4double xih = xik + 0.5*dxi; 941 G4double pdfIK = std::max(GetFSquared(ma 1009 G4double pdfIK = std::max(GetFSquared(mat,xih),0.); 942 pdfih->push_back(pdfIK); 1010 pdfih->push_back(pdfIK); 943 pdfmax = std::max(pdfmax,pdfIK); 1011 pdfmax = std::max(pdfmax,pdfIK); 944 } 1012 } 945 } 1013 } 946 1014 947 //Simpson's integration 1015 //Simpson's integration 948 G4double cons = dxi*0.5*(1./3.); 1016 G4double cons = dxi*0.5*(1./3.); 949 sumi->push_back(0.); 1017 sumi->push_back(0.); 950 for (G4int k=1;k<nip;k++) 1018 for (G4int k=1;k<nip;k++) 951 { 1019 { 952 G4double previous = (*sumi)[k-1]; 1020 G4double previous = (*sumi)[k-1]; 953 G4double next = previous + cons*((*pdf 1021 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); 954 sumi->push_back(next); 1022 sumi->push_back(next); 955 } 1023 } 956 G4double lastIntegral = (*sumi)[sumi->size 1024 G4double lastIntegral = (*sumi)[sumi->size()-1]; 957 (*area)[i] = lastIntegral; 1025 (*area)[i] = lastIntegral; 958 1026 959 //Normalize cumulative function 1027 //Normalize cumulative function 960 G4double factor = 1.0/lastIntegral; 1028 G4double factor = 1.0/lastIntegral; 961 for (std::size_t k=0;k<sumi->size();++k) << 1029 for (size_t k=0;k<sumi->size();k++) 962 (*sumi)[k] *= factor; 1030 (*sumi)[k] *= factor; 963 1031 964 //When the PDF vanishes at one of the inte 1032 //When the PDF vanishes at one of the interval end points, its value is modified 965 if ((*pdfi)[0] < 1e-35) 1033 if ((*pdfi)[0] < 1e-35) 966 (*pdfi)[0] = 1e-5*pdfmax; 1034 (*pdfi)[0] = 1e-5*pdfmax; 967 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1035 if ((*pdfi)[pdfi->size()-1] < 1e-35) 968 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1036 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 969 1037 970 G4double pli = (*pdfi)[0]*factor; 1038 G4double pli = (*pdfi)[0]*factor; 971 G4double pui = (*pdfi)[pdfi->size()-1]*fac 1039 G4double pui = (*pdfi)[pdfi->size()-1]*factor; 972 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal 1040 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal); 973 G4double A_temp = (1.0/(pli*dxLocal))-1.0- 1041 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp; 974 G4double C_temp = 1.0+A_temp+B_temp; 1042 G4double C_temp = 1.0+A_temp+B_temp; 975 if (C_temp < 1e-35) 1043 if (C_temp < 1e-35) 976 { 1044 { 977 (*a)[i]= 0.; 1045 (*a)[i]= 0.; 978 (*b)[i] = 0.; 1046 (*b)[i] = 0.; 979 (*c)[i] = 1; 1047 (*c)[i] = 1; 980 } 1048 } 981 else 1049 else 982 { 1050 { 983 (*a)[i]= A_temp; 1051 (*a)[i]= A_temp; 984 (*b)[i] = B_temp; 1052 (*b)[i] = B_temp; 985 (*c)[i] = C_temp; 1053 (*c)[i] = C_temp; 986 } 1054 } 987 //OK, now get ERR(I), the integral of the 1055 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation 988 //and the true pdf, extended over the inte 1056 //and the true pdf, extended over the interval (X(I),X(I+1)) 989 G4int icase = 1; //loop code 1057 G4int icase = 1; //loop code 990 G4bool reLoop = false; 1058 G4bool reLoop = false; 991 do 1059 do 992 { 1060 { 993 reLoop = false; 1061 reLoop = false; 994 (*err)[i] = 0.; //zero variable 1062 (*err)[i] = 0.; //zero variable 995 for (G4int k=0;k<nip;k++) 1063 for (G4int k=0;k<nip;k++) 996 { 1064 { 997 G4double rr = (*sumi)[k]; 1065 G4double rr = (*sumi)[k]; 998 G4double pap = (*area)[i]*(1.0+((*a)[i]+ 1066 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ 999 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1 1067 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); 1000 if (k == 0 || k == nip-1) 1068 if (k == 0 || k == nip-1) 1001 (*err)[i] += 0.5*std::fabs(pap-(*pdfi 1069 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); 1002 else 1070 else 1003 (*err)[i] += std::fabs(pap-(*pdfi)[k] 1071 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 1004 } 1072 } 1005 (*err)[i] *= dxi; 1073 (*err)[i] *= dxi; 1006 1074 1007 //If err(I) is too large, the pdf is 1075 //If err(I) is too large, the pdf is approximated by a uniform distribution 1008 if ((*err)[i] > 0.1*(*area)[i] && ica 1076 if ((*err)[i] > 0.1*(*area)[i] && icase == 1) 1009 { 1077 { 1010 (*b)[i] = 0; 1078 (*b)[i] = 0; 1011 (*a)[i] = 0; 1079 (*a)[i] = 0; 1012 (*c)[i] = 1.; 1080 (*c)[i] = 1.; 1013 icase = 2; 1081 icase = 2; 1014 reLoop = true; 1082 reLoop = true; 1015 } 1083 } 1016 }while(reLoop); 1084 }while(reLoop); 1017 delete pdfi; 1085 delete pdfi; 1018 delete pdfih; 1086 delete pdfih; 1019 delete sumi; 1087 delete sumi; 1020 } 1088 } 1021 }while(x->size() < np); 1089 }while(x->size() < np); 1022 1090 1023 if (x->size() != np || a->size() != np || 1091 if (x->size() != np || a->size() != np || 1024 err->size() != np || area->size() != np 1092 err->size() != np || area->size() != np) 1025 { 1093 { 1026 G4Exception("G4PenelopeRayleighModel::I 1094 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()", 1027 "em2050",FatalException, 1095 "em2050",FatalException, 1028 "Problem in building the extended Table 1096 "Problem in building the extended Table for Sampling: array dimensions do not match "); 1029 } 1097 } 1030 1098 1031 /****************************************** 1099 /******************************************************************************* 1032 Renormalization 1100 Renormalization 1033 ******************************************* 1101 ********************************************************************************/ 1034 G4double ws = 0; 1102 G4double ws = 0; 1035 for (std::size_t i=0;i<np-1;++i) << 1103 for (size_t i=0;i<np-1;i++) 1036 ws += (*area)[i]; 1104 ws += (*area)[i]; 1037 ws = 1.0/ws; 1105 ws = 1.0/ws; 1038 G4double errMax = 0; 1106 G4double errMax = 0; 1039 for (std::size_t i=0;i<np-1;++i) << 1107 for (size_t i=0;i<np-1;i++) 1040 { 1108 { 1041 (*area)[i] *= ws; 1109 (*area)[i] *= ws; 1042 (*err)[i] *= ws; 1110 (*err)[i] *= ws; 1043 errMax = std::max(errMax,(*err)[i]); 1111 errMax = std::max(errMax,(*err)[i]); 1044 } 1112 } 1045 1113 1046 //Vector with the normalized cumulative dis 1114 //Vector with the normalized cumulative distribution 1047 G4DataVector* PAC = new G4DataVector(); 1115 G4DataVector* PAC = new G4DataVector(); 1048 PAC->push_back(0.); 1116 PAC->push_back(0.); 1049 for (std::size_t i=0;i<np-1;++i) << 1117 for (size_t i=0;i<np-1;i++) 1050 { 1118 { 1051 G4double previous = (*PAC)[i]; 1119 G4double previous = (*PAC)[i]; 1052 PAC->push_back(previous+(*area)[i]); 1120 PAC->push_back(previous+(*area)[i]); 1053 } 1121 } 1054 (*PAC)[PAC->size()-1] = 1.; 1122 (*PAC)[PAC->size()-1] = 1.; 1055 1123 1056 /****************************************** 1124 /******************************************************************************* 1057 Pre-calculated limits for the initial binar 1125 Pre-calculated limits for the initial binary search for subsequent sampling 1058 ******************************************* 1126 ********************************************************************************/ 1059 std::vector<std::size_t> *ITTL = new std::v << 1127 1060 std::vector<std::size_t> *ITTU = new std::v << 1128 //G4DataVector* ITTL = new G4DataVector(); >> 1129 std::vector<size_t> *ITTL = new std::vector<size_t>; >> 1130 std::vector<size_t> *ITTU = new std::vector<size_t>; 1061 1131 1062 //Just create place-holders 1132 //Just create place-holders 1063 for (std::size_t i=0;i<np;++i) << 1133 for (size_t i=0;i<np;i++) 1064 { 1134 { 1065 ITTL->push_back(0); 1135 ITTL->push_back(0); 1066 ITTU->push_back(0); 1136 ITTU->push_back(0); 1067 } 1137 } 1068 1138 1069 G4double bin = 1.0/(np-1); 1139 G4double bin = 1.0/(np-1); 1070 (*ITTL)[0]=0; 1140 (*ITTL)[0]=0; 1071 for (std::size_t i=1;i<(np-1);++i) << 1141 for (size_t i=1;i<(np-1);i++) 1072 { 1142 { 1073 G4double ptst = i*bin; 1143 G4double ptst = i*bin; 1074 G4bool found = false; 1144 G4bool found = false; 1075 for (std::size_t j=(*ITTL)[i-1];j<np && << 1145 for (size_t j=(*ITTL)[i-1];j<np && !found;j++) 1076 { 1146 { 1077 if ((*PAC)[j] > ptst) 1147 if ((*PAC)[j] > ptst) 1078 { 1148 { 1079 (*ITTL)[i] = j-1; 1149 (*ITTL)[i] = j-1; 1080 (*ITTU)[i-1] = j; 1150 (*ITTU)[i-1] = j; 1081 found = true; 1151 found = true; 1082 } 1152 } 1083 } 1153 } 1084 } 1154 } 1085 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; 1155 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; 1086 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; 1156 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; 1087 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; 1157 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; 1088 1158 1089 if (ITTU->size() != np || ITTU->size() != n 1159 if (ITTU->size() != np || ITTU->size() != np) 1090 { 1160 { 1091 G4Exception("G4PenelopeRayleighModel::I 1161 G4Exception("G4PenelopeRayleighModel::InitializeSamplingAlgorithm()", 1092 "em2051",FatalException, 1162 "em2051",FatalException, 1093 "Problem in building the Limit Tables f 1163 "Problem in building the Limit Tables for Sampling: array dimensions do not match"); 1094 } 1164 } 1095 1165 >> 1166 1096 /****************************************** 1167 /******************************************************************************** 1097 Copy tables 1168 Copy tables 1098 ******************************************* 1169 ********************************************************************************/ 1099 G4PenelopeSamplingData* theTable = new G4Pe 1170 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np); 1100 for (std::size_t i=0;i<np;++i) << 1171 for (size_t i=0;i<np;i++) 1101 { 1172 { 1102 theTable->AddPoint((*x)[i],(*PAC)[i],(* 1173 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]); 1103 } 1174 } 1104 1175 1105 if (fVerboseLevel > 2) << 1176 if (verboseLevel > 2) 1106 { 1177 { 1107 G4cout << "**************************** 1178 G4cout << "*************************************************************************" << 1108 G4endl; 1179 G4endl; 1109 G4cout << "Sampling table for Penelope 1180 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl; 1110 theTable->DumpTable(); 1181 theTable->DumpTable(); 1111 } 1182 } 1112 fSamplingTable->insert(std::make_pair(mat,t << 1183 samplingTable->insert(std::make_pair(mat,theTable)); >> 1184 1113 1185 1114 //Clean up temporary vectors 1186 //Clean up temporary vectors 1115 delete x; 1187 delete x; 1116 delete a; 1188 delete a; 1117 delete b; 1189 delete b; 1118 delete c; 1190 delete c; 1119 delete err; 1191 delete err; 1120 delete area; 1192 delete area; 1121 delete PAC; 1193 delete PAC; 1122 delete ITTL; 1194 delete ITTL; 1123 delete ITTU; 1195 delete ITTU; 1124 1196 1125 //DONE! 1197 //DONE! 1126 return; 1198 return; >> 1199 1127 } 1200 } 1128 1201 1129 //....oooOO0OOooo........oooOO0OOooo........o 1202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1130 1203 1131 void G4PenelopeRayleighModel::GetPMaxTable(co 1204 void G4PenelopeRayleighModel::GetPMaxTable(const G4Material* mat) 1132 { 1205 { 1133 if (!fPMaxTable) << 1206 >> 1207 if (!pMaxTable) 1134 { 1208 { 1135 G4cout << "G4PenelopeRayleighModel::Bui 1209 G4cout << "G4PenelopeRayleighModel::BuildPMaxTable" << G4endl; 1136 G4cout << "Going to instanziate the fPM << 1210 G4cout << "Going to instanziate the pMaxTable !" << G4endl; 1137 G4cout << "That should _not_ be here! " 1211 G4cout << "That should _not_ be here! " << G4endl; 1138 fPMaxTable = new std::map<const G4Mater << 1212 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; 1139 } 1213 } 1140 //check if the table is already there 1214 //check if the table is already there 1141 if (fPMaxTable->count(mat)) << 1215 if (pMaxTable->count(mat)) 1142 return; 1216 return; 1143 1217 1144 //otherwise build it 1218 //otherwise build it 1145 if (!fSamplingTable) << 1219 if (!samplingTable) 1146 { 1220 { 1147 G4Exception("G4PenelopeRayleighModel::G 1221 G4Exception("G4PenelopeRayleighModel::GetPMaxTable()", 1148 "em2052",FatalException, 1222 "em2052",FatalException, 1149 "SamplingTable is not properly instanti 1223 "SamplingTable is not properly instantiated"); 1150 return; 1224 return; 1151 } 1225 } 1152 1226 1153 //This should not be: the sampling table is 1227 //This should not be: the sampling table is built before the p-table 1154 if (!fSamplingTable->count(mat)) << 1228 if (!samplingTable->count(mat)) 1155 { 1229 { 1156 G4ExceptionDescription ed; 1230 G4ExceptionDescription ed; 1157 ed << "Sampling table for material " < 1231 ed << "Sampling table for material " << mat->GetName() << " not found"; 1158 G4Exception("G4PenelopeRayleighModel:: 1232 G4Exception("G4PenelopeRayleighModel::GetPMaxTable()", 1159 "em2052",FatalException, 1233 "em2052",FatalException, 1160 ed); 1234 ed); 1161 return; 1235 return; 1162 } 1236 } 1163 1237 1164 G4PenelopeSamplingData *theTable = fSamplin << 1238 G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second; 1165 std::size_t tablePoints = theTable->GetNumb << 1239 size_t tablePoints = theTable->GetNumberOfStoredPoints(); 1166 1240 1167 std::size_t nOfEnergyPoints = fLogEnergyGri << 1241 size_t nOfEnergyPoints = logEnergyGridPMax.size(); 1168 G4PhysicsFreeVector* theVec = new G4Physics 1242 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints); 1169 1243 1170 const std::size_t nip = 51; //hard-coded in << 1244 const size_t nip = 51; //hard-coded in Penelope 1171 1245 1172 for (std::size_t ie=0;ie<fLogEnergyGridPMax << 1246 for (size_t ie=0;ie<logEnergyGridPMax.size();ie++) 1173 { 1247 { 1174 G4double energy = G4Exp(fLogEnergyGridP << 1248 G4double energy = G4Exp(logEnergyGridPMax[ie]); 1175 G4double Qm = 2.0*energy/electron_mass_ 1249 G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now 1176 G4double Qm2 = Qm*Qm; 1250 G4double Qm2 = Qm*Qm; 1177 G4double firstQ2 = theTable->GetX(0); 1251 G4double firstQ2 = theTable->GetX(0); 1178 G4double lastQ2 = theTable->GetX(tableP 1252 G4double lastQ2 = theTable->GetX(tablePoints-1); 1179 G4double thePMax = 0; 1253 G4double thePMax = 0; 1180 1254 1181 if (Qm2 > firstQ2) 1255 if (Qm2 > firstQ2) 1182 { 1256 { 1183 if (Qm2 < lastQ2) 1257 if (Qm2 < lastQ2) 1184 { 1258 { 1185 //bisection to look for the index of 1259 //bisection to look for the index of Qm 1186 std::size_t lowerBound = 0; << 1260 size_t lowerBound = 0; 1187 std::size_t upperBound = tablePoints- << 1261 size_t upperBound = tablePoints-1; 1188 while (lowerBound <= upperBound) 1262 while (lowerBound <= upperBound) 1189 { 1263 { 1190 std::size_t midBin = (lowerBound + uppe << 1264 size_t midBin = (lowerBound + upperBound)/2; 1191 if( Qm2 < theTable->GetX(midBin)) 1265 if( Qm2 < theTable->GetX(midBin)) 1192 { upperBound = midBin-1; } 1266 { upperBound = midBin-1; } 1193 else 1267 else 1194 { lowerBound = midBin+1; } 1268 { lowerBound = midBin+1; } 1195 } 1269 } 1196 //upperBound is the output (but also 1270 //upperBound is the output (but also lowerBounf --> should be the same!) 1197 G4double Q1 = theTable->GetX(upperBou 1271 G4double Q1 = theTable->GetX(upperBound); 1198 G4double Q2 = Qm2; 1272 G4double Q2 = Qm2; 1199 G4double DQ = (Q2-Q1)/((G4double)(nip 1273 G4double DQ = (Q2-Q1)/((G4double)(nip-1)); 1200 G4double theA = theTable->GetA(upperB 1274 G4double theA = theTable->GetA(upperBound); 1201 G4double theB = theTable->GetB(upperB 1275 G4double theB = theTable->GetB(upperBound); 1202 G4double thePAC = theTable->GetPAC(up 1276 G4double thePAC = theTable->GetPAC(upperBound); 1203 G4DataVector* fun = new G4DataVector( 1277 G4DataVector* fun = new G4DataVector(); 1204 for (std::size_t k=0;k<nip;++k) << 1278 for (size_t k=0;k<nip;k++) 1205 { 1279 { 1206 G4double qi = Q1 + k*DQ; 1280 G4double qi = Q1 + k*DQ; 1207 G4double tau = (qi-Q1)/ 1281 G4double tau = (qi-Q1)/ 1208 (theTable->GetX(upperBound+1)-Q1); 1282 (theTable->GetX(upperBound+1)-Q1); 1209 G4double con1 = 2.0*theB*tau; 1283 G4double con1 = 2.0*theB*tau; 1210 G4double ci = 1.0+theA+theB; 1284 G4double ci = 1.0+theA+theB; 1211 G4double con2 = ci-theA*tau; 1285 G4double con2 = ci-theA*tau; 1212 G4double etap = 0; 1286 G4double etap = 0; 1213 if (std::fabs(con1) > 1.0e-16*std::fabs 1287 if (std::fabs(con1) > 1.0e-16*std::fabs(con2)) 1214 etap = con2*(1.0-std::sqrt(1.0-2.0*ta 1288 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1; 1215 else 1289 else 1216 etap = tau/con2; 1290 etap = tau/con2; 1217 G4double theFun = (theTable->GetPAC(upp 1291 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)* 1218 (1.0+(theA+theB*etap)*etap)*(1.0+(the 1292 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/ 1219 ((1.0-theB*etap*etap)*ci*(theTable->G 1293 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1)); 1220 fun->push_back(theFun); 1294 fun->push_back(theFun); 1221 } 1295 } 1222 //Now intergrate numerically the fun 1296 //Now intergrate numerically the fun Cavalieri-Simpson's method 1223 G4DataVector* sum = new G4DataVector; 1297 G4DataVector* sum = new G4DataVector; 1224 G4double CONS = DQ*(1./12.); 1298 G4double CONS = DQ*(1./12.); 1225 G4double HCONS = 0.5*CONS; 1299 G4double HCONS = 0.5*CONS; 1226 sum->push_back(0.); 1300 sum->push_back(0.); 1227 G4double secondPoint = (*sum)[0] + 1301 G4double secondPoint = (*sum)[0] + 1228 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C 1302 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS; 1229 sum->push_back(secondPoint); 1303 sum->push_back(secondPoint); 1230 for (std::size_t hh=2;hh<nip-1;++hh) << 1304 for (size_t hh=2;hh<nip-1;hh++) 1231 { 1305 { 1232 G4double previous = (*sum)[hh-1]; 1306 G4double previous = (*sum)[hh-1]; 1233 G4double next = previous+(13.0*((*fun)[ 1307 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])- 1234 (*fun)[hh+1]-(*fun)[hh-2])*HCON 1308 (*fun)[hh+1]-(*fun)[hh-2])*HCONS; 1235 sum->push_back(next); 1309 sum->push_back(next); 1236 } 1310 } 1237 G4double last = (*sum)[nip-2]+(5.0*(* 1311 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]- 1238 (*fun)[nip-3])*CONS; 1312 (*fun)[nip-3])*CONS; 1239 sum->push_back(last); 1313 sum->push_back(last); 1240 thePMax = thePAC + (*sum)[sum->size() 1314 thePMax = thePAC + (*sum)[sum->size()-1]; //last point 1241 delete fun; 1315 delete fun; 1242 delete sum; 1316 delete sum; 1243 } 1317 } 1244 else 1318 else 1245 { 1319 { 1246 thePMax = 1.0; 1320 thePMax = 1.0; 1247 } 1321 } 1248 } 1322 } 1249 else 1323 else 1250 { 1324 { 1251 thePMax = theTable->GetPAC(0); 1325 thePMax = theTable->GetPAC(0); 1252 } 1326 } 1253 1327 1254 //Write number in the table 1328 //Write number in the table 1255 theVec->PutValue(ie,energy,thePMax); 1329 theVec->PutValue(ie,energy,thePMax); 1256 } 1330 } 1257 1331 1258 fPMaxTable->insert(std::make_pair(mat,theVe << 1332 pMaxTable->insert(std::make_pair(mat,theVec)); 1259 return; 1333 return; >> 1334 1260 } 1335 } 1261 1336 1262 //....oooOO0OOooo........oooOO0OOooo........o 1337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1263 1338 1264 void G4PenelopeRayleighModel::DumpFormFactorT 1339 void G4PenelopeRayleighModel::DumpFormFactorTable(const G4Material* mat) 1265 { 1340 { 1266 G4cout << "******************************** 1341 G4cout << "*****************************************************************" << G4endl; 1267 G4cout << "G4PenelopeRayleighModel: Form Fa 1342 G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl; 1268 //try to use the same format as Penelope-Fo 1343 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F 1269 G4cout << "Q/(m_e*c) F(Q) 1344 G4cout << "Q/(m_e*c) F(Q) " << G4endl; 1270 G4cout << "******************************** 1345 G4cout << "*****************************************************************" << G4endl; 1271 if (!fLogFormFactorTable->count(mat)) << 1346 if (!logFormFactorTable->count(mat)) 1272 BuildFormFactorTable(mat); 1347 BuildFormFactorTable(mat); 1273 1348 1274 G4PhysicsFreeVector* theVec = fLogFormFacto << 1349 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; 1275 for (std::size_t i=0;i<theVec->GetVectorLen << 1350 for (size_t i=0;i<theVec->GetVectorLength();i++) 1276 { 1351 { 1277 G4double logQ2 = theVec->GetLowEdgeEner 1352 G4double logQ2 = theVec->GetLowEdgeEnergy(i); 1278 G4double Q = G4Exp(0.5*logQ2); 1353 G4double Q = G4Exp(0.5*logQ2); 1279 G4double logF2 = (*theVec)[i]; 1354 G4double logF2 = (*theVec)[i]; 1280 G4double F = G4Exp(0.5*logF2); 1355 G4double F = G4Exp(0.5*logF2); 1281 G4cout << Q << " " << F << 1356 G4cout << Q << " " << F << G4endl; 1282 } 1357 } 1283 //DONE 1358 //DONE 1284 return; 1359 return; 1285 } 1360 } 1286 1361 1287 //....oooOO0OOooo........oooOO0OOooo........o 1362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 1288 1363 1289 void G4PenelopeRayleighModel::SetParticle(con 1364 void G4PenelopeRayleighModel::SetParticle(const G4ParticleDefinition* p) 1290 { 1365 { 1291 if(!fParticle) { 1366 if(!fParticle) { 1292 fParticle = p; 1367 fParticle = p; 1293 } 1368 } 1294 } 1369 } 1295 1370