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