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: G4PenelopeGammaConversionModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $ >> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $ 26 // 28 // 27 // Author: Luciano Pandola 29 // Author: Luciano Pandola 28 // 30 // 29 // History: 31 // History: 30 // -------- 32 // -------- 31 // 13 Jan 2010 L Pandola First implementa << 33 // 06 Oct 2008 L Pandola Migration from process to model 32 // 24 May 2011 L Pandola Renamed (make v2 << 33 // 18 Sep 2013 L Pandola Migration to MT << 34 // data and create << 35 // 34 // 36 35 37 #include "G4PenelopeGammaConversionModel.hh" 36 #include "G4PenelopeGammaConversionModel.hh" 38 #include "G4PhysicalConstants.hh" << 39 #include "G4SystemOfUnits.hh" << 40 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh" 41 #include "G4MaterialCutsCouple.hh" 38 #include "G4MaterialCutsCouple.hh" 42 #include "G4ProductionCutsTable.hh" 39 #include "G4ProductionCutsTable.hh" 43 #include "G4DynamicParticle.hh" 40 #include "G4DynamicParticle.hh" 44 #include "G4Element.hh" 41 #include "G4Element.hh" 45 #include "G4Gamma.hh" 42 #include "G4Gamma.hh" 46 #include "G4Electron.hh" 43 #include "G4Electron.hh" 47 #include "G4Positron.hh" 44 #include "G4Positron.hh" 48 #include "G4PhysicsFreeVector.hh" << 45 #include "G4CrossSectionHandler.hh" 49 #include "G4MaterialCutsCouple.hh" << 50 #include "G4AutoLock.hh" << 51 #include "G4Exp.hh" << 52 46 53 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 54 const G4int G4PenelopeGammaConversionModel::fM << 55 G4PhysicsFreeVector* G4PenelopeGammaConversion << 56 G4double G4PenelopeGammaConversionModel::fAtom << 57 1.2281e+02,7.3167e+01,6.9 << 58 6.4696e+01,6.1228e+01,5.7 << 59 5.0787e+01,4.7851e+01,4.6 << 60 4.4503e+01,4.3815e+01,4.3 << 61 4.1586e+01,4.0953e+01,4.0 << 62 3.9756e+01,3.9144e+01,3.8 << 63 3.7174e+01,3.6663e+01,3.5 << 64 3.4688e+01,3.4197e+01,3.3 << 65 3.3068e+01,3.2740e+01,3.2 << 66 3.1884e+01,3.1622e+01,3.1 << 67 3.0950e+01,3.0758e+01,3.0 << 68 3.0097e+01,2.9832e+01,2.9 << 69 2.9247e+01,2.9085e+01,2.8 << 70 2.8580e+01,2.8442e+01,2.8 << 71 2.7973e+01,2.7819e+01,2.7 << 72 2.7285e+01,2.7093e+01,2.6 << 73 2.6516e+01,2.6304e+01,2.6 << 74 2.5730e+01,2.5577e+01,2.5 << 75 2.5100e+01,2.4941e+01,2.4 << 76 2.4506e+01,2.4391e+01,2.4 << 77 2.4039e+01,2.3922e+01,2.3 << 78 2.3621e+01,2.3523e+01,2.3 << 79 2.3238e+01,2.3139e+01,2.3 << 80 2.2833e+01,2.2694e+01,2.2 << 81 2.2446e+01,2.2358e+01,2.2 << 82 48 83 //....oooOO0OOooo........oooOO0OOooo........oo << 84 49 85 G4PenelopeGammaConversionModel::G4PenelopeGamm << 50 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*, 86 const G4String& nam) << 51 const G4String& nam) 87 :G4VEmModel(nam),fParticleChange(nullptr),fP << 52 :G4VEmModel(nam),fTheScreeningRadii(0),crossSectionHandler(0),isInitialised(false) 88 fEffectiveCharge(nullptr),fMaterialInvScree << 89 fScreeningFunction(nullptr),fIsInitialised( << 90 { 53 { 91 fIntrinsicLowEnergyLimit = 2.0*electron_mass 54 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2; 92 fIntrinsicHighEnergyLimit = 100.0*GeV; 55 fIntrinsicHighEnergyLimit = 100.0*GeV; 93 fSmallEnergy = 1.1*MeV; 56 fSmallEnergy = 1.1*MeV; 94 57 95 if (part) << 58 SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 96 SetParticle(part); << 97 << 98 // SetLowEnergyLimit(fIntrinsicLowEnergyLim << 99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 59 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 100 // 60 // 101 fVerboseLevel= 0; << 61 verboseLevel= 0; 102 // Verbosity scale: 62 // Verbosity scale: 103 // 0 = nothing << 63 // 0 = nothing 104 // 1 = warning for energy non-conservation << 64 // 1 = warning for energy non-conservation 105 // 2 = details of energy budget 65 // 2 = details of energy budget 106 // 3 = calculation of cross sections, file o 66 // 3 = calculation of cross sections, file openings, sampling of atoms 107 // 4 = entering in methods 67 // 4 = entering in methods 108 } 68 } 109 69 110 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 111 71 112 G4PenelopeGammaConversionModel::~G4PenelopeGam 72 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel() 113 { << 73 { 114 //Delete shared tables, they exist only in t << 74 if (crossSectionHandler) delete crossSectionHandler; 115 if (IsMaster() || fLocalTable) << 75 if (fTheScreeningRadii) delete fTheScreeningRadii; 116 { << 117 for(G4int i=0; i<=fMaxZ; ++i) << 118 { << 119 if(fLogAtomicCrossSection[i]) { << 120 delete fLogAtomicCrossSection[i]; << 121 fLogAtomicCrossSection[i] = nullptr; << 122 } << 123 } << 124 if (fEffectiveCharge) << 125 delete fEffectiveCharge; << 126 if (fMaterialInvScreeningRadius) << 127 delete fMaterialInvScreeningRadius; << 128 if (fScreeningFunction) << 129 delete fScreeningFunction; << 130 } << 131 } 76 } 132 77 133 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 134 79 135 void G4PenelopeGammaConversionModel::Initialis << 80 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*, 136 const G4DataVector&) << 81 const G4DataVector& ) 137 { 82 { 138 if (fVerboseLevel > 3) << 83 if (verboseLevel > 3) 139 G4cout << "Calling G4PenelopeGammaConvers 84 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl; 140 85 141 SetParticle(part); << 86 //Delete the old cross section handler, if necessary 142 << 87 if (crossSectionHandler) 143 //Only the master model creates/fills/destro << 144 if (IsMaster() && part == fParticle) << 145 { 88 { 146 //delete old material data... << 89 crossSectionHandler->Clear(); 147 if (fEffectiveCharge) << 90 delete crossSectionHandler; 148 { << 91 } 149 delete fEffectiveCharge; << 92 150 fEffectiveCharge = nullptr; << 93 //Check energy limits 151 } << 94 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit) 152 if (fMaterialInvScreeningRadius) << 95 { 153 { << 96 G4cout << "G4PenelopeGammaConversionModel: low energy limit increased from " << 154 delete fMaterialInvScreeningRadius; << 97 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl; 155 fMaterialInvScreeningRadius = nullptr; << 98 SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 156 } << 99 } 157 if (fScreeningFunction) << 100 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit) 158 { << 101 { 159 delete fScreeningFunction; << 102 G4cout << "G4PenelopeGammaConversionModel: high energy limit decreased from " << 160 fScreeningFunction = nullptr; << 103 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl; 161 } << 104 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 162 //and create new ones << 163 fEffectiveCharge = new std::map<const G4 << 164 fMaterialInvScreeningRadius = new std::m << 165 fScreeningFunction = new std::map<const << 166 << 167 G4ProductionCutsTable* theCoupleTable = << 168 G4ProductionCutsTable::GetProductionCutsTabl << 169 << 170 for (G4int i=0;i<(G4int)theCoupleTable-> << 171 { << 172 const G4Material* material = << 173 theCoupleTable->GetMaterialCutsCouple(i) << 174 const G4ElementVector* theElementVector = << 175 << 176 for (std::size_t j=0;j<material->GetNumber << 177 { << 178 G4int iZ = theElementVector->at(j)->Ge << 179 //read data files only in the master << 180 if (iZ <= fMaxZ && !fLogAtomicCrossSe << 181 ReadDataFile(iZ); << 182 } << 183 << 184 //check if material data are available << 185 if (!fEffectiveCharge->count(material)) << 186 InitializeScreeningFunctions(material); << 187 } << 188 if (fVerboseLevel > 0) { << 189 G4cout << "Penelope Gamma Conversion model v << 190 << "Energy range: " << 191 << LowEnergyLimit() / MeV << " MeV - << 192 << HighEnergyLimit() / GeV << " GeV" << 193 << G4endl; << 194 } << 195 } 105 } 196 if(fIsInitialised) return; << 197 fParticleChange = GetParticleChangeForGamma( << 198 fIsInitialised = true; << 199 } << 200 106 201 //....oooOO0OOooo........oooOO0OOooo........oo << 107 //Re-initialize cross section handler >> 108 crossSectionHandler = new G4CrossSectionHandler(); >> 109 crossSectionHandler->Initialise(0,LowEnergyLimit(),HighEnergyLimit(),400); >> 110 crossSectionHandler->Clear(); >> 111 G4String crossSectionFile = "penelope/pp-cs-pen-"; >> 112 crossSectionHandler->LoadData(crossSectionFile); >> 113 //This is used to retrieve cross section values later on >> 114 crossSectionHandler->BuildMeanFreePathForMaterials(); 202 115 203 void G4PenelopeGammaConversionModel::Initialis << 116 if (verboseLevel > 2) 204 G4VEmModel *masterModel) << 117 G4cout << "Loaded cross section files for PenelopeGammaConversion" << G4endl; 205 { << 206 if (fVerboseLevel > 3) << 207 G4cout << "Calling G4PenelopeGammaConvers << 208 // << 209 //Check that particle matches: one might hav << 210 //for e+ and e-). << 211 // << 212 if (part == fParticle) << 213 { << 214 //Get the const table pointers from the << 215 const G4PenelopeGammaConversionModel* th << 216 static_cast<G4PenelopeGammaConversionModel*> << 217 118 218 //Copy pointers to the data tables << 119 G4cout << "Penelope Gamma Conversion model is initialized " << G4endl 219 fEffectiveCharge = theModel->fEffectiveC << 120 << "Energy range: " 220 fMaterialInvScreeningRadius = theModel-> << 121 << LowEnergyLimit() / MeV << " MeV - " 221 fScreeningFunction = theModel->fScreenin << 122 << HighEnergyLimit() / GeV << " GeV" 222 for(G4int i=0; i<=fMaxZ; ++i) << 123 << G4endl; 223 fLogAtomicCrossSection[i] = theModel->fLogAt << 224 124 225 //Same verbosity for all workers, as the << 125 if(isInitialised) return; 226 fVerboseLevel = theModel->fVerboseLevel; << 227 } << 228 126 229 return; << 127 if(pParticleChange) >> 128 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); >> 129 else >> 130 fParticleChange = new G4ParticleChangeForGamma(); >> 131 isInitialised = true; 230 } 132 } 231 133 232 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 233 namespace { G4Mutex PenelopeGammaConversionMo << 234 135 235 G4double G4PenelopeGammaConversionModel::Compu 136 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom( 236 const G4ParticleDefinition << 137 const G4ParticleDefinition*, 237 G4double energy, << 138 G4double energy, 238 G4double Z, G4double, << 139 G4double Z, G4double, 239 G4double, G4double) << 140 G4double, G4double) 240 { 141 { 241 // 142 // 242 // Penelope model v2008. << 143 // Penelope model. 243 // Cross section (including triplet producti << 144 // Cross section (including triplet production) read from database and managed 244 // through the G4CrossSectionHandler utility 145 // through the G4CrossSectionHandler utility. Cross section data are from 245 // M.J. Berger and J.H. Hubbel (XCOM), Repor 146 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598 246 // 147 // >> 148 >> 149 if (verboseLevel > 3) >> 150 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl; 247 151 248 if (energy < fIntrinsicLowEnergyLimit) << 152 G4int iZ = (G4int) Z; 249 return 0; << 153 if (!crossSectionHandler) 250 << 154 { 251 G4int iZ = G4int(Z); << 155 G4cout << "G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom" << G4endl; 252 << 156 G4cout << "The cross section handler is not correctly initialized" << G4endl; 253 if (!fLogAtomicCrossSection[iZ]) << 157 G4Exception(); 254 { << 158 } 255 //If we are here, it means that Initial << 159 G4double cs = crossSectionHandler->FindValue(iZ,energy); 256 //not filled up. This can happen in a U << 257 if (fVerboseLevel > 0) << 258 { << 259 //Issue a G4Exception (warning) only in v << 260 G4ExceptionDescription ed; << 261 ed << "Unable to retrieve the cross secti << 262 ed << "This can happen only in Unit Tests << 263 G4Exception("G4PenelopeGammaConversionMod << 264 "em2018",JustWarning,ed); << 265 } << 266 //protect file reading via autolock << 267 G4AutoLock lock(&PenelopeGammaConversio << 268 ReadDataFile(iZ); << 269 lock.unlock(); << 270 fLocalTable = true; << 271 } << 272 G4double cs = 0; << 273 G4double logene = G4Log(energy); << 274 G4PhysicsFreeVector* theVec = fLogAtomicCros << 275 G4double logXS = theVec->Value(logene); << 276 cs = G4Exp(logXS); << 277 160 278 if (fVerboseLevel > 2) << 161 if (verboseLevel > 2) 279 G4cout << "Gamma conversion cross section << 162 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z << 280 " = " << cs/barn << " barn" << G4endl; 163 " = " << cs/barn << " barn" << G4endl; 281 return cs; 164 return cs; 282 } 165 } 283 166 284 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 285 168 286 void << 169 void G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 287 G4PenelopeGammaConversionModel::SampleSecondar << 170 const G4MaterialCutsCouple* couple, 288 const G4MaterialCutsCouple* coup << 171 const G4DynamicParticle* aDynamicGamma, 289 const G4DynamicParticle* aDynami << 172 G4double, 290 G4double, << 173 G4double) 291 G4double) << 292 { 174 { 293 // 175 // 294 // Penelope model v2008. << 176 // Penelope model. 295 // Final state is sampled according to the B 177 // Final state is sampled according to the Bethe-Heitler model with Coulomb 296 // corrections, according to the semi-empiri 178 // corrections, according to the semi-empirical model of 297 // J. Baro' et al., Radiat. Phys. Chem. 44 179 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531. 298 // 180 // 299 // The model uses the high energy Coulomb co << 181 // The model uses the high energy Coulomb correction from 300 // H. Davies et al., Phys. Rev. 93 (1954) 7 182 // H. Davies et al., Phys. Rev. 93 (1954) 788 301 // and atomic screening radii tabulated from << 183 // and atomic screening radii tabulated from 302 // J.H. Hubbel et al., J. Phys. Chem. Ref. 184 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023 303 // for Z= 1 to 92. << 185 // for Z= 1 to 92. This managed in this model by the method >> 186 // GetScreeningRadius(). 304 // 187 // 305 if (fVerboseLevel > 3) << 188 if (verboseLevel > 3) 306 G4cout << "Calling SamplingSecondaries() o 189 G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl; 307 190 308 G4double photonEnergy = aDynamicGamma->GetKi 191 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 309 192 310 // Always kill primary << 193 if (photonEnergy <= LowEnergyLimit()) 311 fParticleChange->ProposeTrackStatus(fStopAnd << 194 { 312 fParticleChange->SetProposedKineticEnergy(0. << 195 fParticleChange->ProposeTrackStatus(fStopAndKill); 313 << 196 fParticleChange->SetProposedKineticEnergy(0.); 314 if (photonEnergy <= fIntrinsicLowEnergyLimit << 315 { << 316 fParticleChange->ProposeLocalEnergyDepos 197 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy); 317 return ; 198 return ; 318 } << 199 } 319 200 320 G4ParticleMomentum photonDirection = aDynami 201 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); 321 const G4Material* mat = couple->GetMaterial( << 322 202 323 //Either Initialize() was not called, or we << 203 G4double eps ; 324 //not invoked << 204 G4double eki = electron_mass_c2 / photonEnergy ; 325 if (!fEffectiveCharge) << 326 { << 327 //create a **thread-local** version of t << 328 //Unit Tests << 329 fLocalTable = true; << 330 fEffectiveCharge = new std::map<const G4 << 331 fMaterialInvScreeningRadius = new std::m << 332 fScreeningFunction = new std::map<const << 333 } << 334 205 335 if (!fEffectiveCharge->count(mat)) << 206 // Do it fast if photon energy < 1.1 MeV >> 207 if (photonEnergy < fSmallEnergy ) 336 { 208 { 337 //If we are here, it means that Initiali << 209 eps = eki + (1-2*eki) * G4UniformRand(); 338 //not filled up. This can happen in a Un << 339 if (fVerboseLevel > 0) << 340 { << 341 //Issue a G4Exception (warning) only in ve << 342 G4ExceptionDescription ed; << 343 ed << "Unable to allocate the EffectiveCha << 344 mat->GetName() << G4endl; << 345 ed << "This can happen only in Unit Tests" << 346 G4Exception("G4PenelopeGammaConversionMode << 347 "em2019",JustWarning,ed); << 348 } << 349 //protect file reading via autolock << 350 G4AutoLock lock(&PenelopeGammaConversion << 351 InitializeScreeningFunctions(mat); << 352 lock.unlock(); << 353 } 210 } 354 << 355 // eps is the fraction of the photon energy << 356 G4double eps = 0; << 357 G4double eki = electron_mass_c2/photonEnergy << 358 << 359 //Do it fast for photon energy < 1.1 MeV (cl << 360 if (photonEnergy < fSmallEnergy) << 361 eps = eki + (1.0-2.0*eki)*G4UniformRand(); << 362 else 211 else 363 { 212 { 364 //Complete calculation << 213 // Select randomly one element in the current material 365 G4double effC = fEffectiveCharge->find(m << 214 if (verboseLevel > 2) 366 G4double alz = effC*fine_structure_const << 215 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl; 367 G4double T = std::sqrt(2.0*eki); << 216 //use crossSectionHandler instead of G4EmElementSelector because in this case 368 G4double F00=(-1.774-1.210e1*alz+1.118e1 << 217 //the dimension of the table is equal to the dimension of the database 369 +(8.523+7.326e1*alz-4.441e1*alz*alz)* << 218 //(less interpolation errors) 370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz << 219 G4int Z_int = crossSectionHandler->SelectRandomAtom(couple,photonEnergy); 371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T << 220 if (verboseLevel > 2) 372 << 221 G4cout << "Selected Z = " << Z_int << G4endl; 373 G4double F0b = fScreeningFunction->find( << 222 374 G4double g0 = F0b + F00; << 223 //Low energy and Coulomb corrections 375 G4double invRad = fMaterialInvScreeningR << 224 G4double Z=(G4double) Z_int; 376 G4double bmin = 4.0*eki/invRad; << 225 G4double ZAlpha = Z*fine_structure_const; 377 std::pair<G4double,G4double> scree = Ge << 226 G4double ScreenRadius = GetScreeningRadius(Z); 378 G4double g1 = scree.first; << 227 G4double funct1=0,g0=0; 379 G4double g2 = scree.second; << 228 G4double g1min=0,g2min=0; 380 G4double g1min = g1+g0; << 229 funct1 = 4.0*std::log(ScreenRadius); 381 G4double g2min = g2+g0; << 230 g0 = funct1-4*CoulombCorrection(ZAlpha)+LowEnergyCorrection(ZAlpha,eki); 382 G4double xr = 0.5-eki; << 231 G4double bmin = 2*eki*ScreenRadius; 383 G4double a1 = 2.*g1min*xr*xr/3.; << 232 std::vector<G4double> ScreenFunctionValues = ScreenFunction(bmin); 384 G4double p1 = a1/(a1+g2min); << 233 if (ScreenFunctionValues.size() != 2) >> 234 { >> 235 G4cout << "G4PenelopeGammaConversionModel::SampleSecondaries" << G4endl; >> 236 G4cout << "ScreenFunction did not return 2 values! Something wrong! " << G4endl; >> 237 G4Exception(); >> 238 } >> 239 g1min=g0+ScreenFunctionValues[0]; >> 240 g2min=g0+ScreenFunctionValues[1]; >> 241 G4double xr,a1,p1; >> 242 xr=0.5-eki; >> 243 a1=(2.0/3.0)*g1min*xr*xr; >> 244 p1=a1/(a1+g2min); 385 245 386 G4bool loopAgain = false; << 387 //Random sampling of eps 246 //Random sampling of eps >> 247 G4double rand1,rand2,rand3,b; >> 248 G4double g1; >> 249 388 do{ 250 do{ 389 loopAgain = false; << 251 rand1 = G4UniformRand(); 390 if (G4UniformRand() <= p1) << 252 if (rand1 < p1) { 391 { << 253 rand2 = 2.0*G4UniformRand()-1.0; 392 G4double ru2m1 = 2.0*G4UniformRand()-1. << 254 if (rand2 < 0) { 393 if (ru2m1 < 0) << 255 eps = 0.5 - xr*std::pow(std::abs(rand2),(1./3.)); 394 eps = 0.5-xr*std::pow(std::abs(ru2m1), << 395 else << 396 eps = 0.5+xr*std::pow(ru2m1,1./3.); << 397 G4double B = eki/(invRad*eps*(1.0-eps)); << 398 scree = GetScreeningFunctions(B); << 399 g1 = scree.first; << 400 g1 = std::max(g1+g0,0.); << 401 if (G4UniformRand()*g1min > g1) << 402 loopAgain = true; << 403 } 256 } 404 else << 257 else >> 258 { >> 259 eps = 0.5 + xr*std::pow(rand2,(1./3.)); >> 260 } >> 261 b = (eki*ScreenRadius)/(2*eps*(1.0-eps)); >> 262 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b); >> 263 g1 = g0+ScreenFunctionSampling[0]; >> 264 if (g1 < 0) g1=0; >> 265 rand3 = G4UniformRand()*g1min; >> 266 } >> 267 else 405 { 268 { 406 eps = eki+2.0*xr*G4UniformRand(); 269 eps = eki+2.0*xr*G4UniformRand(); 407 G4double B = eki/(invRad*eps*(1.0-eps)); << 270 b = (eki*ScreenRadius)/(2*eps*(1.0-eps)); 408 scree = GetScreeningFunctions(B); << 271 std::vector<G4double> ScreenFunctionSampling = ScreenFunction(b); 409 g2 = scree.second; << 272 g1 = g0+ScreenFunctionSampling[1]; 410 g2 = std::max(g2+g0,0.); << 273 if (g1 < 0) g1=0; 411 if (G4UniformRand()*g2min > g2) << 274 rand3 = G4UniformRand()*g2min; 412 loopAgain = true; << 275 } 413 } << 276 } while (rand3>g1); 414 }while(loopAgain); << 277 } //End of eps sampling 415 } << 416 if (fVerboseLevel > 4) << 417 G4cout << "Sampled eps = " << eps << G4end << 418 278 419 G4double electronTotEnergy = eps*photonEnerg 279 G4double electronTotEnergy = eps*photonEnergy; 420 G4double positronTotEnergy = (1.0-eps)*photo 280 G4double positronTotEnergy = (1.0-eps)*photonEnergy; 421 << 281 422 // Scattered electron (positron) angles. ( Z 282 // Scattered electron (positron) angles. ( Z - axis along the parent photon) 423 283 424 //electron kinematics 284 //electron kinematics 425 G4double electronKineEnergy = std::max(0.,el << 285 G4double costheta_el,costheta_po; 426 G4double costheta_el = G4UniformRand()*2.0-1 << 286 G4double phi_el,phi_po; >> 287 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; >> 288 costheta_el = G4UniformRand()*2.0-1.0; 427 G4double kk = std::sqrt(electronKineEnergy*( 289 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2)); 428 costheta_el = (costheta_el*electronTotEnergy 290 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk); 429 G4double phi_el = twopi * G4UniformRand() ; << 291 phi_el = twopi * G4UniformRand() ; 430 G4double dirX_el = std::sqrt(1.-costheta_el* 292 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el); 431 G4double dirY_el = std::sqrt(1.-costheta_el* 293 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el); 432 G4double dirZ_el = costheta_el; 294 G4double dirZ_el = costheta_el; 433 295 434 //positron kinematics 296 //positron kinematics 435 G4double positronKineEnergy = std::max(0.,po 297 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; 436 G4double costheta_po = G4UniformRand()*2.0-1 << 298 costheta_po = G4UniformRand()*2.0-1.0; 437 kk = std::sqrt(positronKineEnergy*(positronK 299 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2)); 438 costheta_po = (costheta_po*positronTotEnergy 300 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk); 439 G4double phi_po = twopi * G4UniformRand() ; << 301 phi_po = twopi * G4UniformRand() ; 440 G4double dirX_po = std::sqrt(1.-costheta_po* 302 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po); 441 G4double dirY_po = std::sqrt(1.-costheta_po* 303 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po); 442 G4double dirZ_po = costheta_po; 304 G4double dirZ_po = costheta_po; 443 305 444 // Kinematics of the created pair: 306 // Kinematics of the created pair: 445 // the electron and positron are assumed to 307 // the electron and positron are assumed to have a symetric angular 446 // distribution with respect to the Z axis a 308 // distribution with respect to the Z axis along the parent photon 447 G4double localEnergyDeposit = 0. ; 309 G4double localEnergyDeposit = 0. ; 448 310 449 if (electronKineEnergy > 0.0) << 311 //Generate explicitely the electron in the pair, only if it is > threshold >> 312 const G4ProductionCutsTable* theCoupleTable= >> 313 G4ProductionCutsTable::GetProductionCutsTable(); >> 314 size_t indx = couple->GetIndex(); >> 315 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx]; >> 316 //G4double cutP = (*(theCoupleTable->GetEnergyCutsVector(2)))[indx]; >> 317 >> 318 if (electronKineEnergy > cutE) 450 { 319 { 451 G4ThreeVector electronDirection ( dirX_e 320 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el); 452 electronDirection.rotateUz(photonDirecti 321 electronDirection.rotateUz(photonDirection); 453 G4DynamicParticle* electron = new G4Dyna 322 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 454 electronDirection, 323 electronDirection, 455 electronKineEnergy); 324 electronKineEnergy); 456 fvect->push_back(electron); 325 fvect->push_back(electron); 457 } 326 } 458 else 327 else 459 { 328 { 460 localEnergyDeposit += electronKineEnergy 329 localEnergyDeposit += electronKineEnergy; 461 electronKineEnergy = 0; 330 electronKineEnergy = 0; 462 } 331 } 463 332 464 //Generate the positron. Real particle in an 333 //Generate the positron. Real particle in any case, because it will annihilate. If below 465 //threshold, produce it at rest 334 //threshold, produce it at rest 466 if (positronKineEnergy < 0.0) << 335 if (positronKineEnergy < cutE) 467 { 336 { 468 localEnergyDeposit += positronKineEnergy 337 localEnergyDeposit += positronKineEnergy; 469 positronKineEnergy = 0; //produce it at 338 positronKineEnergy = 0; //produce it at rest 470 } 339 } 471 G4ThreeVector positronDirection(dirX_po,dirY 340 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po); 472 positronDirection.rotateUz(photonDirection); 341 positronDirection.rotateUz(photonDirection); 473 G4DynamicParticle* positron = new G4DynamicP 342 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(), 474 positronDirection, positronK 343 positronDirection, positronKineEnergy); 475 fvect->push_back(positron); 344 fvect->push_back(positron); 476 345 477 //Add rest of energy to the local energy dep << 346 //Update the status of the primary gamma (kill it) >> 347 fParticleChange->SetProposedKineticEnergy(0.); 478 fParticleChange->ProposeLocalEnergyDeposit(l 348 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); >> 349 fParticleChange->ProposeTrackStatus(fStopAndKill); 479 350 480 if (fVerboseLevel > 1) << 351 >> 352 if (verboseLevel > 1) 481 { 353 { 482 G4cout << "----------------------------- 354 G4cout << "-----------------------------------------------------------" << G4endl; 483 G4cout << "Energy balance from G4Penelop 355 G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl; 484 G4cout << "Incoming photon energy: " << 356 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl; 485 G4cout << "----------------------------- 357 G4cout << "-----------------------------------------------------------" << G4endl; 486 if (electronKineEnergy) 358 if (electronKineEnergy) 487 G4cout << "Electron (explicitly produced) " << 359 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" << G4endl; 488 << G4endl; << 489 if (positronKineEnergy) 360 if (positronKineEnergy) 490 G4cout << "Positron (not at rest) " << posit 361 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl; 491 G4cout << "Rest masses of e+/- " << 2.0* 362 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl; 492 if (localEnergyDeposit) 363 if (localEnergyDeposit) 493 G4cout << "Local energy deposit " << localEn 364 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 494 G4cout << "Total final state: " << (elec 365 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+ 495 localEnergyDeposit+2.0*electron_ma 366 localEnergyDeposit+2.0*electron_mass_c2)/keV << 496 " keV" << G4endl; 367 " keV" << G4endl; 497 G4cout << "----------------------------- 368 G4cout << "-----------------------------------------------------------" << G4endl; 498 } 369 } 499 if (fVerboseLevel > 0) << 370 if (verboseLevel > 0) 500 { 371 { 501 G4double energyDiff = std::fabs(electron 372 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+ 502 localEnergyDeposit+2.0*electron_ 373 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy); 503 if (energyDiff > 0.05*keV) 374 if (energyDiff > 0.05*keV) 504 G4cout << "Warning from G4PenelopeGammaConve << 375 G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " << 505 << (electronKineEnergy+positronKineEn << 376 (electronKineEnergy+positronKineEnergy+ 506 localEnergyDeposit+2.0*electron_mass_c2 << 377 localEnergyDeposit+2.0*electron_mass_c2)/keV << " keV (final) vs. " << 507 << " keV (final) vs. " << photonEnerg << 378 photonEnergy/keV << " keV (initial)" << G4endl; 508 } 379 } >> 380 509 } 381 } 510 382 511 //....oooOO0OOooo........oooOO0OOooo........oo 383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 512 384 513 void G4PenelopeGammaConversionModel::ReadDataF << 385 std::vector<G4double> G4PenelopeGammaConversionModel::ScreenFunction(G4double b) 514 { 386 { 515 if (!IsMaster()) << 387 std::vector<G4double> result; 516 //Should not be here! << 388 result.clear(); 517 G4Exception("G4PenelopeGammaConversionMode << 389 G4double bsquare=b*b; 518 "em0100",FatalException,"Worker thread in << 390 G4double a0,f1,f2; 519 << 391 f1=2.0-2*std::log(1+bsquare); 520 if (fVerboseLevel > 2) << 392 f2=f1-(2.0/3.0); 521 { << 393 if (b < 1.0e-10) 522 G4cout << "G4PenelopeGammaConversionMode << 523 G4cout << "Going to read Gamma Conversio << 524 } << 525 << 526 const char* path = G4FindDataDir("G4LEDATA << 527 if(!path) << 528 { 394 { 529 G4String excep = << 395 f1=f1-twopi*b; 530 "G4PenelopeGammaConversionModel - G4LEDATA e << 531 G4Exception("G4PenelopeGammaConversionMo << 532 "em0006",FatalException,excep); << 533 return; << 534 } 396 } 535 << 536 /* << 537 Read the cross section file << 538 */ << 539 std::ostringstream ost; << 540 if (Z>9) << 541 ost << path << "/penelope/pairproduction/p << 542 else 397 else 543 ost << path << "/penelope/pairproduction/p << 398 { 544 std::ifstream file(ost.str().c_str()); << 399 a0 = 4*b*std::atan(1.0/b); 545 if (!file.is_open()) << 400 f1 = f1 - a0; 546 { << 401 f2 = f2+2*bsquare*(4.0-a0-3*std::log((1+bsquare)/bsquare)); 547 G4String excep = "G4PenelopeGammaConvers << 548 G4String(ost.str()) + " not found!"; << 549 G4Exception("G4PenelopeGammaConversionMo << 550 "em0003",FatalException,excep); << 551 } << 552 << 553 //I have to know in advance how many points << 554 //to initialize the G4PhysicsFreeVector() << 555 std::size_t ndata=0; << 556 G4String line; << 557 while( getline(file, line) ) << 558 ndata++; << 559 ndata -= 1; //remove one header line << 560 << 561 file.clear(); << 562 file.close(); << 563 file.open(ost.str().c_str()); << 564 G4int readZ =0; << 565 file >> readZ; << 566 << 567 if (fVerboseLevel > 3) << 568 G4cout << "Element Z=" << Z << G4endl; << 569 << 570 //check the right file is opened. << 571 if (readZ != Z) << 572 { << 573 G4ExceptionDescription ed; << 574 ed << "Corrupted data file for Z=" << Z << 575 G4Exception("G4PenelopeGammaConversionMo << 576 "em0005",FatalException,ed); << 577 } << 578 << 579 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 580 G4double ene=0,xs=0; << 581 for (std::size_t i=0;i<ndata;++i) << 582 { << 583 file >> ene >> xs; << 584 //dimensional quantities << 585 ene *= eV; << 586 xs *= barn; << 587 if (xs < 1e-40*cm2) //protection against << 588 xs = 1e-40*cm2; << 589 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 590 } 402 } 591 file.close(); << 403 result.push_back(0.5*(3*f1-f2)); 592 << 404 result.push_back(0.25*(3*f1+f2)); 593 return; << 405 return result; 594 } 406 } 595 407 596 //....oooOO0OOooo........oooOO0OOooo........oo 408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 597 409 598 void G4PenelopeGammaConversionModel::Initializ << 410 G4double G4PenelopeGammaConversionModel::GetScreeningRadius(G4double Z) 599 { 411 { 600 // This is subroutine GPPa0 of Penelope << 412 G4double result = 0; 601 // << 413 G4bool foundElement = false; 602 // 1) calculate the effective Z for the purp << 414 G4int iZ = (G4int) Z; 603 // << 415 if (!fTheScreeningRadii) 604 G4double zeff = 0; << 416 fTheScreeningRadii = new std::map<G4int,G4double>; 605 G4int intZ = 0; << 417 606 G4int nElements = (G4int)material->GetNumber << 418 if (fTheScreeningRadii->count(iZ)) 607 const G4ElementVector* elementVector = mater << 419 { 608 << 420 //The element is already loaded: just return it 609 //avoid calculations if only one building el << 421 result = fTheScreeningRadii->find(iZ)->second; 610 if (nElements == 1) << 422 return result; 611 { << 612 zeff = (*elementVector)[0]->GetZ(); << 613 intZ = (G4int) zeff; << 614 } 423 } 615 else // many elements...let's do the calcula << 424 else //retrieve all from file 616 { 425 { 617 const G4double* fractionVector = materia << 426 char* path = getenv("G4LEDATA"); 618 << 427 if (!path) 619 G4double atot = 0; << 620 for (G4int i=0;i<nElements;i++) << 621 { 428 { 622 G4double Zelement = (*elementVector)[i]->G << 429 G4String excep = "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!"; 623 G4double Aelement = (*elementVector)[i]->G << 430 G4Exception(excep); 624 atot += Aelement*fractionVector[i]; << 431 } 625 zeff += Zelement*Aelement*fractionVector[i << 432 G4String pathString(path); >> 433 G4String pathFile = pathString + "/penelope/pp-pen.dat"; >> 434 std::ifstream file(pathFile); >> 435 >> 436 if (!(file.is_open())) >> 437 { >> 438 G4String excep = "G4PenelopeGammaConversionModel - data file " + pathFile + "not found!"; >> 439 G4Exception(excep); >> 440 } >> 441 G4int k; >> 442 G4double a1,a2; >> 443 while(!file.eof()) { >> 444 file >> k >> a1 >> a2; >> 445 fTheScreeningRadii->insert(std::make_pair(k,a1)); >> 446 if ((G4double) k == Z) >> 447 { >> 448 result = a1; >> 449 foundElement = true; >> 450 } >> 451 } >> 452 file.close(); >> 453 if (verboseLevel > 2) >> 454 G4cout << "Read file pp-pen.dat" << G4endl; >> 455 if (foundElement) >> 456 return result; >> 457 else >> 458 { >> 459 G4String excep = "G4PenelopeGammaConversionModel - Screening Radius for not found in the data file"; >> 460 G4Exception(excep); >> 461 return 0; 626 } 462 } 627 atot /= material->GetTotNbOfAtomsPerVolu << 628 zeff /= (material->GetTotNbOfAtomsPerVol << 629 << 630 intZ = (G4int) (zeff+0.25); << 631 if (intZ <= 0) << 632 intZ = 1; << 633 if (intZ > fMaxZ) << 634 intZ = fMaxZ; << 635 } << 636 << 637 if (fEffectiveCharge) << 638 fEffectiveCharge->insert(std::make_pair(ma << 639 << 640 // << 641 // 2) Calculate Coulomb Correction << 642 // << 643 G4double alz = fine_structure_const*zeff; << 644 G4double alzSquared = alz*alz; << 645 G4double fc = alzSquared*(0.202059-alzSquar << 646 (0.03693-alzSquared* << 647 (0.00835-alzSquared*(0.00201-alzSq << 648 (0.00049-alzSquared* << 649 (0.00012-alzSquared*0.00003))) << 650 +1.0/(alzSquared+1.0)); << 651 // << 652 // 3) Screening functions and low-energy cor << 653 // << 654 G4double matRadius = 2.0/ fAtomicScreeningRa << 655 if (fMaterialInvScreeningRadius) << 656 fMaterialInvScreeningRadius->insert(std::m << 657 << 658 std::pair<G4double,G4double> myPair(0,0); << 659 G4double f0a = 4.0*G4Log(fAtomicScreeningRad << 660 G4double f0b = f0a - 4.0*fc; << 661 myPair.first = f0a; << 662 myPair.second = f0b; << 663 << 664 if (fScreeningFunction) << 665 fScreeningFunction->insert(std::make_pair( << 666 << 667 if (fVerboseLevel > 2) << 668 { << 669 G4cout << "Average Z for material " << m << 670 zeff << G4endl; << 671 G4cout << "Effective radius for material << 672 fAtomicScreeningRadius[intZ] << " m_e*c/hbar << 673 matRadius << G4endl; << 674 G4cout << "Screening parameters F0 for m << 675 f0a << "," << f0b << G4endl; << 676 } 463 } 677 return; << 678 } 464 } 679 465 680 //....oooOO0OOooo........oooOO0OOooo........oo 466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 681 467 682 std::pair<G4double,G4double> << 468 G4double G4PenelopeGammaConversionModel::CoulombCorrection(G4double a) 683 G4PenelopeGammaConversionModel::GetScreeningFu << 684 { 469 { 685 // This is subroutine SCHIFF of Penelope << 470 G4double fc=0; 686 // << 471 G4double b[7] = {0.202059,-0.03693,0.00835,-0.00201,0.00049,-0.00012,0.00003}; 687 // Screening Functions F1(B) and F2(B) in th << 472 G4double aSquared = a*a; 688 // section for pair production << 473 G4double aFourth = aSquared*aSquared; 689 // << 474 G4double aEighth = aFourth*aFourth; 690 std::pair<G4double,G4double> result(0.,0.); << 475 691 G4double BSquared = B*B; << 476 fc = ((1.0/(1.0+a*a))+b[0]+b[1]*aSquared+b[2]*aFourth+b[3]*(aSquared*aFourth)+ 692 G4double f1 = 2.0-2.0*G4Log(1.0+BSquared); << 477 b[4]*aEighth+b[5]*(aEighth*aSquared)+b[6]*(aEighth*aFourth)); 693 G4double f2 = f1 - 6.66666666e-1; // (-2/3) << 478 fc=aSquared*fc; 694 if (B < 1.0e-10) << 479 return fc; 695 f1 = f1-twopi*B; << 696 else << 697 { << 698 G4double a0 = 4.0*B*std::atan(1./B); << 699 f1 = f1 - a0; << 700 f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1. << 701 } << 702 G4double g1 = 0.5*(3.0*f1-f2); << 703 G4double g2 = 0.25*(3.0*f1+f2); << 704 << 705 result.first = g1; << 706 result.second = g2; << 707 << 708 return result; << 709 } 480 } 710 481 711 //....oooOO0OOooo........oooOO0OOooo........oo << 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 712 483 713 void G4PenelopeGammaConversionModel::SetPartic << 484 G4double G4PenelopeGammaConversionModel::LowEnergyCorrection(G4double a,G4double eki) 714 { 485 { 715 if(!fParticle) { << 486 G4double f0=0,t=0; 716 fParticle = p; << 487 G4double b[12] = {-1.744,-12.10,11.18,8.523,73.26,-41.41,-13.52,-121.1,94.41,8.946,62.05,-63.41}; 717 } << 488 t=std::sqrt(2.0*eki); >> 489 G4double tSq = t*t; >> 490 f0=(b[0]+b[1]*a+b[2]*a*a)*t+(b[3]+b[4]*a+b[5]*a*a)*(tSq)+(b[6]+b[7]*a+b[8]*a*a)*(tSq*t)+ >> 491 (b[9]+b[10]*a+b[11]*a*a)*(tSq*tSq); >> 492 return f0; >> 493 718 } 494 } 719 495