| 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" 43 #include "G4OpMieHG.hh"
42 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
43 #include "G4OpticalParameters.hh" <<
44 #include "G4OpProcessSubType.hh" 45 #include "G4OpProcessSubType.hh"
45 46
46 //....oooOO0OOooo........oooOO0OOooo........oo <<
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 53
58 //....oooOO0OOooo........oooOO0OOooo........oo << 54 SetProcessSubType(fOpMieHG);
59 G4OpMieHG::~G4OpMieHG() = default; <<
60 <<
61 //....oooOO0OOooo........oooOO0OOooo........oo <<
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 gg;
100 } << 92 G4int direction;
101 else << 93 if (G4UniformRand()<=ForwardRatio){
102 { << 94 gg = forward_g;
103 gg = MPT->GetConstProperty(kMIEHG_B << 95 direction = 1;
104 direction = -1; << 96 } else {
105 } << 97 gg = 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 (gg!=0) {
114 ((1. - gg + 2. * gg * << 106 Theta = std::acos(2*r*(1+gg)*(1+gg)*(1-gg+gg*r)/((1-gg+2*gg*r)*(1-gg+2*gg*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 Value(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