Changeset 1340 for trunk/source/processes/optical
- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- Location:
- trunk/source/processes/optical
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/optical/History
r1337 r1340 17 17 ---------------------------------------------------------- 18 18 19 14th June 2010 Peter Gumplinger (op-V09-03-03) 19 29th Oct 2010 Peter Gumplinger (op-V09-03-06) 20 Fix minor Coverty Errors 21 22 22th Oct 2010 Peter Gumplinger (op-V09-03-05) 23 Add Mie Scattering of Optical Photons 24 Thanks to Xin Qian (Kellogg Radiation Lab of Caltech) 25 Based on work from Vlasios Vasileiou (University of Maryland) 26 Henyey-Greenstein phase functions 27 28 21th Oct 2010 Peter Gumplinger (op-V09-03-04) 29 Fix Bug 207; Thanks to Xin Qian (Kellogg Radiation Lab of Caltech) 30 31 14th Jun 2010 Peter Gumplinger (op-V09-03-03) 20 32 ProposeLocalEnergyDeposit(thePhotonMomentum) every time the 21 33 photon is fStopAndKill because of NoRINDEX. This will trigger -
trunk/source/processes/optical/include/G4OpProcessSubType.hh
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4OpProcessSubType.hh,v 1. 3 2008/12/18 13:02:46 gunterExp $28 // GEANT4 tag $Name: geant4-09-04-beta-01$27 // $Id: G4OpProcessSubType.hh,v 1.4 2010/10/22 01:01:37 gum Exp $ 28 // GEANT4 tag $Name: op-V09-03-06 $ 29 29 // 30 30 //--------------------------------------------------------------- … … 48 48 fOpBoundary = 32, 49 49 fOpRayleigh = 33, 50 fOpWLS = 34 50 fOpWLS = 34, 51 fOpMieHG = 35 51 52 }; 52 53 -
trunk/source/processes/optical/src/G4OpBoundaryProcess.cc
r1337 r1340 120 120 ->GetSurfaceTolerance(); 121 121 122 iTE = iTM = 0; 123 thePhotonMomentum = 0.; 124 Rindex1 = Rindex2 = cost1 = cost2 = sint1 = sint2 = 0.; 122 125 } 123 126 -
trunk/source/processes/optical/src/G4OpRayleigh.cc
r1337 r1340 25 25 // 26 26 // 27 // $Id: G4OpRayleigh.cc,v 1.1 7 2008/10/24 19:51:12gum Exp $28 // GEANT4 tag $Name: geant4-09-04-beta-01$27 // $Id: G4OpRayleigh.cc,v 1.19 2010/10/29 23:18:35 gum Exp $ 28 // GEANT4 tag $Name: op-V09-03-06 $ 29 29 // 30 30 // … … 39 39 // Created: 1996-05-31 40 40 // Author: Juliet Armstrong 41 // Updated: 2005-07-28 - add G4ProcessType to constructor 41 // Updated: 2010-06-11 - Fix Bug 207; Thanks to Xin Qian 42 // (Kellogg Radiation Lab of Caltech) 43 // 2005-07-28 - add G4ProcessType to constructor 42 44 // 2001-10-18 by Peter Gumplinger 43 45 // eliminate unused variable warning on Linux (gcc-2.95.2) … … 116 118 // ------------- 117 119 // 118 G4VParticleChange* 120 G4VParticleChange* 119 121 G4OpRayleigh::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) 120 122 { … … 124 126 125 127 if (verboseLevel>0) { 126 G4cout << "Scattering Photon!" << G4endl; 127 G4cout << "Old Momentum Direction: " 128 << aParticle->GetMomentumDirection() << G4endl; 129 G4cout << "Old Polarization: " 130 << aParticle->GetPolarization() << G4endl; 131 } 132 133 // find polar angle w.r.t. old polarization vector 134 135 G4double rand = G4UniformRand(); 136 137 G4double CosTheta = std::pow(rand, 1./3.); 138 G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta); 139 140 if(G4UniformRand() < 0.5)CosTheta = -CosTheta; 141 142 // find azimuthal angle w.r.t old polarization vector 143 144 rand = G4UniformRand(); 145 146 G4double Phi = twopi*rand; 147 G4double SinPhi = std::sin(Phi); 148 G4double CosPhi = std::cos(Phi); 149 150 G4double unit_x = SinTheta * CosPhi; 151 G4double unit_y = SinTheta * SinPhi; 152 G4double unit_z = CosTheta; 153 154 G4ThreeVector NewPolarization (unit_x,unit_y,unit_z); 155 156 // Rotate new polarization direction into global reference system 157 158 G4ThreeVector OldPolarization = aParticle->GetPolarization(); 159 OldPolarization = OldPolarization.unit(); 160 161 NewPolarization.rotateUz(OldPolarization); 162 NewPolarization = NewPolarization.unit(); 163 164 // -- new momentum direction is normal to the new 165 // polarization vector and in the same plane as the 166 // old and new polarization vectors -- 167 168 G4ThreeVector NewMomentumDirection = 169 OldPolarization - NewPolarization * CosTheta; 170 171 if(G4UniformRand() < 0.5)NewMomentumDirection = -NewMomentumDirection; 172 NewMomentumDirection = NewMomentumDirection.unit(); 173 174 aParticleChange.ProposePolarization(NewPolarization); 175 176 aParticleChange.ProposeMomentumDirection(NewMomentumDirection); 128 G4cout << "Scattering Photon!" << G4endl; 129 G4cout << "Old Momentum Direction: " 130 << aParticle->GetMomentumDirection() << G4endl; 131 G4cout << "Old Polarization: " 132 << aParticle->GetPolarization() << G4endl; 133 } 134 135 G4double cosTheta; 136 G4ThreeVector OldMomentumDirection, NewMomentumDirection; 137 G4ThreeVector OldPolarization, NewPolarization; 138 139 do { 140 // Try to simulate the scattered photon momentum direction 141 // w.r.t. the initial photon momentum direction 142 143 G4double CosTheta = G4UniformRand(); 144 G4double SinTheta = std::sqrt(1.-CosTheta*CosTheta); 145 // consider for the angle 90-180 degrees 146 if (G4UniformRand() < 0.5) CosTheta = -CosTheta; 147 148 // simulate the phi angle 149 G4double rand = twopi*G4UniformRand(); 150 G4double SinPhi = std::sin(rand); 151 G4double CosPhi = std::cos(rand); 152 153 // start constructing the new momentum direction 154 G4double unit_x = SinTheta * CosPhi; 155 G4double unit_y = SinTheta * SinPhi; 156 G4double unit_z = CosTheta; 157 NewMomentumDirection.set (unit_x,unit_y,unit_z); 158 159 // Rotate the new momentum direction into global reference system 160 OldMomentumDirection = aParticle->GetMomentumDirection(); 161 OldMomentumDirection = OldMomentumDirection.unit(); 162 NewMomentumDirection.rotateUz(OldMomentumDirection); 163 NewMomentumDirection = NewMomentumDirection.unit(); 164 165 // calculate the new polarization direction 166 // The new polarization needs to be in the same plane as the new 167 // momentum direction and the old polarization direction 168 OldPolarization = aParticle->GetPolarization(); 169 G4double constant = -1./NewMomentumDirection.dot(OldPolarization); 170 171 NewPolarization = NewMomentumDirection + constant*OldPolarization; 172 NewPolarization = NewPolarization.unit(); 173 174 // There is a corner case, where the Newmomentum direction 175 // is the same as oldpolariztion direction: 176 // random generate the azimuthal angle w.r.t. Newmomentum direction 177 if (NewPolarization.mag() == 0.) { 178 rand = G4UniformRand()*twopi; 179 NewPolarization.set(std::cos(rand),std::sin(rand),0.); 180 NewPolarization.rotateUz(NewMomentumDirection); 181 } else { 182 // There are two directions which are perpendicular 183 // to the new momentum direction 184 if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization; 185 } 186 187 // simulate according to the distribution cos^2(theta) 188 cosTheta = NewPolarization.dot(OldPolarization); 189 } while (std::pow(cosTheta,2) < G4UniformRand()); 190 191 aParticleChange.ProposePolarization(NewPolarization); 192 aParticleChange.ProposeMomentumDirection(NewMomentumDirection); 177 193 178 194 if (verboseLevel>0) { 179 180 181 182 183 184 185 186 187 195 G4cout << "New Polarization: " 196 << NewPolarization << G4endl; 197 G4cout << "Polarization Change: " 198 << *(aParticleChange.GetPolarization()) << G4endl; 199 G4cout << "New Momentum Direction: " 200 << NewMomentumDirection << G4endl; 201 G4cout << "Momentum Change: " 202 << *(aParticleChange.GetMomentumDirection()) << G4endl; 203 } 188 204 189 205 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); … … 211 227 for (G4int i=0 ; i < numOfMaterials; i++) 212 228 { 213 G4PhysicsOrderedFreeVector* ScatteringLengths = 214 new G4PhysicsOrderedFreeVector(); 229 G4PhysicsOrderedFreeVector* ScatteringLengths = NULL; 215 230 216 231 G4MaterialPropertiesTable *aMaterialPropertiesTable = … … 231 246 DefaultWater = true; 232 247 233 248 ScatteringLengths = 234 249 RayleighAttenuationLengthGenerator(aMaterialPropertiesTable); 235 250 }
Note: See TracChangeset
for help on using the changeset viewer.