読者です 読者をやめる 読者になる 読者になる

Geant4での可視光の散乱(1)

Geant4のoptical photonでは、G4OpticalSurfaceの設定によって散乱の仕方が変わります。"unified"の場合は、境界面を"polished"、"polishedackpainted"、"polishedfrontpainted"にしている場合、\sigma_\alphaという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から設定できる、\sigma_\alphaの値からそのずれを計算します。

まず

    G4double sigma_alpha = 0.0;
    if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();

\sigma_\alphaを取出し、

    G4double f_max = std::min(1.0,4.*sigma_alpha);

で最大値を14\sigma_\alphaのどちらか小さいほうに設定します。

次に

      do {
        alpha = G4RandGauss::shoot(0.0,sigma_alpha);
      } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );

散乱角\alpha標準偏差\sigma_\alphaのGaussianで発生させ、
\alpha \lt \frac{\pi}{2}
\sin(\alpha) \lt 0からf_{max}の間の乱数
の2条件を満たすときのみ、\alphaを採用します。

結局これは、\alpha\lt0.25のとき\alphaの確率分布が
\exp(-\frac{\alpha^2}{2\sigma_\alpha^2})\sin(\alpha)
に従うような分布になります(最初のコメント文に書いてありますが)。

完全鏡面だった場合の法線ベクトルに対して、得られた散乱角のベクトル分だけ回転させて、最終的な散乱面の法線を決定しています。

"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θに比例するため、どの角度から見えても同じ反射光量に見える理想的な反射。