Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. 1 // 3 // 2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work 3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement, 4 // * << 6 // and all its terms. 5 // * The Geant4 software is copyright of th << 6 // * the Geant4 Collaboration. It is provided << 7 // * conditions of the Geant4 Software License << 8 // * LICENSE and available at http://cern.ch/ << 9 // * include a list of copyright holders. << 10 // * << 11 // * Neither the authors of this software syst << 12 // * institutes,nor the agencies providing fin << 13 // * work make any representation or warran << 14 // * regarding this software system or assum << 15 // * use. Please see the license in the file << 16 // * for the full disclaimer and the limitatio << 17 // * << 18 // * This code implementation is the result << 19 // * technical work of the GEANT4 collaboratio << 20 // * By using, copying, modifying or distri << 21 // * any work based on the software) you ag << 22 // * use in resulting scientific publicati << 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* << 25 // << 26 // << 27 // 7 // >> 8 // $Id: G4OpRayleigh.cc,v 1.5 2001/01/31 04:08:12 gum Exp $ >> 9 // GEANT4 tag $Name: geant4-03-01 $ 28 // 10 // >> 11 // 29 ////////////////////////////////////////////// 12 //////////////////////////////////////////////////////////////////////// 30 // Optical Photon Rayleigh Scattering Class Im 13 // Optical Photon Rayleigh Scattering Class Implementation 31 ////////////////////////////////////////////// 14 //////////////////////////////////////////////////////////////////////// 32 // 15 // 33 // File: G4OpRayleigh.cc << 16 // File: G4OpRayleigh.cc 34 // Description: Discrete Process -- Rayleigh s << 17 // Description: Discrete Process -- Rayleigh scattering of optical 35 // photons << 18 // photons 36 // Version: 1.0 19 // Version: 1.0 37 // Created: 1996-05-31 << 20 // Created: 1996-05-31 38 // Author: Juliet Armstrong 21 // Author: Juliet Armstrong 39 // Updated: 2014-10-10 - This version cal << 22 // Updated: 2001-01-30 by Peter Gumplinger 40 // length for more materials than << 41 // default is kept). To do this t << 42 // ISOTHERMAL_COMPRESSIBILITY as << 43 // optionally an RS_SCALE_LENGTH << 44 // from Philip Graham (Queen Mary << 45 // 2010-06-11 - Fix Bug 207; Than << 46 // (Kellogg Radiation Lab of Calt << 47 // 2005-07-28 - add G4ProcessType << 48 // 2001-10-18 by Peter Gumplinger << 49 // eliminate unused variable warn << 50 // 2001-09-18 by mma << 51 // >numOfMaterials=G4Material::Ge << 52 // 2001-01-30 by Peter Gumplinger << 53 // > allow for positiv and negati 23 // > allow for positiv and negative CosTheta and force the 54 // > new momentum direction to be 24 // > new momentum direction to be in the same plane as the 55 // > new and old polarization vec 25 // > new and old polarization vectors 56 // 2001-01-29 by Peter Gumplinger 26 // 2001-01-29 by Peter Gumplinger 57 // > fix calculation of SinTheta 27 // > fix calculation of SinTheta (from CosTheta) 58 // 1997-04-09 by Peter Gumplinger 28 // 1997-04-09 by Peter Gumplinger 59 // > new physics/tracking scheme 29 // > new physics/tracking scheme >> 30 // mail: gum@triumf.ca 60 // 31 // 61 ////////////////////////////////////////////// 32 //////////////////////////////////////////////////////////////////////// 62 33 63 #include "G4OpRayleigh.hh" << 64 #include "G4ios.hh" 34 #include "G4ios.hh" 65 #include "G4PhysicalConstants.hh" << 35 #include "G4OpRayleigh.hh" 66 #include "G4SystemOfUnits.hh" << 36 67 #include "G4OpticalParameters.hh" << 37 ///////////////////////// 68 #include "G4OpProcessSubType.hh" << 38 // Class Implementation 69 << 39 ///////////////////////// 70 //....oooOO0OOooo........oooOO0OOooo........oo << 40 71 G4OpRayleigh::G4OpRayleigh(const G4String& pro << 41 ////////////// 72 : G4VDiscreteProcess(processName, type) << 42 // Operators >> 43 ////////////// >> 44 >> 45 // G4OpRayleigh::operator=(const G4OpRayleigh &right) >> 46 // { >> 47 // } >> 48 >> 49 ///////////////// >> 50 // Constructors >> 51 ///////////////// >> 52 >> 53 G4OpRayleigh::G4OpRayleigh(const G4String& processName) >> 54 : G4VDiscreteProcess(processName) 73 { 55 { 74 Initialise(); << 56 75 SetProcessSubType(fOpRayleigh); << 57 thePhysicsTable = NULL; 76 thePhysicsTable = nullptr; << 58 77 << 59 if (verboseLevel>0) { 78 if(verboseLevel > 0) << 60 G4cout << GetProcessName() << " is created " << G4endl; 79 { << 61 } 80 G4cout << GetProcessName() << " is created << 62 81 } << 63 BuildThePhysicsTable(); 82 } 64 } 83 65 84 //....oooOO0OOooo........oooOO0OOooo........oo << 66 // G4OpRayleigh::G4OpRayleigh(const G4OpRayleigh &right) >> 67 // { >> 68 // } >> 69 >> 70 //////////////// >> 71 // Destructors >> 72 //////////////// >> 73 85 G4OpRayleigh::~G4OpRayleigh() 74 G4OpRayleigh::~G4OpRayleigh() 86 { 75 { 87 // VI: inside this PhysicsTable all properti << 76 if (thePhysicsTable!= NULL) { 88 // it is not possible to destroy << 77 thePhysicsTable->clearAndDestroy(); 89 delete thePhysicsTable; << 78 delete thePhysicsTable; >> 79 } 90 } 80 } 91 81 92 //....oooOO0OOooo........oooOO0OOooo........oo << 82 //////////// 93 void G4OpRayleigh::PreparePhysicsTable(const G << 83 // Methods 94 { << 84 //////////// 95 Initialise(); << 96 } << 97 85 98 //....oooOO0OOooo........oooOO0OOooo........oo << 86 // PostStepDoIt 99 void G4OpRayleigh::Initialise() << 87 // ------------- >> 88 // >> 89 G4VParticleChange* >> 90 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) 100 { 91 { 101 SetVerboseLevel(G4OpticalParameters::Instanc << 92 aParticleChange.Initialize(aTrack); 102 } << 103 93 104 //....oooOO0OOooo........oooOO0OOooo........oo << 94 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 105 G4VParticleChange* G4OpRayleigh::PostStepDoIt( << 95 const G4Material* aMaterial = aTrack.GetMaterial(); 106 << 107 { << 108 aParticleChange.Initialize(aTrack); << 109 const G4DynamicParticle* aParticle = aTrack. << 110 96 111 if(verboseLevel > 1) << 97 if (verboseLevel>0) { 112 { << 98 G4cout << "Scattering Photon!" << G4endl; 113 G4cout << "OpRayleigh: Scattering Photon!" << 99 G4cout << "Old Momentum Direction: " 114 << "Old Momentum Direction: " << aP << 100 << aParticle->GetMomentumDirection() << G4endl; 115 << G4endl << "Old Polarization: " < << 101 G4cout << "Old Polarization: " 116 << G4endl; << 102 << aParticle->GetPolarization() << G4endl; 117 } << 103 } 118 << 119 G4double cosTheta; << 120 G4ThreeVector oldMomDir, newMomDir; << 121 G4ThreeVector oldPol, newPol; << 122 G4double rand; << 123 G4double cost, sint, sinphi, cosphi; << 124 << 125 do << 126 { << 127 // Try to simulate the scattered photon mo << 128 // w.r.t. the initial photon momentum dire << 129 cost = G4UniformRand(); << 130 sint = std::sqrt(1. - cost * cost); << 131 // consider for the angle 90-180 degrees << 132 if(G4UniformRand() < 0.5) << 133 cost = -cost; << 134 << 135 // simulate the phi angle << 136 rand = twopi * G4UniformRand(); << 137 sinphi = std::sin(rand); << 138 cosphi = std::cos(rand); << 139 << 140 // construct the new momentum direction << 141 newMomDir.set(sint * cosphi, sint * sinphi << 142 oldMomDir = aParticle->GetMomentumDirectio << 143 newMomDir.rotateUz(oldMomDir); << 144 << 145 // calculate the new polarization directio << 146 // The new polarization needs to be in the << 147 // momentum direction and the old polariza << 148 oldPol = aParticle->GetPolarization(); << 149 newPol = (oldPol - newMomDir.dot(oldPol) * << 150 << 151 // There is a corner case, where the new m << 152 // is the same as old polarization directi << 153 // random generate the azimuthal angle w.r << 154 if(newPol.mag() == 0.) << 155 { << 156 rand = G4UniformRand() * twopi; << 157 newPol.set(std::cos(rand), std::sin(rand << 158 newPol.rotateUz(newMomDir); << 159 } << 160 else << 161 { << 162 // There are two directions perpendicula << 163 if(G4UniformRand() < 0.5) << 164 newPol = -newPol; << 165 } << 166 << 167 // simulate according to the distribution << 168 cosTheta = newPol.dot(oldPol); << 169 // Loop checking, 13-Aug-2015, Peter Gumpl << 170 } while(std::pow(cosTheta, 2) < G4UniformRan << 171 << 172 aParticleChange.ProposePolarization(newPol); << 173 aParticleChange.ProposeMomentumDirection(new << 174 << 175 if(verboseLevel > 1) << 176 { << 177 G4cout << "New Polarization: " << newPol < << 178 << "Polarization Change: " << *(aPa << 179 << G4endl << "New Momentum Directio << 180 << "Momentum Change: " << *(aPartic << 181 << G4endl; << 182 } << 183 104 184 return G4VDiscreteProcess::PostStepDoIt(aTra << 105 // find polar angle w.r.t. old polarization vector 185 } << 186 106 187 //....oooOO0OOooo........oooOO0OOooo........oo << 107 G4double rand = G4UniformRand(); 188 void G4OpRayleigh::BuildPhysicsTable(const G4P << 108 189 { << 109 G4double CosTheta = pow(rand, 1./3.); 190 if(thePhysicsTable) << 110 G4double SinTheta = sqrt(1.-CosTheta*CosTheta); 191 { << 111 192 // thePhysicsTable->clearAndDestroy(); << 112 if(G4UniformRand() < 0.5)CosTheta = -CosTheta; 193 delete thePhysicsTable; << 113 194 thePhysicsTable = nullptr; << 114 // find azimuthal angle w.r.t old polarization vector 195 } << 115 196 << 116 rand = G4UniformRand(); 197 const G4MaterialTable* theMaterialTable = G4 << 117 198 const size_t numOfMaterials = G4 << 118 G4double Phi = twopi*rand; 199 thePhysicsTable = ne << 119 G4double SinPhi = sin(Phi); 200 << 120 G4double CosPhi = cos(Phi); 201 for(size_t i = 0; i < numOfMaterials; ++i) << 121 202 { << 122 G4double unit_x = SinTheta * CosPhi; 203 G4Material* material = (*the << 123 G4double unit_y = SinTheta * SinPhi; 204 G4MaterialPropertiesTable* matProp = mater << 124 G4double unit_z = CosTheta; 205 G4PhysicsFreeVector* rayleigh = nullptr; << 125 206 if(matProp) << 126 G4ThreeVector NewPolarization (unit_x,unit_y,unit_z); 207 { << 127 208 rayleigh = matProp->GetProperty(kRAYLEIG << 128 // Rotate new polarization direction into global reference system 209 if(rayleigh == nullptr) << 129 210 rayleigh = CalculateRayleighMeanFreePa << 130 G4ThreeVector OldPolarization = aParticle->GetPolarization(); 211 } << 131 OldPolarization = OldPolarization.unit(); 212 thePhysicsTable->insertAt(i, rayleigh); << 132 213 } << 133 NewPolarization.rotateUz(OldPolarization); >> 134 NewPolarization = NewPolarization.unit(); >> 135 >> 136 // -- new momentum direction is normal to the new >> 137 // polarization vector and in the same plane as the >> 138 // old and new polarization vectors -- >> 139 >> 140 G4ThreeVector NewMomentumDirection = >> 141 OldPolarization - NewPolarization * CosTheta; >> 142 >> 143 if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection; >> 144 NewMomentumDirection = NewMomentumDirection.unit(); >> 145 >> 146 aParticleChange.SetPolarizationChange(NewPolarization); >> 147 >> 148 aParticleChange.SetMomentumChange(NewMomentumDirection); >> 149 >> 150 if (verboseLevel>0) { >> 151 G4cout << "New Polarization: " >> 152 << NewPolarization << G4endl; >> 153 G4cout << "Polarization Change: " >> 154 << *(aParticleChange.GetPolarizationChange()) << G4endl; >> 155 G4cout << "New Momentum Direction: " >> 156 << NewMomentumDirection << G4endl; >> 157 G4cout << "Momentum Change: " >> 158 << *(aParticleChange.GetMomentumChange()) << G4endl; >> 159 } >> 160 >> 161 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 214 } 162 } 215 163 216 //....oooOO0OOooo........oooOO0OOooo........oo << 164 // BuildThePhysicsTable for the Rayleigh Scattering process 217 G4double G4OpRayleigh::GetMeanFreePath(const G << 165 // -------------------------------------------------------- 218 G4Force << 166 // >> 167 void G4OpRayleigh::BuildThePhysicsTable() 219 { 168 { 220 auto rayleigh = static_cast<G4PhysicsFreeVec << 169 // Builds a table of scattering lengths for each material 221 (*thePhysicsTable)(aTrack.GetMaterial()- << 170 >> 171 if (thePhysicsTable) return; >> 172 >> 173 const G4MaterialTable* theMaterialTable= >> 174 G4Material::GetMaterialTable(); >> 175 G4int numOfMaterials = theMaterialTable->length(); >> 176 >> 177 // create a new physics table >> 178 >> 179 thePhysicsTable = new G4PhysicsTable(numOfMaterials); >> 180 >> 181 // loop for materials >> 182 >> 183 for (G4int i=0 ; i < numOfMaterials; i++) >> 184 { >> 185 G4PhysicsOrderedFreeVector* ScatteringLengths = >> 186 new G4PhysicsOrderedFreeVector(); 222 187 223 G4double rsLength = DBL_MAX; << 188 if ((*theMaterialTable)[i]->GetName() == "Water") 224 if(rayleigh) << 189 { 225 { << 190 G4MaterialPropertiesTable *MaterialPT = 226 rsLength = rayleigh->Value(aTrack.GetDynam << 191 (*theMaterialTable)[i]->GetMaterialPropertiesTable(); 227 idx_rslength); << 192 // Call utility routine to Generate 228 } << 193 // Rayleigh Scattering Lengths 229 return rsLength; << 194 ScatteringLengths = >> 195 RayleighAttenuationLengthGenerator(MaterialPT); >> 196 } >> 197 >> 198 thePhysicsTable->insertAt(i,ScatteringLengths); >> 199 } 230 } 200 } 231 201 232 //....oooOO0OOooo........oooOO0OOooo........oo << 202 // GetMeanFreePath() 233 G4PhysicsFreeVector* G4OpRayleigh::CalculateRa << 203 // ----------------- 234 const G4Material* material) const << 204 // >> 205 G4double G4OpRayleigh::GetMeanFreePath(const G4Track& aTrack, >> 206 G4double , >> 207 G4ForceCondition* ) 235 { 208 { 236 G4MaterialPropertiesTable* MPT = material->G << 209 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); >> 210 const G4Material* aMaterial = aTrack.GetMaterial(); >> 211 >> 212 G4double thePhotonMomentum = aParticle->GetTotalMomentum(); >> 213 >> 214 G4double AttenuationLength = DBL_MAX; 237 215 238 // Retrieve the beta_T or isothermal compres << 216 if (aMaterial->GetName() == "Water") { 239 // compatibility use a constant if the mater << 240 // doesn't have an ISOTHERMAL_COMPRESSIBILIT << 241 G4double betat; << 242 if(material->GetName() == "Water") << 243 { << 244 betat = 7.658e-23 * m3 / MeV; << 245 } << 246 else if(MPT->ConstPropertyExists(kISOTHERMAL << 247 { << 248 betat = MPT->GetConstProperty(kISOTHERMAL_ << 249 } << 250 else << 251 { << 252 return nullptr; << 253 } << 254 << 255 // If the material doesn't have a RINDEX pro << 256 G4MaterialPropertyVector* rIndex = MPT->GetP << 257 if(rIndex == nullptr) << 258 return nullptr; << 259 << 260 // Retrieve the optional scale factor (scale << 261 G4double scaleFactor = 1.0; << 262 if(MPT->ConstPropertyExists(kRS_SCALE_FACTOR << 263 { << 264 scaleFactor = MPT->GetConstProperty(kRS_SC << 265 } << 266 << 267 // Retrieve the material temperature. For ba << 268 // constant if the material is "Water" << 269 G4double temperature; << 270 if(material->GetName() == "Water") << 271 { << 272 temperature = << 273 283.15 * kelvin; // Temperature of wate << 274 } << 275 else << 276 { << 277 temperature = material->GetTemperature(); << 278 } << 279 << 280 auto rayleighMFPs = new G4PhysicsFreeVector( << 281 // This calculates the meanFreePath via the << 282 const G4double c1 = << 283 scaleFactor * betat * temperature * k_Bolt << 284 << 285 for(size_t uRIndex = 0; uRIndex < rIndex->Ge << 286 { << 287 const G4double energy = rIndex->Ene << 288 const G4double rIndexSquared = (*rIndex)[u << 289 const G4double xlambda = h_Planck * << 290 const G4double c2 = std::pow(tw << 291 const G4double c3 = << 292 std::pow(((rIndexSquared - 1.0) * (rInde << 293 << 294 const G4double meanFreePath = 1.0 / (c1 * << 295 << 296 if(verboseLevel > 0) << 297 { << 298 G4cout << energy << "MeV\t" << meanFreeP << 299 } << 300 217 301 rayleighMFPs->InsertValues(energy, meanFre << 218 G4bool isOutRange; 302 } << 219 >> 220 AttenuationLength = >> 221 (*thePhysicsTable)(aMaterial->GetIndex())-> >> 222 GetValue(thePhotonMomentum, isOutRange); >> 223 } >> 224 else { >> 225 >> 226 G4MaterialPropertiesTable* aMaterialPropertyTable = >> 227 aMaterial->GetMaterialPropertiesTable(); >> 228 >> 229 if(aMaterialPropertyTable){ >> 230 G4MaterialPropertyVector* AttenuationLengthVector = >> 231 aMaterialPropertyTable->GetProperty("RAYLEIGH"); >> 232 if(AttenuationLengthVector){ >> 233 AttenuationLength = AttenuationLengthVector -> >> 234 GetProperty(thePhotonMomentum); >> 235 } >> 236 else{ >> 237 // G4cout << "No Rayleigh scattering length specified" << G4endl; >> 238 } >> 239 } >> 240 else{ >> 241 // G4cout << "No Rayleigh scattering length specified" << G4endl; >> 242 } >> 243 } 303 244 304 return rayleighMFPs; << 245 return AttenuationLength; 305 } 246 } 306 247 307 //....oooOO0OOooo........oooOO0OOooo........oo << 248 // RayleighAttenuationLengthGenerator() 308 void G4OpRayleigh::SetVerboseLevel(G4int verbo << 249 // ------------------------------------ >> 250 // Private method to compute Rayleigh Scattering Lengths (for water) >> 251 // >> 252 G4PhysicsOrderedFreeVector* >> 253 G4OpRayleigh::RayleighAttenuationLengthGenerator(G4MaterialPropertiesTable *aMPT) 309 { 254 { 310 verboseLevel = verbose; << 255 // Physical Constants 311 G4OpticalParameters::Instance()->SetRayleigh << 256 >> 257 // isothermal compressibility of water >> 258 G4double betat = 7.658e-23*m3/MeV; >> 259 >> 260 // K Boltzman >> 261 G4double kboltz = 8.61739e-11*MeV/kelvin; >> 262 >> 263 // Temperature of water is 10 degrees celsius >> 264 // conversion to kelvin: >> 265 // TCelsius = TKelvin - 273.15 => 273.15 + 10 = 283.15 >> 266 G4double temp = 283.15*kelvin; >> 267 >> 268 // Retrieve vectors for refraction index >> 269 // and photon momentum from the material properties table >> 270 >> 271 G4MaterialPropertyVector* Rindex = aMPT->GetProperty("RINDEX"); >> 272 >> 273 G4double refsq; >> 274 G4double e; >> 275 G4double xlambda; >> 276 G4double c1, c2, c3, c4; >> 277 G4double Dist; >> 278 G4double refraction_index; >> 279 >> 280 G4double no_unit = 1.0; >> 281 >> 282 G4PhysicsOrderedFreeVector *RayleighScatteringLengths = >> 283 new G4PhysicsOrderedFreeVector(); >> 284 Rindex->ResetIterator(); >> 285 >> 286 while (++(*Rindex)) { >> 287 >> 288 e = (Rindex->GetPhotonMomentum()); >> 289 >> 290 refraction_index = Rindex->GetProperty(); >> 291 refsq = refraction_index*refraction_index; >> 292 xlambda = h_Planck*c_light/e; >> 293 >> 294 if (verboseLevel>0) { >> 295 G4cout << Rindex->GetPhotonMomentum() << " MeV\t"; >> 296 G4cout << xlambda << " mm\t"; >> 297 } >> 298 >> 299 c1 = 1 / (6.0 * pi); >> 300 c2 = pow((2.0 * pi / xlambda), 4); >> 301 c3 = pow( ( (refsq - 1.0) * (refsq + 2.0) / 3.0 ), 2); >> 302 c4 = betat * temp * kboltz; >> 303 >> 304 Dist = 1.0 / (c1*c2*c3*c4); >> 305 >> 306 if (verboseLevel>0) { >> 307 G4cout << Dist << " mm" << G4endl; >> 308 } >> 309 RayleighScatteringLengths-> >> 310 InsertValues(Rindex->GetPhotonMomentum(), Dist); >> 311 } >> 312 >> 313 return RayleighScatteringLengths; 312 } 314 } 313 315