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 // File G4OpMieHG.hh 29 // File G4OpMieHG.hh 30 // Description: Discrete Process -- Mie Scatte 30 // Description: Discrete Process -- Mie Scattering of Optical Photons 31 // Created: 2010-07-03 31 // Created: 2010-07-03 32 // Author: Xin Qian 32 // Author: Xin Qian 33 // Based on work from Vlasios Vasileiou 33 // Based on work from Vlasios Vasileiou 34 // 34 // 35 // This subroutine will mimic the Mie scatteri << 35 // This subroutine will mimic the Mie scattering based on 36 // Henyey-Greenstein phase function 36 // Henyey-Greenstein phase function 37 // Forward and backward angles are treated sep 37 // Forward and backward angles are treated separately. 38 // 38 // >> 39 // mail: gum@triumf.ca >> 40 // 39 ////////////////////////////////////////////// 41 //////////////////////////////////////////////////////////////////////// 40 42 41 #include "G4OpMieHG.hh" << 42 #include "G4PhysicalConstants.hh" << 43 #include "G4OpticalParameters.hh" << 44 #include "G4OpProcessSubType.hh" 43 #include "G4OpProcessSubType.hh" 45 44 46 //....oooOO0OOooo........oooOO0OOooo........oo << 45 #include "G4OpMieHG.hh" >> 46 47 G4OpMieHG::G4OpMieHG(const G4String& processNa 47 G4OpMieHG::G4OpMieHG(const G4String& processName, G4ProcessType type) 48 : G4VDiscreteProcess(processName, type) << 48 : G4VDiscreteProcess(processName, type) 49 { 49 { 50 Initialise(); << 50 if (verboseLevel>0) { 51 if(verboseLevel > 0) << 51 G4cout << GetProcessName() << " is created " << G4endl; 52 { << 52 } 53 G4cout << GetProcessName() << " is created << 54 } << 55 SetProcessSubType(fOpMieHG); << 56 } << 57 << 58 //....oooOO0OOooo........oooOO0OOooo........oo << 59 G4OpMieHG::~G4OpMieHG() = default; << 60 53 61 //....oooOO0OOooo........oooOO0OOooo........oo << 54 SetProcessSubType(fOpMieHG); 62 void G4OpMieHG::PreparePhysicsTable(const G4Pa << 63 { << 64 Initialise(); << 65 } 55 } 66 56 67 //....oooOO0OOooo........oooOO0OOooo........oo << 57 G4OpMieHG::~G4OpMieHG(){} 68 void G4OpMieHG::Initialise() << 58 69 { << 59 //////////// 70 SetVerboseLevel(G4OpticalParameters::Instanc << 60 // Methods 71 } << 61 //////////// 72 62 73 //....oooOO0OOooo........oooOO0OOooo........oo << 63 // PostStepDoIt 74 G4VParticleChange* G4OpMieHG::PostStepDoIt(con << 64 // ------------- 75 con << 65 // >> 66 G4VParticleChange* >> 67 G4OpMieHG::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) 76 { 68 { 77 aParticleChange.Initialize(aTrack); << 69 aParticleChange.Initialize(aTrack); 78 70 79 const G4DynamicParticle* aParticle = aTrack. << 71 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 80 const G4MaterialPropertiesTable* MPT = << 72 const G4Material* aMaterial = aTrack.GetMaterial(); 81 aTrack.GetMaterial()->GetMaterialPropertie << 73 G4MaterialPropertiesTable* aMaterialPropertyTable = 82 << 74 aMaterial->GetMaterialPropertiesTable(); 83 G4double forwardRatio = MPT->GetConstPropert << 75 84 << 76 G4double forward_g = 85 if(verboseLevel > 1) << 77 aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD"); 86 { << 78 G4double backward_g = 87 G4cout << "OpMie Scattering Photon!" << G4 << 79 aMaterialPropertyTable->GetConstProperty("MIEHG_BACKWARD"); 88 << " Old Momentum Direction: " << a << 80 G4double ForwardRatio = 89 << G4endl << 81 aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD_RATIO"); 90 << " MIE Old Polarization: " << aPa << 82 91 << G4endl; << 83 if (verboseLevel>0) { 92 } << 84 G4cout << "MIE Scattering Photon!" << G4endl; 93 << 85 G4cout << "MIE Old Momentum Direction: " 94 G4double gg; << 86 << aParticle->GetMomentumDirection() << G4endl; 95 G4int direction; << 87 G4cout << "MIE Old Polarization: " 96 if(G4UniformRand() <= forwardRatio) << 88 << aParticle->GetPolarization() << G4endl; 97 { << 89 } 98 gg = MPT->GetConstProperty(kMIEHG_F << 90 99 direction = 1; << 91 G4double g; 100 } << 92 G4int direction; 101 else << 93 if (G4UniformRand()<=ForwardRatio){ 102 { << 94 g = forward_g; 103 gg = MPT->GetConstProperty(kMIEHG_B << 95 direction = 1; 104 direction = -1; << 96 } else { 105 } << 97 g = backward_g; 106 << 98 direction = -1; 107 G4double r = G4UniformRand(); << 99 } 108 << 100 109 // sample the direction << 101 G4double r = G4UniformRand(); 110 G4double theta; << 102 111 if(gg != 0.) << 103 G4double Theta; 112 { << 104 //sample the direction 113 theta = std::acos(2. * r * (1. + gg) * (1. << 105 if (g!=0) { 114 ((1. - gg + 2. * gg * << 106 Theta = std::acos(2*r*(1+g)*(1+g)*(1-g+g*r)/((1-g+2*g*r)*(1-g+2*g*r)) -1); 115 1.); << 107 } else { 116 } << 108 Theta = std::acos(2*r-1.); 117 else << 109 } 118 { << 110 G4double Phi = G4UniformRand()*2*pi; 119 theta = std::acos(2. * r - 1.); << 111 120 } << 112 if (direction==-1) Theta = pi - Theta; //backward scattering 121 G4double phi = G4UniformRand() * twopi; << 113 122 << 114 G4ThreeVector NewMomentumDirection, OldMomentumDirection; 123 if(direction == -1) << 115 G4ThreeVector OldPolarization, NewPolarization; 124 theta = pi - theta; // backward scatterin << 116 125 << 117 NewMomentumDirection.set 126 G4ThreeVector newMomDir, oldMomDir; << 118 (std::sin(Theta)*std::cos(Phi), std::sin(Theta)*std::sin(Phi), std::cos(Theta)); 127 G4ThreeVector newPol, oldPol; << 119 OldMomentumDirection = aParticle->GetMomentumDirection(); 128 << 120 NewMomentumDirection.rotateUz(OldMomentumDirection); 129 G4double sinth = std::sin(theta); << 121 NewMomentumDirection = NewMomentumDirection.unit(); 130 newMomDir.set(sinth * std::cos(phi), sinth * << 122 131 oldMomDir = aParticle->GetMomentumDirection( << 123 OldPolarization = aParticle->GetPolarization(); 132 newMomDir.rotateUz(oldMomDir); << 124 G4double constant = -1./NewMomentumDirection.dot(OldPolarization); 133 newMomDir = newMomDir.unit(); << 125 134 << 126 NewPolarization = NewMomentumDirection + constant*OldPolarization; 135 oldPol = aParticle->GetPolarization(); << 127 NewPolarization = NewPolarization.unit(); 136 newPol = newMomDir - oldPol / newMomDir.dot( << 128 137 newPol = newPol.unit(); << 129 if (NewPolarization.mag()==0) { 138 << 130 r = G4UniformRand()*twopi; 139 if(newPol.mag() == 0.) << 131 NewPolarization.set(std::cos(r),std::sin(r),0.); 140 { << 132 NewPolarization.rotateUz(NewMomentumDirection); 141 r = G4UniformRand() * twopi; << 133 } else { 142 newPol.set(std::cos(r), std::sin(r), 0.); << 134 // There are two directions which perpendicular 143 newPol.rotateUz(newMomDir); << 135 // new momentum direction 144 } << 136 if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization; 145 else << 137 } 146 { << 138 147 // There are two directions perpendicular << 139 aParticleChange.ProposePolarization(NewPolarization); 148 if(G4UniformRand() < 0.5) << 140 aParticleChange.ProposeMomentumDirection(NewMomentumDirection); 149 newPol = -newPol; << 141 150 } << 142 if (verboseLevel>0) { 151 << 143 G4cout << "MIE New Polarization: " 152 aParticleChange.ProposePolarization(newPol); << 144 << NewPolarization << G4endl; 153 aParticleChange.ProposeMomentumDirection(new << 145 G4cout << "MIE Polarization Change: " 154 << 146 << *(aParticleChange.GetPolarization()) << G4endl; 155 if(verboseLevel > 1) << 147 G4cout << "MIE New Momentum Direction: " 156 { << 148 << NewMomentumDirection << G4endl; 157 G4cout << "OpMie New Polarization: " << ne << 149 G4cout << "MIE Momentum Change: " 158 << " Polarization Change: " << *(aP << 150 << *(aParticleChange.GetMomentumDirection()) << G4endl; 159 << G4endl << " New Momentum Directi << 151 } 160 << " Momentum Change: " << *(aParti << 161 << G4endl; << 162 } << 163 152 164 return G4VDiscreteProcess::PostStepDoIt(aTra << 153 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 165 } 154 } 166 155 167 //....oooOO0OOooo........oooOO0OOooo........oo << 156 // GetMeanFreePath() 168 G4double G4OpMieHG::GetMeanFreePath(const G4Tr << 157 // ----------------- 169 G4ForceCon << 158 // 170 { << 159 G4double G4OpMieHG::GetMeanFreePath(const G4Track& aTrack, 171 G4double attLength = DBL_MAX; << 160 G4double , 172 G4MaterialPropertiesTable* MPT = << 161 G4ForceCondition* ) 173 aTrack.GetMaterial()->GetMaterialPropertie << 162 { 174 if(MPT) << 163 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 175 { << 164 const G4Material* aMaterial = aTrack.GetMaterial(); 176 G4MaterialPropertyVector* attVector = MPT- << 165 177 if(attVector) << 166 G4double thePhotonEnergy = aParticle->GetTotalEnergy(); 178 { << 167 179 attLength = attVector->Value( << 168 G4double AttenuationLength = DBL_MAX; 180 aTrack.GetDynamicParticle()->GetTotalE << 169 181 } << 170 G4MaterialPropertiesTable* aMaterialPropertyTable = 182 } << 171 aMaterial->GetMaterialPropertiesTable(); 183 return attLength; << 172 184 } << 173 if (aMaterialPropertyTable) { >> 174 G4MaterialPropertyVector* AttenuationLengthVector = >> 175 aMaterialPropertyTable->GetProperty("MIEHG"); >> 176 if (AttenuationLengthVector) { >> 177 AttenuationLength = AttenuationLengthVector -> >> 178 GetProperty(thePhotonEnergy); >> 179 } else { >> 180 // G4cout << "No Mie scattering length specified" << G4endl; >> 181 } >> 182 } else { >> 183 // G4cout << "No Mie scattering length specified" << G4endl; >> 184 } 185 185 186 //....oooOO0OOooo........oooOO0OOooo........oo << 186 // G4cout << thePhotonEnergy/GeV << " \t" << AttenuationLength/m << G4endl; 187 void G4OpMieHG::SetVerboseLevel(G4int verbose) << 187 188 { << 188 return AttenuationLength; 189 verboseLevel = verbose; << 190 G4OpticalParameters::Instance()->SetMieVerbo << 191 } 189 } 192 190