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);

で最大値を $1$ か $4\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になっている）を元の法線ベクトルに加えて、