Geant4での可視光の散乱(1)
Geant4のoptical photonでは、G4OpticalSurfaceの設定によって散乱の仕方が変わります。"unified"の場合は、境界面を"polished"、"polishedackpainted"、"polishedfrontpainted"にしている場合、というparameter(単位はrad)を導入することで、その散乱「具合」を変更可能です。また"glisur"の場合は、"polish"の数値(0〜1)で変わります。
また"unified"では、"ground"、"groundbackpainted"、"groundfrontpainted"の場合に"SPECULARLOBECONSTANT"、"SPECULARSPIKECONSTANT"、"BACKSCATTERCONSTANT"を全く入れないと、問答無用でランバート反射*1が使われます。
実際に計算しているのはG4OpBoundaryProcess.ccの中の、G4OpBoundaryProcess::GetFacetNormalです(以下は見やすくするために若干修正)。
G4ThreeVector G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum, const G4ThreeVector& Normal ) const { G4ThreeVector FacetNormal; if (theModel == unified) { G4double alpha; G4double sigma_alpha = OpticalSurface ? OpticalSurface->GetSigmaAlpha() : 0.; G4double f_max = std::min(1.0,4.*sigma_alpha); do { do { alpha = G4RandGauss::shoot(0.0,sigma_alpha); } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi ); G4double phi = G4UniformRand()*twopi; G4double SinAlpha = std::sin(alpha); G4double CosAlpha = std::cos(alpha); G4double SinPhi = std::sin(phi); G4double CosPhi = std::cos(phi); G4double unit_x = SinAlpha * CosPhi; G4double unit_y = SinAlpha * SinPhi; G4double unit_z = CosAlpha; FacetNormal.setX(unit_x); FacetNormal.setY(unit_y); FacetNormal.setZ(unit_z); G4ThreeVector tmpNormal = Normal; FacetNormal.rotateUz(tmpNormal); } while (Momentum * FacetNormal >= 0.0); } else { G4double polish = OpticalSurface ? OpticalSurface->GetPolish() : 1.; if (polish < 1.0) { do { G4ThreeVector smear; do { smear.setX(2.*G4UniformRand()-1.0); smear.setY(2.*G4UniformRand()-1.0); smear.setZ(2.*G4UniformRand()-1.0); } while (smear.mag()>1.0); smear = (1.-polish) * smear; FacetNormal = Normal + smear; } while (Momentum * FacetNormal >= 0.0); FacetNormal = FacetNormal.unit(); } else { FacetNormal = Normal; } // if } // if return FacetNormal; }
"unified" + "polished"の場合
もともとの境界面が持つ法線ベクトルに対して、SetSigmaAlphaから設定できる、の値からそのずれを計算します。
まず
G4double sigma_alpha = 0.0; if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
でを取出し、
G4double f_max = std::min(1.0,4.*sigma_alpha);
で最大値をかのどちらか小さいほうに設定します。
次に
do { alpha = G4RandGauss::shoot(0.0,sigma_alpha); } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
散乱角を標準偏差のGaussianで発生させ、
からの間の乱数
の2条件を満たすときのみ、を採用します。
結局これは、のときの確率分布が
に従うような分布になります(最初のコメント文に書いてありますが)。
完全鏡面だった場合の法線ベクトルに対して、得られた散乱角のベクトル分だけ回転させて、最終的な散乱面の法線を決定しています。
"glisur" + "polish"の場合
"unified"と同様に、SetPolishで設定できるpolishの値から法線を決定します。
まずはここで、
G4ThreeVector smear; do { smear.setX(2.*G4UniformRand()-1.0); smear.setY(2.*G4UniformRand()-1.0); smear.setZ(2.*G4UniformRand()-1.0); } while (smear.mag()>1.0); smear = (1.-polish) * smear;
一様等方なベクトルを乱数で生成します。この得られたベクトル(長さは1 - polishになっている)を元の法線ベクトルに加えて、
FacetNormal = Normal + smear;
後は数学的におかしくならないように微調整をして終わりです。polish = 0の場合が完全な散乱のように一瞬勘違いしますが、FacetNormalにはNormalの成分が半分くらい残っていますので、polish = -∞にしない限り完全な等方散乱にはならないように見えます。
*1:反射面をどの角度に傾けても(つまり1/cosθだけ反射光束の密度が増えても)、反射角の分布がcosθに比例するため、どの角度から見えても同じ反射光量に見える理想的な反射。