G4UCNMicroRoughnessHelper.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
 
// This file contains the source code of various functions all related to the
// calculation of microroughness.
//
// see A. Steyerl, Z. Physik 254 (1972) 169.
//
// A description of the functions can be found in the corresponding header file
 
#include "G4UCNMicroRoughnessHelper.hh"
 
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "globals.hh"
 
G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper() { fpInstance = nullptr; }
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::GetInstance()
{
  if (fpInstance == nullptr) {
    fpInstance = new G4UCNMicroRoughnessHelper;
  }
  return fpInstance;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::S2(G4double costheta2, G4double klk2) const
{
  // case 1: radicand is positive,
  // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
 
  if (costheta2 >= klk2) {
    return 4 * costheta2 / (2 * costheta2 - klk2 + 2 * std::sqrt(costheta2 * (costheta2 - klk2)));
  }
 
  return std::norm(2 * std::sqrt(costheta2) /
                   (std::sqrt(costheta2) + std::sqrt(std::complex<G4double>(costheta2 - klk2))));
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::SS2(G4double costhetas2, G4double klks2) const
{
  // cf. p. 175 of the Steyerl paper
 
  return 4 * costhetas2 /
         (2 * costhetas2 + klks2 + 2 * std::sqrt(costhetas2 * (costhetas2 + klks2)));
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::Fmu(G4double k2, G4double thetai, G4double thetao,
  G4double phio, G4double b2, G4double w2, G4double AngCut) const
{
  G4double mu_squared;
 
  // Checks if the distribution is peaked around the specular direction
 
  if ((std::fabs(thetai - thetao) < AngCut) && (std::fabs(phio) < AngCut)) {
    mu_squared = 0.;
  }
  else {
    // cf. p. 177 of the Steyerl paper
 
    G4double sinthetai = std::sin(thetai);
    G4double sinthetao = std::sin(thetao);
    mu_squared = k2 * (sinthetai * sinthetai + sinthetao * sinthetao -
                        2. * sinthetai * sinthetao * std::cos(phio));
  }
 
  // cf. p. 177 of the Steyerl paper
 
  return b2 * w2 / twopi * std::exp(-mu_squared * w2 / 2);
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::FmuS(G4double k, G4double kS, G4double thetai, G4double thetaSo,
  G4double phiSo, G4double b2, G4double w2, G4double AngCut, G4double thetarefract) const
{
  G4double mu_squared;
 
  // Checks if the distribution is peaked around the direction of
  // unperturbed refraction
 
  if ((std::fabs(thetarefract - thetaSo) < AngCut) && (std::fabs(phiSo) < AngCut)) {
    mu_squared = 0.;
  }
  else {
    G4double sinthetai = std::sin(thetai);
    G4double sinthetaSo = std::sin(thetaSo);
 
    // cf. p. 177 of the Steyerl paper
    mu_squared = k * k * sinthetai * sinthetai + kS * kS * sinthetaSo * sinthetaSo -
                 2. * k * kS * sinthetai * sinthetaSo * std::cos(phiSo);
  }
 
  // cf. p. 177 of the Steyerl paper
 
  return b2 * w2 / twopi * std::exp(-mu_squared * w2 / 2);
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::IntIplus(G4double E, G4double fermipot, G4double theta_i,
  G4int AngNoTheta, G4int AngNoPhi, G4double b2, G4double w2, G4double* max, G4double AngCut) const
{
  *max = 0.;
 
  // max_theta_o saves the theta-position of the max probability,
  // the previous value is saved in a_max_theta_o
 
  G4double a_max_theta_o, max_theta_o = theta_i, a_max_phi_o, max_phi_o = 0.;
 
  // max_phi_o saves the phi-position of the max probability,
  // the previous value is saved in a_max_phi_o
 
  // Definition of the stepsizes in theta_o and phi_o
 
  G4double theta_o;
  G4double phi_o;
  G4double Intens;
  G4double ang_steptheta = 90. * degree / (AngNoTheta - 1);
  G4double ang_stepphi = 360. * degree / (AngNoPhi - 1);
  G4double costheta_i = std::cos(theta_i);
  G4double costheta_i_squared = costheta_i * costheta_i;
 
  // (k_l/k)^2
  G4double kl4d4 =
    neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
 
  // (k_l/k)^2
  G4double klk2 = fermipot / E;
 
  G4double costheta_o_squared;
 
  // k^2
  G4double k2 = 2 * neutron_mass_c2 * E / hbarc_squared;
 
  G4double wkeit = 0.;
 
  // Loop through theta_o
 
  for (theta_o = 0. * degree; theta_o <= 90. * degree + 1e-6; theta_o += ang_steptheta) {
    costheta_o_squared = std::cos(theta_o) * std::cos(theta_o);
 
    // Loop through phi_o
 
    for (phi_o = -180. * degree; phi_o <= 180. * degree + 1e-6; phi_o += ang_stepphi) {
      // calculates probability for a certain theta_o,phi_o pair
 
      Intens = kl4d4 / costheta_i * S2(costheta_i_squared, klk2) * S2(costheta_o_squared, klk2) *
               Fmu(k2, theta_i, theta_o, phi_o, b2, w2, AngCut) * std::sin(theta_o);
 
      if (Intens > *max) {
        *max = Intens;
        max_theta_o = theta_o;
        max_phi_o = phi_o;
      }
 
      // Adds the newly found probability to the integral probability
 
      wkeit += Intens * ang_steptheta * ang_stepphi;
    }
  }
 
  // Fine-Iteration to find maximum of distribution
  // only if the energy is not zero
 
  if (E > 1e-10 * eV) {
    // Break-condition for refining
 
    // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    while ((ang_stepphi >= AngCut * AngCut) || (ang_steptheta >= AngCut * AngCut)) {
      a_max_theta_o = max_theta_o;
      a_max_phi_o = max_phi_o;
      ang_stepphi /= 2.;
      ang_steptheta /= 2.;
 
      for (theta_o = a_max_theta_o - ang_steptheta; theta_o <= a_max_theta_o - ang_steptheta + 1e-6;
           theta_o += ang_steptheta)
      {
        // G4cout << "theta_o: " << theta_o/degree << G4endl;
        costheta_o_squared = std::cos(theta_o) * std::cos(theta_o);
        for (phi_o = a_max_phi_o - ang_stepphi; phi_o <= a_max_phi_o + ang_stepphi + 1e-6;
             phi_o += ang_stepphi)
        {
          // G4cout << "phi_o: " << phi_o/degree << G4endl;
          Intens = kl4d4 / costheta_i * S2(costheta_i_squared, klk2) *
                   S2(costheta_o_squared, klk2) * Fmu(k2, theta_i, theta_o, phi_o, b2, w2, AngCut) *
                   std::sin(theta_o);
          if (Intens > *max) {
            *max = Intens;
            max_theta_o = theta_o;
            max_phi_o = phi_o;
          }
        }
      }
    }
  }
  return wkeit;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::IntIminus(G4double E, G4double fermipot, G4double theta_i,
  G4int AngNoTheta, G4int AngNoPhi, G4double b2, G4double w2, G4double* max, G4double AngCut) const
{
  G4double a_max_thetas_o, max_thetas_o = theta_i;
  G4double a_max_phis_o, max_phis_o = 0.;
  G4double thetas_o;
  G4double phis_o;
  G4double IntensS;
  G4double ang_steptheta = 180. * degree / (AngNoTheta - 1);
  G4double ang_stepphi = 180. * degree / (AngNoPhi - 1);
  G4double costheta_i = std::cos(theta_i);
  G4double costheta_i_squared = costheta_i * costheta_i;
 
  *max = 0.;
  G4double wkeit = 0.;
 
  if (E < fermipot) {
    return wkeit;
  }
 
  // k_l^4/4
  G4double kl4d4 =
    neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
  // (k_l/k)^2
  G4double klk2 = fermipot / E;
 
  // (k_l/k')^2
  G4double klks2 = fermipot / (E - fermipot);
 
  // k'/k
  G4double ksdk = std::sqrt((E - fermipot) / E);
 
  G4double costhetas_o_squared;
 
  // k
  G4double k = std::sqrt(2 * neutron_mass_c2 * E / hbarc_squared);
 
  // k'
  G4double kS = ksdk * k;
 
  for (thetas_o = 0. * degree; thetas_o <= 90. * degree + 1e-6; thetas_o += ang_steptheta) {
    costhetas_o_squared = std::cos(thetas_o) * std::cos(thetas_o);
 
    for (phis_o = -180. * degree; phis_o <= 180. * degree + 1e-6; phis_o += ang_stepphi) {
      // cf. Steyerl-paper p. 176, eq. 21
      if (costhetas_o_squared >= -klks2) {
        // calculates probability for a certain theta'_o, phi'_o pair
 
        G4double thetarefract = thetas_o;
        if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
          thetarefract = std::asin(std::sin(theta_i) / ksdk);
        }
 
        IntensS = kl4d4 / costheta_i * ksdk * S2(costheta_i_squared, klk2) *
                  SS2(costhetas_o_squared, klks2) *
                  FmuS(k, kS, theta_i, thetas_o, phis_o, b2, w2, AngCut, thetarefract) *
                  std::sin(thetas_o);
      }
      else {
        IntensS = 0.;
      }
      // checks if the new probability is larger than
      // the maximum found so far
      if (IntensS > *max) {
        *max = IntensS;
      }
      wkeit += IntensS * ang_steptheta * ang_stepphi;
    }
  }
 
  // Fine-Iteration to find maximum of distribution
 
  if (E > 1e-10 * eV) {
    // Break-condition for refining
 
    while (ang_stepphi >= AngCut * AngCut || ang_steptheta >= AngCut * AngCut) {
      a_max_thetas_o = max_thetas_o;
      a_max_phis_o = max_phis_o;
      ang_stepphi /= 2.;
      ang_steptheta /= 2.;
      
      for (thetas_o = a_max_thetas_o - ang_steptheta;
           thetas_o <= a_max_thetas_o - ang_steptheta + 1e-6; thetas_o += ang_steptheta)
      {
        costhetas_o_squared = std::cos(thetas_o) * std::cos(thetas_o);
        for (phis_o = a_max_phis_o - ang_stepphi; phis_o <= a_max_phis_o + ang_stepphi + 1e-6;
             phis_o += ang_stepphi)
        {
          G4double thetarefract = thetas_o;
          if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
            thetarefract = std::asin(std::sin(theta_i) / ksdk);
          }
 
          IntensS = kl4d4 / costheta_i * ksdk * S2(costheta_i_squared, klk2) *
                    SS2(costhetas_o_squared, klks2) *
                    FmuS(k, kS, theta_i, thetas_o, phis_o, b2, w2, AngCut, thetarefract) *
                    std::sin(thetas_o);
          if (IntensS > *max) {
            *max = IntensS;
            max_thetas_o = thetas_o;
            max_phis_o = phis_o;
          }
        }
      }
    }
  }
  return wkeit;
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::ProbIplus(G4double E, G4double fermipot, G4double theta_i,
  G4double theta_o, G4double phi_o, G4double b, G4double w, G4double AngCut) const
{
  // k_l^4/4
  G4double kl4d4 =
    neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
 
  // (k_l/k)^2
  G4double klk2 = fermipot / E;
 
  G4double costheta_i = std::cos(theta_i);
  G4double costheta_o = std::cos(theta_o);
 
  // eq. 20 on page 176 in the steyerl paper
 
  return kl4d4 / costheta_i * S2(costheta_i * costheta_i, klk2) *
         S2(costheta_o * costheta_o, klk2) *
         Fmu(
           2 * neutron_mass_c2 * E / hbarc_squared, theta_i, theta_o, phi_o, b * b, w * w, AngCut) *
         std::sin(theta_o);
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4double G4UCNMicroRoughnessHelper::ProbIminus(G4double E, G4double fermipot, G4double theta_i,
  G4double thetas_o, G4double phis_o, G4double b, G4double w, G4double AngCut) const
{
  // k_l^4/4
  G4double kl4d4 =
    neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
  // (k_l/k)^2
  G4double klk2 = fermipot / E;
 
  // (k_l/k')^2
  G4double klks2 = fermipot / (E - fermipot);
 
  if (E < fermipot) {
    G4cout << " ProbIminus E < fermipot " << G4endl;
    return 0.;
  }
 
  // k'/k
  G4double ksdk = std::sqrt((E - fermipot) / E);
 
  // k
  G4double k = std::sqrt(2 * neutron_mass_c2 * E / hbarc_squared);
 
  // k'/k
  G4double kS = ksdk * k;
 
  G4double costheta_i = std::cos(theta_i);
  G4double costhetas_o = std::cos(thetas_o);
 
  // eq. 20 on page 176 in the steyerl paper
 
  G4double thetarefract = thetas_o;
  if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
    thetarefract = std::asin(std::sin(theta_i) / ksdk);
  }
 
  return kl4d4 / costheta_i * ksdk * S2(costheta_i * costheta_i, klk2) *
         SS2(costhetas_o * costhetas_o, klks2) *
         FmuS(k, kS, theta_i, thetas_o, phis_o, b * b, w * w, AngCut, thetarefract) *
         std::sin(thetas_o);
}