| 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 ////////////////////////////////////////////// 39 ////////////////////////////////////////////////////////////////////////
40 40
41 #include "G4OpMieHG.hh" 41 #include "G4OpMieHG.hh"
42 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh"
43 #include "G4OpticalParameters.hh" 43 #include "G4OpticalParameters.hh"
44 #include "G4OpProcessSubType.hh" 44 #include "G4OpProcessSubType.hh"
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 Initialise();
51 if(verboseLevel > 0) 51 if(verboseLevel > 0)
52 { 52 {
53 G4cout << GetProcessName() << " is created 53 G4cout << GetProcessName() << " is created " << G4endl;
54 } 54 }
55 SetProcessSubType(fOpMieHG); 55 SetProcessSubType(fOpMieHG);
56 } 56 }
57 57
58 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 G4OpMieHG::~G4OpMieHG() = default; << 59 G4OpMieHG::~G4OpMieHG() {}
60 60
61 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 void G4OpMieHG::PreparePhysicsTable(const G4Pa 62 void G4OpMieHG::PreparePhysicsTable(const G4ParticleDefinition&)
63 { 63 {
64 Initialise(); 64 Initialise();
65 } 65 }
66 66
67 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 void G4OpMieHG::Initialise() 68 void G4OpMieHG::Initialise()
69 { 69 {
70 SetVerboseLevel(G4OpticalParameters::Instanc 70 SetVerboseLevel(G4OpticalParameters::Instance()->GetMieVerboseLevel());
71 } 71 }
72 72
73 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 G4VParticleChange* G4OpMieHG::PostStepDoIt(con 74 G4VParticleChange* G4OpMieHG::PostStepDoIt(const G4Track& aTrack,
75 con 75 const G4Step& aStep)
76 { 76 {
77 aParticleChange.Initialize(aTrack); 77 aParticleChange.Initialize(aTrack);
78 78
79 const G4DynamicParticle* aParticle = aTrack. 79 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
80 const G4MaterialPropertiesTable* MPT = 80 const G4MaterialPropertiesTable* MPT =
81 aTrack.GetMaterial()->GetMaterialPropertie 81 aTrack.GetMaterial()->GetMaterialPropertiesTable();
82 82
83 G4double forwardRatio = MPT->GetConstPropert 83 G4double forwardRatio = MPT->GetConstProperty(kMIEHG_FORWARD_RATIO);
84 84
85 if(verboseLevel > 1) 85 if(verboseLevel > 1)
86 { 86 {
87 G4cout << "OpMie Scattering Photon!" << G4 87 G4cout << "OpMie Scattering Photon!" << G4endl
88 << " Old Momentum Direction: " << a 88 << " Old Momentum Direction: " << aParticle->GetMomentumDirection()
89 << G4endl 89 << G4endl
90 << " MIE Old Polarization: " << aPa 90 << " MIE Old Polarization: " << aParticle->GetPolarization()
91 << G4endl; 91 << G4endl;
92 } 92 }
93 93
94 G4double gg; 94 G4double gg;
95 G4int direction; 95 G4int direction;
96 if(G4UniformRand() <= forwardRatio) 96 if(G4UniformRand() <= forwardRatio)
97 { 97 {
98 gg = MPT->GetConstProperty(kMIEHG_F 98 gg = MPT->GetConstProperty(kMIEHG_FORWARD);
99 direction = 1; 99 direction = 1;
100 } 100 }
101 else 101 else
102 { 102 {
103 gg = MPT->GetConstProperty(kMIEHG_B 103 gg = MPT->GetConstProperty(kMIEHG_BACKWARD);
104 direction = -1; 104 direction = -1;
105 } 105 }
106 106
107 G4double r = G4UniformRand(); 107 G4double r = G4UniformRand();
108 108
109 // sample the direction 109 // sample the direction
110 G4double theta; 110 G4double theta;
111 if(gg != 0.) 111 if(gg != 0.)
112 { 112 {
113 theta = std::acos(2. * r * (1. + gg) * (1. 113 theta = std::acos(2. * r * (1. + gg) * (1. + gg) * (1. - gg + gg * r) /
114 ((1. - gg + 2. * gg * 114 ((1. - gg + 2. * gg * r) * (1. - gg + 2. * gg * r)) -
115 1.); 115 1.);
116 } 116 }
117 else 117 else
118 { 118 {
119 theta = std::acos(2. * r - 1.); 119 theta = std::acos(2. * r - 1.);
120 } 120 }
121 G4double phi = G4UniformRand() * twopi; 121 G4double phi = G4UniformRand() * twopi;
122 122
123 if(direction == -1) 123 if(direction == -1)
124 theta = pi - theta; // backward scatterin 124 theta = pi - theta; // backward scattering
125 125
126 G4ThreeVector newMomDir, oldMomDir; 126 G4ThreeVector newMomDir, oldMomDir;
127 G4ThreeVector newPol, oldPol; 127 G4ThreeVector newPol, oldPol;
128 128
129 G4double sinth = std::sin(theta); 129 G4double sinth = std::sin(theta);
130 newMomDir.set(sinth * std::cos(phi), sinth * 130 newMomDir.set(sinth * std::cos(phi), sinth * std::sin(phi), std::cos(theta));
131 oldMomDir = aParticle->GetMomentumDirection( 131 oldMomDir = aParticle->GetMomentumDirection();
132 newMomDir.rotateUz(oldMomDir); 132 newMomDir.rotateUz(oldMomDir);
133 newMomDir = newMomDir.unit(); 133 newMomDir = newMomDir.unit();
134 134
135 oldPol = aParticle->GetPolarization(); 135 oldPol = aParticle->GetPolarization();
136 newPol = newMomDir - oldPol / newMomDir.dot( 136 newPol = newMomDir - oldPol / newMomDir.dot(oldPol);
137 newPol = newPol.unit(); 137 newPol = newPol.unit();
138 138
139 if(newPol.mag() == 0.) 139 if(newPol.mag() == 0.)
140 { 140 {
141 r = G4UniformRand() * twopi; 141 r = G4UniformRand() * twopi;
142 newPol.set(std::cos(r), std::sin(r), 0.); 142 newPol.set(std::cos(r), std::sin(r), 0.);
143 newPol.rotateUz(newMomDir); 143 newPol.rotateUz(newMomDir);
144 } 144 }
145 else 145 else
146 { 146 {
147 // There are two directions perpendicular 147 // There are two directions perpendicular to new momentum direction
148 if(G4UniformRand() < 0.5) 148 if(G4UniformRand() < 0.5)
149 newPol = -newPol; 149 newPol = -newPol;
150 } 150 }
151 151
152 aParticleChange.ProposePolarization(newPol); 152 aParticleChange.ProposePolarization(newPol);
153 aParticleChange.ProposeMomentumDirection(new 153 aParticleChange.ProposeMomentumDirection(newMomDir);
154 154
155 if(verboseLevel > 1) 155 if(verboseLevel > 1)
156 { 156 {
157 G4cout << "OpMie New Polarization: " << ne 157 G4cout << "OpMie New Polarization: " << newPol << G4endl
158 << " Polarization Change: " << *(aP 158 << " Polarization Change: " << *(aParticleChange.GetPolarization())
159 << G4endl << " New Momentum Directi 159 << G4endl << " New Momentum Direction: " << newMomDir << G4endl
160 << " Momentum Change: " << *(aParti 160 << " Momentum Change: " << *(aParticleChange.GetMomentumDirection())
161 << G4endl; 161 << G4endl;
162 } 162 }
163 163
164 return G4VDiscreteProcess::PostStepDoIt(aTra 164 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
165 } 165 }
166 166
167 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
168 G4double G4OpMieHG::GetMeanFreePath(const G4Tr 168 G4double G4OpMieHG::GetMeanFreePath(const G4Track& aTrack, G4double,
169 G4ForceCon 169 G4ForceCondition*)
170 { 170 {
171 G4double attLength = DBL_MAX; 171 G4double attLength = DBL_MAX;
172 G4MaterialPropertiesTable* MPT = 172 G4MaterialPropertiesTable* MPT =
173 aTrack.GetMaterial()->GetMaterialPropertie 173 aTrack.GetMaterial()->GetMaterialPropertiesTable();
174 if(MPT) 174 if(MPT)
175 { 175 {
176 G4MaterialPropertyVector* attVector = MPT- 176 G4MaterialPropertyVector* attVector = MPT->GetProperty(kMIEHG);
177 if(attVector) 177 if(attVector)
178 { 178 {
179 attLength = attVector->Value( 179 attLength = attVector->Value(
180 aTrack.GetDynamicParticle()->GetTotalE 180 aTrack.GetDynamicParticle()->GetTotalEnergy(), idx_mie);
181 } 181 }
182 } 182 }
183 return attLength; 183 return attLength;
184 } 184 }
185 185
186 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
187 void G4OpMieHG::SetVerboseLevel(G4int verbose) 187 void G4OpMieHG::SetVerboseLevel(G4int verbose)
188 { 188 {
189 verboseLevel = verbose; 189 verboseLevel = verbose;
190 G4OpticalParameters::Instance()->SetMieVerbo 190 G4OpticalParameters::Instance()->SetMieVerboseLevel(verboseLevel);
191 } 191 }
192 192