Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 //////////////////////////////////////////////////////////////////////// 28 // 29 // File G4OpMieHG.hh 30 // Description: Discrete Process -- Mie Scattering of Optical Photons 31 // Created: 2010-07-03 32 // Author: Xin Qian 33 // Based on work from Vlasios Vasileiou 34 // 35 // This subroutine will mimic the Mie scattering based on 36 // Henyey-Greenstein phase function 37 // Forward and backward angles are treated separately. 38 // 39 //////////////////////////////////////////////////////////////////////// 40 41 #include "G4OpMieHG.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4OpticalParameters.hh" 44 #include "G4OpProcessSubType.hh" 45 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 G4OpMieHG::G4OpMieHG(const G4String& processName, G4ProcessType type) 48 : G4VDiscreteProcess(processName, type) 49 { 50 Initialise(); 51 if(verboseLevel > 0) 52 { 53 G4cout << GetProcessName() << " is created " << G4endl; 54 } 55 SetProcessSubType(fOpMieHG); 56 } 57 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 59 G4OpMieHG::~G4OpMieHG() = default; 60 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 62 void G4OpMieHG::PreparePhysicsTable(const G4ParticleDefinition&) 63 { 64 Initialise(); 65 } 66 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 void G4OpMieHG::Initialise() 69 { 70 SetVerboseLevel(G4OpticalParameters::Instance()->GetMieVerboseLevel()); 71 } 72 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 74 G4VParticleChange* G4OpMieHG::PostStepDoIt(const G4Track& aTrack, 75 const G4Step& aStep) 76 { 77 aParticleChange.Initialize(aTrack); 78 79 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 80 const G4MaterialPropertiesTable* MPT = 81 aTrack.GetMaterial()->GetMaterialPropertiesTable(); 82 83 G4double forwardRatio = MPT->GetConstProperty(kMIEHG_FORWARD_RATIO); 84 85 if(verboseLevel > 1) 86 { 87 G4cout << "OpMie Scattering Photon!" << G4endl 88 << " Old Momentum Direction: " << aParticle->GetMomentumDirection() 89 << G4endl 90 << " MIE Old Polarization: " << aParticle->GetPolarization() 91 << G4endl; 92 } 93 94 G4double gg; 95 G4int direction; 96 if(G4UniformRand() <= forwardRatio) 97 { 98 gg = MPT->GetConstProperty(kMIEHG_FORWARD); 99 direction = 1; 100 } 101 else 102 { 103 gg = MPT->GetConstProperty(kMIEHG_BACKWARD); 104 direction = -1; 105 } 106 107 G4double r = G4UniformRand(); 108 109 // sample the direction 110 G4double theta; 111 if(gg != 0.) 112 { 113 theta = std::acos(2. * r * (1. + gg) * (1. + gg) * (1. - gg + gg * r) / 114 ((1. - gg + 2. * gg * r) * (1. - gg + 2. * gg * r)) - 115 1.); 116 } 117 else 118 { 119 theta = std::acos(2. * r - 1.); 120 } 121 G4double phi = G4UniformRand() * twopi; 122 123 if(direction == -1) 124 theta = pi - theta; // backward scattering 125 126 G4ThreeVector newMomDir, oldMomDir; 127 G4ThreeVector newPol, oldPol; 128 129 G4double sinth = std::sin(theta); 130 newMomDir.set(sinth * std::cos(phi), sinth * std::sin(phi), std::cos(theta)); 131 oldMomDir = aParticle->GetMomentumDirection(); 132 newMomDir.rotateUz(oldMomDir); 133 newMomDir = newMomDir.unit(); 134 135 oldPol = aParticle->GetPolarization(); 136 newPol = newMomDir - oldPol / newMomDir.dot(oldPol); 137 newPol = newPol.unit(); 138 139 if(newPol.mag() == 0.) 140 { 141 r = G4UniformRand() * twopi; 142 newPol.set(std::cos(r), std::sin(r), 0.); 143 newPol.rotateUz(newMomDir); 144 } 145 else 146 { 147 // There are two directions perpendicular to new momentum direction 148 if(G4UniformRand() < 0.5) 149 newPol = -newPol; 150 } 151 152 aParticleChange.ProposePolarization(newPol); 153 aParticleChange.ProposeMomentumDirection(newMomDir); 154 155 if(verboseLevel > 1) 156 { 157 G4cout << "OpMie New Polarization: " << newPol << G4endl 158 << " Polarization Change: " << *(aParticleChange.GetPolarization()) 159 << G4endl << " New Momentum Direction: " << newMomDir << G4endl 160 << " Momentum Change: " << *(aParticleChange.GetMomentumDirection()) 161 << G4endl; 162 } 163 164 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 165 } 166 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 168 G4double G4OpMieHG::GetMeanFreePath(const G4Track& aTrack, G4double, 169 G4ForceCondition*) 170 { 171 G4double attLength = DBL_MAX; 172 G4MaterialPropertiesTable* MPT = 173 aTrack.GetMaterial()->GetMaterialPropertiesTable(); 174 if(MPT) 175 { 176 G4MaterialPropertyVector* attVector = MPT->GetProperty(kMIEHG); 177 if(attVector) 178 { 179 attLength = attVector->Value( 180 aTrack.GetDynamicParticle()->GetTotalEnergy(), idx_mie); 181 } 182 } 183 return attLength; 184 } 185 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 187 void G4OpMieHG::SetVerboseLevel(G4int verbose) 188 { 189 verboseLevel = verbose; 190 G4OpticalParameters::Instance()->SetMieVerboseLevel(verboseLevel); 191 } 192