00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00042
00043 #include "G4OpMieHG.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4OpProcessSubType.hh"
00046
00047 G4OpMieHG::G4OpMieHG(const G4String& processName, G4ProcessType type)
00048 : G4VDiscreteProcess(processName, type)
00049 {
00050 if (verboseLevel>0) {
00051 G4cout << GetProcessName() << " is created " << G4endl;
00052 }
00053
00054 SetProcessSubType(fOpMieHG);
00055 }
00056
00057 G4OpMieHG::~G4OpMieHG(){}
00058
00060
00062
00063
00064
00065
00066 G4VParticleChange*
00067 G4OpMieHG::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00068 {
00069 aParticleChange.Initialize(aTrack);
00070
00071 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00072 const G4Material* aMaterial = aTrack.GetMaterial();
00073 G4MaterialPropertiesTable* aMaterialPropertyTable =
00074 aMaterial->GetMaterialPropertiesTable();
00075
00076 G4double forward_g =
00077 aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD");
00078 G4double backward_g =
00079 aMaterialPropertyTable->GetConstProperty("MIEHG_BACKWARD");
00080 G4double ForwardRatio =
00081 aMaterialPropertyTable->GetConstProperty("MIEHG_FORWARD_RATIO");
00082
00083 if (verboseLevel>0) {
00084 G4cout << "MIE Scattering Photon!" << G4endl;
00085 G4cout << "MIE Old Momentum Direction: "
00086 << aParticle->GetMomentumDirection() << G4endl;
00087 G4cout << "MIE Old Polarization: "
00088 << aParticle->GetPolarization() << G4endl;
00089 }
00090
00091 G4double gg;
00092 G4int direction;
00093 if (G4UniformRand()<=ForwardRatio){
00094 gg = forward_g;
00095 direction = 1;
00096 } else {
00097 gg = backward_g;
00098 direction = -1;
00099 }
00100
00101 G4double r = G4UniformRand();
00102
00103 G4double Theta;
00104
00105 if (gg!=0) {
00106 Theta = std::acos(2*r*(1+gg)*(1+gg)*(1-gg+gg*r)/((1-gg+2*gg*r)*(1-gg+2*gg*r)) -1);
00107 } else {
00108 Theta = std::acos(2*r-1.);
00109 }
00110 G4double Phi = G4UniformRand()*2*pi;
00111
00112 if (direction==-1) Theta = pi - Theta;
00113
00114 G4ThreeVector NewMomentumDirection, OldMomentumDirection;
00115 G4ThreeVector OldPolarization, NewPolarization;
00116
00117 NewMomentumDirection.set
00118 (std::sin(Theta)*std::cos(Phi), std::sin(Theta)*std::sin(Phi), std::cos(Theta));
00119 OldMomentumDirection = aParticle->GetMomentumDirection();
00120 NewMomentumDirection.rotateUz(OldMomentumDirection);
00121 NewMomentumDirection = NewMomentumDirection.unit();
00122
00123 OldPolarization = aParticle->GetPolarization();
00124 G4double constant = -1./NewMomentumDirection.dot(OldPolarization);
00125
00126 NewPolarization = NewMomentumDirection + constant*OldPolarization;
00127 NewPolarization = NewPolarization.unit();
00128
00129 if (NewPolarization.mag()==0) {
00130 r = G4UniformRand()*twopi;
00131 NewPolarization.set(std::cos(r),std::sin(r),0.);
00132 NewPolarization.rotateUz(NewMomentumDirection);
00133 } else {
00134
00135
00136 if (G4UniformRand() < 0.5) NewPolarization = -NewPolarization;
00137 }
00138
00139 aParticleChange.ProposePolarization(NewPolarization);
00140 aParticleChange.ProposeMomentumDirection(NewMomentumDirection);
00141
00142 if (verboseLevel>0) {
00143 G4cout << "MIE New Polarization: "
00144 << NewPolarization << G4endl;
00145 G4cout << "MIE Polarization Change: "
00146 << *(aParticleChange.GetPolarization()) << G4endl;
00147 G4cout << "MIE New Momentum Direction: "
00148 << NewMomentumDirection << G4endl;
00149 G4cout << "MIE Momentum Change: "
00150 << *(aParticleChange.GetMomentumDirection()) << G4endl;
00151 }
00152
00153 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00154 }
00155
00156
00157
00158
00159 G4double G4OpMieHG::GetMeanFreePath(const G4Track& aTrack,
00160 G4double ,
00161 G4ForceCondition* )
00162 {
00163 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00164 const G4Material* aMaterial = aTrack.GetMaterial();
00165
00166 G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00167
00168 G4double AttenuationLength = DBL_MAX;
00169
00170 G4MaterialPropertiesTable* aMaterialPropertyTable =
00171 aMaterial->GetMaterialPropertiesTable();
00172
00173 if (aMaterialPropertyTable) {
00174 G4MaterialPropertyVector* AttenuationLengthVector =
00175 aMaterialPropertyTable->GetProperty("MIEHG");
00176 if (AttenuationLengthVector) {
00177 AttenuationLength = AttenuationLengthVector ->
00178 Value(thePhotonEnergy);
00179 } else {
00180
00181 }
00182 } else {
00183
00184 }
00185
00186
00187
00188 return AttenuationLength;
00189 }