#include <G4UCNMicroRoughnessHelper.hh>

Public Member Functions
G4double	S2 (G4double, G4double) const

G4double	SS2 (G4double, G4double) const

G4double	Fmu (G4double, G4double, G4double, G4double, G4double, G4double, G4double) const

G4double	FmuS (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const

G4double	IntIplus (G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const

G4double	ProbIplus (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const

G4double	IntIminus (G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const

G4double	ProbIminus (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const

Static Public Member Functions
static G4UCNMicroRoughnessHelper *	GetInstance ()

Protected Member Functions
	G4UCNMicroRoughnessHelper ()

	~G4UCNMicroRoughnessHelper ()

Detailed Description

Definition at line 56 of file G4UCNMicroRoughnessHelper.hh.

Constructor & Destructor Documentation

G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper ( )

protected

Definition at line 53 of file G4UCNMicroRoughnessHelper.cc.

53 {;}

Here is the caller graph for this function:

G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper ( )

protected

Definition at line 57 of file G4UCNMicroRoughnessHelper.cc.

 {
   delete fpInstance;
   fpInstance = nullptr;
 }

Member Function Documentation

G4double G4UCNMicroRoughnessHelper::Fmu	(	G4double	k2,
		G4double	thetai,
		G4double	thetao,
		G4double	phio,
		G4double	b2,
		G4double	w2,
		G4double	AngCut
	)		const

Definition at line 98 of file G4UCNMicroRoughnessHelper.cc.

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

Here is the caller graph for this function:

G4double G4UCNMicroRoughnessHelper::FmuS	(	G4double	k,
		G4double	kS,
		G4double	thetai,
		G4double	thetaSo,
		G4double	phiSo,
		G4double	b2,
		G4double	w2,
		G4double	AngCut,
		G4double	thetarefract
	)		const

Definition at line 128 of file G4UCNMicroRoughnessHelper.cc.

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

Here is the caller graph for this function:

G4UCNMicroRoughnessHelper * G4UCNMicroRoughnessHelper::GetInstance ( void )

static

Definition at line 65 of file G4UCNMicroRoughnessHelper.cc.

 {
   if (fpInstance == nullptr) fpInstance = new G4UCNMicroRoughnessHelper;
   return fpInstance;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4UCNMicroRoughnessHelper::IntIminus	(	G4double	E,
		G4double	fermipot,
		G4double	theta_i,
		G4int	AngNoTheta,
		G4int	AngNoPhi,
		G4double	b2,
		G4double	w2,
		G4double *	max,
		G4double	AngCut
	)		const

Definition at line 283 of file G4UCNMicroRoughnessHelper.cc.

 {
   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.;
       //G4cout << ang_stepphi/degree << ", " << ang_steptheta/degree 
       //       << ", " << AngCut/degree << G4endl;
       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;
 }

Here is the call graph for this function:

G4double G4UCNMicroRoughnessHelper::IntIplus	(	G4double	E,
		G4double	fermipot,
		G4double	theta_i,
		G4int	AngNoTheta,
		G4int	AngNoPhi,
		G4double	b2,
		G4double	w2,
		G4double *	max,
		G4double	AngCut
	)		const

Definition at line 159 of file G4UCNMicroRoughnessHelper.cc.

 {
   *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);
 
           //G4cout << "S2:  " << S2(costheta_i_squared,klk2) << G4endl;
           //G4cout << "S2:  " << S2(costheta_o_squared,klk2) << G4endl;
           //G4cout << "Fmu: " << Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*sin(theta_o) << G4endl;
           // checks if the new probability is larger than the
           // maximum found so far
 
           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.;
 
       //G4cout << ang_stepphi/degree << ", "
       //       << ang_steptheta/degree << ","
       //       << AngCut/degree << G4endl;
 
       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;
 }

Here is the call graph for this function:

G4double G4UCNMicroRoughnessHelper::ProbIminus	(	G4double	E,
		G4double	fermipot,
		G4double	theta_i,
		G4double	thetas_o,
		G4double	phis_o,
		G4double	b,
		G4double	w,
		G4double	AngCut
	)		const

Definition at line 431 of file G4UCNMicroRoughnessHelper.cc.

 {
   //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);
 }

Here is the call graph for this function:

G4double G4UCNMicroRoughnessHelper::ProbIplus	(	G4double	E,
		G4double	fermipot,
		G4double	theta_i,
		G4double	theta_o,
		G4double	phi_o,
		G4double	b,
		G4double	w,
		G4double	AngCut
	)		const

Definition at line 404 of file G4UCNMicroRoughnessHelper.cc.

 {
   //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);
 }

Here is the call graph for this function:

G4double G4UCNMicroRoughnessHelper::S2	(	G4double	costheta2,
		G4double	klk2
	)		const

Definition at line 74 of file G4UCNMicroRoughnessHelper.cc.

 {
   // 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)));
   else
      return std::norm(2*std::sqrt(costheta2)/(std::sqrt(costheta2) + std::sqrt(std::complex<G4double>(costheta2 - klk2))));
 }

Here is the caller graph for this function:

G4double G4UCNMicroRoughnessHelper::SS2	(	G4double	costhetas2,
		G4double	klks2
	)		const

Definition at line 88 of file G4UCNMicroRoughnessHelper.cc.

 {
   // cf. p. 175 of the Steyerl paper
 
   return 4*costhetas2/
                    (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
 }

Here is the caller graph for this function:

The documentation for this class was generated from the following files:

source/geant4.10.03.p02/source/materials/include/G4UCNMicroRoughnessHelper.hh
source/geant4.10.03.p02/source/materials/src/G4UCNMicroRoughnessHelper.cc

Public Member Functions

Static Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation