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 // 26 // 27 // 27 // 28 // 28 // 29 ////////////////////////////////////////////// 29 //////////////////////////////////////////////////////////////////////// 30 // Optical Photon Rayleigh Scattering Class Im 30 // Optical Photon Rayleigh Scattering Class Implementation 31 ////////////////////////////////////////////// 31 //////////////////////////////////////////////////////////////////////// 32 // 32 // 33 // File: G4OpRayleigh.cc 33 // File: G4OpRayleigh.cc 34 // Description: Discrete Process -- Rayleigh s 34 // Description: Discrete Process -- Rayleigh scattering of optical 35 // photons 35 // photons 36 // Version: 1.0 36 // Version: 1.0 37 // Created: 1996-05-31 37 // Created: 1996-05-31 38 // Author: Juliet Armstrong 38 // Author: Juliet Armstrong 39 // Updated: 2014-10-10 - This version cal 39 // Updated: 2014-10-10 - This version calculates the Rayleigh scattering 40 // length for more materials than 40 // length for more materials than just Water (although the Water 41 // default is kept). To do this t 41 // default is kept). To do this the user would need to specify the 42 // ISOTHERMAL_COMPRESSIBILITY as 42 // ISOTHERMAL_COMPRESSIBILITY as a material property and 43 // optionally an RS_SCALE_LENGTH 43 // optionally an RS_SCALE_LENGTH (useful for testing). Code comes 44 // from Philip Graham (Queen Mary 44 // from Philip Graham (Queen Mary University of London). 45 // 2010-06-11 - Fix Bug 207; Than 45 // 2010-06-11 - Fix Bug 207; Thanks to Xin Qian 46 // (Kellogg Radiation Lab of Calt 46 // (Kellogg Radiation Lab of Caltech) 47 // 2005-07-28 - add G4ProcessType 47 // 2005-07-28 - add G4ProcessType to constructor 48 // 2001-10-18 by Peter Gumplinger 48 // 2001-10-18 by Peter Gumplinger 49 // eliminate unused variable warn 49 // eliminate unused variable warning on Linux (gcc-2.95.2) 50 // 2001-09-18 by mma 50 // 2001-09-18 by mma 51 // >numOfMaterials=G4Material::Ge 51 // >numOfMaterials=G4Material::GetNumberOfMaterials() in BuildPhy 52 // 2001-01-30 by Peter Gumplinger 52 // 2001-01-30 by Peter Gumplinger 53 // > allow for positiv and negati 53 // > allow for positiv and negative CosTheta and force the 54 // > new momentum direction to be 54 // > new momentum direction to be in the same plane as the 55 // > new and old polarization vec 55 // > new and old polarization vectors 56 // 2001-01-29 by Peter Gumplinger 56 // 2001-01-29 by Peter Gumplinger 57 // > fix calculation of SinTheta 57 // > fix calculation of SinTheta (from CosTheta) 58 // 1997-04-09 by Peter Gumplinger 58 // 1997-04-09 by Peter Gumplinger 59 // > new physics/tracking scheme 59 // > new physics/tracking scheme 60 // 60 // 61 ////////////////////////////////////////////// 61 //////////////////////////////////////////////////////////////////////// 62 62 63 #include "G4OpRayleigh.hh" 63 #include "G4OpRayleigh.hh" 64 #include "G4ios.hh" 64 #include "G4ios.hh" 65 #include "G4PhysicalConstants.hh" 65 #include "G4PhysicalConstants.hh" 66 #include "G4SystemOfUnits.hh" 66 #include "G4SystemOfUnits.hh" 67 #include "G4OpticalParameters.hh" 67 #include "G4OpticalParameters.hh" 68 #include "G4OpProcessSubType.hh" 68 #include "G4OpProcessSubType.hh" 69 69 70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 71 G4OpRayleigh::G4OpRayleigh(const G4String& pro 71 G4OpRayleigh::G4OpRayleigh(const G4String& processName, G4ProcessType type) 72 : G4VDiscreteProcess(processName, type) 72 : G4VDiscreteProcess(processName, type) 73 { 73 { 74 Initialise(); 74 Initialise(); 75 SetProcessSubType(fOpRayleigh); 75 SetProcessSubType(fOpRayleigh); 76 thePhysicsTable = nullptr; 76 thePhysicsTable = nullptr; 77 77 78 if(verboseLevel > 0) 78 if(verboseLevel > 0) 79 { 79 { 80 G4cout << GetProcessName() << " is created 80 G4cout << GetProcessName() << " is created " << G4endl; 81 } 81 } 82 } 82 } 83 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 85 G4OpRayleigh::~G4OpRayleigh() 85 G4OpRayleigh::~G4OpRayleigh() 86 { 86 { 87 // VI: inside this PhysicsTable all properti 87 // VI: inside this PhysicsTable all properties are unique 88 // it is not possible to destroy 88 // it is not possible to destroy 89 delete thePhysicsTable; << 89 if(thePhysicsTable) >> 90 { >> 91 delete thePhysicsTable; >> 92 } 90 } 93 } 91 94 92 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 93 void G4OpRayleigh::PreparePhysicsTable(const G 96 void G4OpRayleigh::PreparePhysicsTable(const G4ParticleDefinition&) 94 { 97 { 95 Initialise(); 98 Initialise(); 96 } 99 } 97 100 98 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 99 void G4OpRayleigh::Initialise() 102 void G4OpRayleigh::Initialise() 100 { 103 { 101 SetVerboseLevel(G4OpticalParameters::Instanc 104 SetVerboseLevel(G4OpticalParameters::Instance()->GetRayleighVerboseLevel()); 102 } 105 } 103 106 104 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( 108 G4VParticleChange* G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, 106 109 const G4Step& aStep) 107 { 110 { 108 aParticleChange.Initialize(aTrack); 111 aParticleChange.Initialize(aTrack); 109 const G4DynamicParticle* aParticle = aTrack. 112 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 110 113 111 if(verboseLevel > 1) 114 if(verboseLevel > 1) 112 { 115 { 113 G4cout << "OpRayleigh: Scattering Photon!" 116 G4cout << "OpRayleigh: Scattering Photon!" << G4endl 114 << "Old Momentum Direction: " << aP 117 << "Old Momentum Direction: " << aParticle->GetMomentumDirection() 115 << G4endl << "Old Polarization: " < 118 << G4endl << "Old Polarization: " << aParticle->GetPolarization() 116 << G4endl; 119 << G4endl; 117 } 120 } 118 121 119 G4double cosTheta; 122 G4double cosTheta; 120 G4ThreeVector oldMomDir, newMomDir; 123 G4ThreeVector oldMomDir, newMomDir; 121 G4ThreeVector oldPol, newPol; 124 G4ThreeVector oldPol, newPol; 122 G4double rand; 125 G4double rand; 123 G4double cost, sint, sinphi, cosphi; 126 G4double cost, sint, sinphi, cosphi; 124 127 125 do 128 do 126 { 129 { 127 // Try to simulate the scattered photon mo 130 // Try to simulate the scattered photon momentum direction 128 // w.r.t. the initial photon momentum dire 131 // w.r.t. the initial photon momentum direction 129 cost = G4UniformRand(); 132 cost = G4UniformRand(); 130 sint = std::sqrt(1. - cost * cost); 133 sint = std::sqrt(1. - cost * cost); 131 // consider for the angle 90-180 degrees 134 // consider for the angle 90-180 degrees 132 if(G4UniformRand() < 0.5) 135 if(G4UniformRand() < 0.5) 133 cost = -cost; 136 cost = -cost; 134 137 135 // simulate the phi angle 138 // simulate the phi angle 136 rand = twopi * G4UniformRand(); 139 rand = twopi * G4UniformRand(); 137 sinphi = std::sin(rand); 140 sinphi = std::sin(rand); 138 cosphi = std::cos(rand); 141 cosphi = std::cos(rand); 139 142 140 // construct the new momentum direction 143 // construct the new momentum direction 141 newMomDir.set(sint * cosphi, sint * sinphi 144 newMomDir.set(sint * cosphi, sint * sinphi, cost); 142 oldMomDir = aParticle->GetMomentumDirectio 145 oldMomDir = aParticle->GetMomentumDirection(); 143 newMomDir.rotateUz(oldMomDir); 146 newMomDir.rotateUz(oldMomDir); 144 147 145 // calculate the new polarization directio 148 // calculate the new polarization direction 146 // The new polarization needs to be in the 149 // The new polarization needs to be in the same plane as the new 147 // momentum direction and the old polariza 150 // momentum direction and the old polarization direction 148 oldPol = aParticle->GetPolarization(); 151 oldPol = aParticle->GetPolarization(); 149 newPol = (oldPol - newMomDir.dot(oldPol) * 152 newPol = (oldPol - newMomDir.dot(oldPol) * newMomDir).unit(); 150 153 151 // There is a corner case, where the new m 154 // There is a corner case, where the new momentum direction 152 // is the same as old polarization directi 155 // is the same as old polarization direction: 153 // random generate the azimuthal angle w.r 156 // random generate the azimuthal angle w.r.t. new momentum direction 154 if(newPol.mag() == 0.) 157 if(newPol.mag() == 0.) 155 { 158 { 156 rand = G4UniformRand() * twopi; 159 rand = G4UniformRand() * twopi; 157 newPol.set(std::cos(rand), std::sin(rand 160 newPol.set(std::cos(rand), std::sin(rand), 0.); 158 newPol.rotateUz(newMomDir); 161 newPol.rotateUz(newMomDir); 159 } 162 } 160 else 163 else 161 { 164 { 162 // There are two directions perpendicula 165 // There are two directions perpendicular to the new momentum direction 163 if(G4UniformRand() < 0.5) 166 if(G4UniformRand() < 0.5) 164 newPol = -newPol; 167 newPol = -newPol; 165 } 168 } 166 169 167 // simulate according to the distribution 170 // simulate according to the distribution cos^2(theta) 168 cosTheta = newPol.dot(oldPol); 171 cosTheta = newPol.dot(oldPol); 169 // Loop checking, 13-Aug-2015, Peter Gumpl 172 // Loop checking, 13-Aug-2015, Peter Gumplinger 170 } while(std::pow(cosTheta, 2) < G4UniformRan 173 } while(std::pow(cosTheta, 2) < G4UniformRand()); 171 174 172 aParticleChange.ProposePolarization(newPol); 175 aParticleChange.ProposePolarization(newPol); 173 aParticleChange.ProposeMomentumDirection(new 176 aParticleChange.ProposeMomentumDirection(newMomDir); 174 177 175 if(verboseLevel > 1) 178 if(verboseLevel > 1) 176 { 179 { 177 G4cout << "New Polarization: " << newPol < 180 G4cout << "New Polarization: " << newPol << G4endl 178 << "Polarization Change: " << *(aPa 181 << "Polarization Change: " << *(aParticleChange.GetPolarization()) 179 << G4endl << "New Momentum Directio 182 << G4endl << "New Momentum Direction: " << newMomDir << G4endl 180 << "Momentum Change: " << *(aPartic 183 << "Momentum Change: " << *(aParticleChange.GetMomentumDirection()) 181 << G4endl; 184 << G4endl; 182 } 185 } 183 186 184 return G4VDiscreteProcess::PostStepDoIt(aTra 187 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 185 } 188 } 186 189 187 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 188 void G4OpRayleigh::BuildPhysicsTable(const G4P 191 void G4OpRayleigh::BuildPhysicsTable(const G4ParticleDefinition&) 189 { 192 { 190 if(thePhysicsTable) 193 if(thePhysicsTable) 191 { 194 { 192 // thePhysicsTable->clearAndDestroy(); 195 // thePhysicsTable->clearAndDestroy(); 193 delete thePhysicsTable; 196 delete thePhysicsTable; 194 thePhysicsTable = nullptr; 197 thePhysicsTable = nullptr; 195 } 198 } 196 199 197 const G4MaterialTable* theMaterialTable = G4 200 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 198 const size_t numOfMaterials = G4 201 const size_t numOfMaterials = G4Material::GetNumberOfMaterials(); 199 thePhysicsTable = ne 202 thePhysicsTable = new G4PhysicsTable(numOfMaterials); 200 203 201 for(size_t i = 0; i < numOfMaterials; ++i) 204 for(size_t i = 0; i < numOfMaterials; ++i) 202 { 205 { 203 G4Material* material = (*the 206 G4Material* material = (*theMaterialTable)[i]; 204 G4MaterialPropertiesTable* matProp = mater 207 G4MaterialPropertiesTable* matProp = material->GetMaterialPropertiesTable(); 205 G4PhysicsFreeVector* rayleigh = nullptr; 208 G4PhysicsFreeVector* rayleigh = nullptr; 206 if(matProp) 209 if(matProp) 207 { 210 { 208 rayleigh = matProp->GetProperty(kRAYLEIG 211 rayleigh = matProp->GetProperty(kRAYLEIGH); 209 if(rayleigh == nullptr) 212 if(rayleigh == nullptr) 210 rayleigh = CalculateRayleighMeanFreePa 213 rayleigh = CalculateRayleighMeanFreePaths(material); 211 } 214 } 212 thePhysicsTable->insertAt(i, rayleigh); 215 thePhysicsTable->insertAt(i, rayleigh); 213 } 216 } 214 } 217 } 215 218 216 //....oooOO0OOooo........oooOO0OOooo........oo 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 217 G4double G4OpRayleigh::GetMeanFreePath(const G 220 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack, G4double, 218 G4Force 221 G4ForceCondition*) 219 { 222 { 220 auto rayleigh = static_cast<G4PhysicsFreeVec << 223 G4PhysicsFreeVector* rayleigh = >> 224 static_cast<G4PhysicsFreeVector*>( 221 (*thePhysicsTable)(aTrack.GetMaterial()- 225 (*thePhysicsTable)(aTrack.GetMaterial()->GetIndex())); 222 226 223 G4double rsLength = DBL_MAX; 227 G4double rsLength = DBL_MAX; 224 if(rayleigh) 228 if(rayleigh) 225 { 229 { 226 rsLength = rayleigh->Value(aTrack.GetDynam 230 rsLength = rayleigh->Value(aTrack.GetDynamicParticle()->GetTotalMomentum(), 227 idx_rslength); 231 idx_rslength); 228 } 232 } 229 return rsLength; 233 return rsLength; 230 } 234 } 231 235 232 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa 237 G4PhysicsFreeVector* G4OpRayleigh::CalculateRayleighMeanFreePaths( 234 const G4Material* material) const 238 const G4Material* material) const 235 { 239 { 236 G4MaterialPropertiesTable* MPT = material->G 240 G4MaterialPropertiesTable* MPT = material->GetMaterialPropertiesTable(); 237 241 238 // Retrieve the beta_T or isothermal compres 242 // Retrieve the beta_T or isothermal compressibility value. For backwards 239 // compatibility use a constant if the mater 243 // compatibility use a constant if the material is "Water". If the material 240 // doesn't have an ISOTHERMAL_COMPRESSIBILIT 244 // doesn't have an ISOTHERMAL_COMPRESSIBILITY constant then return 241 G4double betat; 245 G4double betat; 242 if(material->GetName() == "Water") 246 if(material->GetName() == "Water") 243 { 247 { 244 betat = 7.658e-23 * m3 / MeV; 248 betat = 7.658e-23 * m3 / MeV; 245 } 249 } 246 else if(MPT->ConstPropertyExists(kISOTHERMAL 250 else if(MPT->ConstPropertyExists(kISOTHERMAL_COMPRESSIBILITY)) 247 { 251 { 248 betat = MPT->GetConstProperty(kISOTHERMAL_ 252 betat = MPT->GetConstProperty(kISOTHERMAL_COMPRESSIBILITY); 249 } 253 } 250 else 254 else 251 { 255 { 252 return nullptr; 256 return nullptr; 253 } 257 } 254 258 255 // If the material doesn't have a RINDEX pro 259 // If the material doesn't have a RINDEX property vector then return 256 G4MaterialPropertyVector* rIndex = MPT->GetP 260 G4MaterialPropertyVector* rIndex = MPT->GetProperty(kRINDEX); 257 if(rIndex == nullptr) 261 if(rIndex == nullptr) 258 return nullptr; 262 return nullptr; 259 263 260 // Retrieve the optional scale factor (scale 264 // Retrieve the optional scale factor (scales the scattering length) 261 G4double scaleFactor = 1.0; 265 G4double scaleFactor = 1.0; 262 if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR 266 if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR)) 263 { 267 { 264 scaleFactor = MPT->GetConstProperty(kRS_SC 268 scaleFactor = MPT->GetConstProperty(kRS_SCALE_FACTOR); 265 } 269 } 266 270 267 // Retrieve the material temperature. For ba 271 // Retrieve the material temperature. For backwards compatibility use a 268 // constant if the material is "Water" 272 // constant if the material is "Water" 269 G4double temperature; 273 G4double temperature; 270 if(material->GetName() == "Water") 274 if(material->GetName() == "Water") 271 { 275 { 272 temperature = 276 temperature = 273 283.15 * kelvin; // Temperature of wate 277 283.15 * kelvin; // Temperature of water is 10 degrees celsius 274 } 278 } 275 else 279 else 276 { 280 { 277 temperature = material->GetTemperature(); 281 temperature = material->GetTemperature(); 278 } 282 } 279 283 280 auto rayleighMFPs = new G4PhysicsFreeVector( << 284 G4PhysicsFreeVector* rayleighMFPs = new G4PhysicsFreeVector(); 281 // This calculates the meanFreePath via the 285 // This calculates the meanFreePath via the Einstein-Smoluchowski formula 282 const G4double c1 = 286 const G4double c1 = 283 scaleFactor * betat * temperature * k_Bolt 287 scaleFactor * betat * temperature * k_Boltzmann / (6.0 * pi); 284 288 285 for(size_t uRIndex = 0; uRIndex < rIndex->Ge 289 for(size_t uRIndex = 0; uRIndex < rIndex->GetVectorLength(); ++uRIndex) 286 { 290 { 287 const G4double energy = rIndex->Ene 291 const G4double energy = rIndex->Energy(uRIndex); 288 const G4double rIndexSquared = (*rIndex)[u 292 const G4double rIndexSquared = (*rIndex)[uRIndex] * (*rIndex)[uRIndex]; 289 const G4double xlambda = h_Planck * 293 const G4double xlambda = h_Planck * c_light / energy; 290 const G4double c2 = std::pow(tw 294 const G4double c2 = std::pow(twopi / xlambda, 4); 291 const G4double c3 = 295 const G4double c3 = 292 std::pow(((rIndexSquared - 1.0) * (rInde 296 std::pow(((rIndexSquared - 1.0) * (rIndexSquared + 2.0) / 3.0), 2); 293 297 294 const G4double meanFreePath = 1.0 / (c1 * 298 const G4double meanFreePath = 1.0 / (c1 * c2 * c3); 295 299 296 if(verboseLevel > 0) 300 if(verboseLevel > 0) 297 { 301 { 298 G4cout << energy << "MeV\t" << meanFreeP 302 G4cout << energy << "MeV\t" << meanFreePath << "mm" << G4endl; 299 } 303 } 300 304 301 rayleighMFPs->InsertValues(energy, meanFre 305 rayleighMFPs->InsertValues(energy, meanFreePath); 302 } 306 } 303 307 304 return rayleighMFPs; 308 return rayleighMFPs; 305 } 309 } 306 310 307 //....oooOO0OOooo........oooOO0OOooo........oo 311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 308 void G4OpRayleigh::SetVerboseLevel(G4int verbo 312 void G4OpRayleigh::SetVerboseLevel(G4int verbose) 309 { 313 { 310 verboseLevel = verbose; 314 verboseLevel = verbose; 311 G4OpticalParameters::Instance()->SetRayleigh 315 G4OpticalParameters::Instance()->SetRayleighVerboseLevel(verboseLevel); 312 } 316 } 313 317