119 const G4Step* pStep = &aStep;
123 if (hStep) pStep = hStep;
144 GetMaterialPropertiesTable();
146 GetMaterialPropertiesTable();
153 G4cout <<
" vol1: " << volnam1 <<
", vol2: " << volnam2 <<
G4endl;
173 std::vector<G4Navigator*>::iterator iNav =
175 GetActiveNavigatorsIterator();
178 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
181 theGlobalNormal = -theGlobalNormal;
186 ed <<
" G4UCNBoundaryProcess/PostStepDoIt(): "
187 <<
" The Navigator reports that it returned an invalid normal"
189 G4Exception(
"G4UCNBoundaryProcess::PostStepDoIt",
"UCNBoun01",
191 "Invalid Surface Normal - Geometry must return valid surface normal");
198 if (OldMomentum * theGlobalNormal > 0.0) {
199 #ifdef G4OPTICAL_DEBUG
201 ed <<
" G4UCNBoundaryProcess/PostStepDoIt(): "
202 <<
" theGlobalNormal points in a wrong direction. "
204 ed <<
" The momentum of the photon arriving at interface (oldMomentum)"
205 <<
" must exit the volume cross in the step. " <<
G4endl;
206 ed <<
" So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." <<
G4endl;
207 ed <<
" >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal <<
G4endl;
208 ed <<
" Old Momentum (during step) = " << OldMomentum <<
G4endl;
209 ed <<
" Global Normal (Exiting New Vol) = " << theGlobalNormal <<
G4endl;
211 G4Exception(
"G4UCNBoundaryProcess::PostStepDoIt",
"UCNBoun02",
214 "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
216 theGlobalNormal = -theGlobalNormal;
222 G4double theMomentumNormal = theNeutronMomentum*theGlobalNormal;
224 (OldMomentum * theGlobalNormal);
245 G4double FermiPotDiff = FermiPot2 - FermiPot1;
248 G4cout <<
"UCNBoundaryProcess: new FermiPot: " << FermiPot2/
neV
249 <<
"neV, old FermiPot:" << FermiPot1/
neV <<
"neV" <<
G4endl;
278 theta_i = OldMomentum.
angle(-theGlobalNormal);
283 ConditionsValid(Energy, FermiPotDiff, theta_i)) {
303 GetMRIntProbability(theta_i, Energy);
308 GetMRIntTransProbability(theta_i, Energy);
312 G4cout <<
"Energy: " << Energy/
neV <<
"neV"
313 <<
", Enormal: " << Enormal/
neV <<
"neV" <<
G4endl;
315 G4cout <<
"Reflection_prob: " << MRpDiffuse
316 <<
", Transmission_prob: " << MRpDiffuseTrans <<
G4endl;
320 if (!
High(Enormal, FermiPotDiff)) {
328 if (
Loss(pUpScatter, theVelocityNormal, FermiPotDiff)) {
359 NewMomentum =
MRreflect(MRpDiffuse, OldMomentum, theGlobalNormal,
360 Energy, FermiPotDiff);
365 NewMomentum =
Reflect(pDiffuse, OldMomentum, theGlobalNormal);
389 theGlobalNormal, Energy, FermiPotDiff, Enew);
407 << reflectivity <<
G4endl;
413 NewMomentum =
Reflect(pDiffuse, OldMomentum, theGlobalNormal);
425 G4double mass = -std::sqrt(theMomentumNormal*theMomentumNormal -
431 theNeutronMomentum - (theMomentumNormal-
mass)*theGlobalNormal;
443 G4cout <<
"Energy: " << Energy/
neV <<
"neV, Enormal: "
444 << Enormal/
neV <<
"neV, fpdiff " << FermiPotDiff/
neV
445 <<
"neV, Enew " << Enew/
neV <<
"neV" <<
G4endl;
446 G4cout <<
"UCNBoundaryProcess: Transmit and set the new Energy "
476 G4double vRatio = theVelocityNormal/vBound;
478 G4double pLoss = (2*pUpScatter*vRatio)/(std::sqrt(1-(vRatio*vRatio)));
490 pLoss *= std::sqrt(1+2*b*b*vBound*vBound/
491 (hdm*hdm+0.85*hdm*vBound*w+2*vBound*vBound*w*w));
505 G4double r = (std::sqrt(Enormal) - std::sqrt(Enormal - FermiPot)) /
506 (std::sqrt(Enormal) + std::sqrt(Enormal - FermiPot));
515 G4double PdotN = OldMomentum * Normal;
522 if (NewMomentum == OldMomentum ||
G4UniformRand() < pDiffuse) {
561 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
573 G4double PdotN = OldMomentum * Normal;
575 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
597 G4double costheta = OldMomentum*Normal;
599 G4double Enormal = Energy * (costheta*costheta);
603 (1.-pDiffuse-pDiffuseTrans-pLoss);
609 if (decide < pSpecular) {
613 G4double PdotN = OldMomentum * Normal;
614 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
626 if (decide < pSpecular + pDiffuse) {
633 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
636 <<
", " << NewMomentum <<
G4endl;
646 if (decide < pSpecular + pDiffuse + pDiffuseTrans) {
653 MRDiffTrans(Normal, Energy, FermiPot, OldMomentum, pDiffuseTrans);
655 Enew = Energy - FermiPot;
664 if (decide < pSpecular + pDiffuse + pDiffuseTrans + pLoss) {
679 Enew = Energy - FermiPot;
682 G4double produ = OldMomentum * Normal;
684 NewMomentum = palt*OldMomentum-
692 return NewMomentum.
unit();
722 GetMRMaxProbability(theta_i, Energy)/pDiffuse <=
724 GetMRProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/pDiffuse)
731 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
733 GetMRMaxProbability(theta_i, Energy)) > 1) {
734 G4cout <<
"MRMax Wahrscheinlichkeitsueberschreitung!" <<
G4endl;
736 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
738 GetMRMaxProbability(theta_i, Energy)) << G4endl;
740 SetMRMaxProbability(theta_i, Energy,
742 GetMRProbability(theta_i, Energy,
743 FermiPot, theta_o, phi_o));
746 if(++count > 10000) { accepted =
true; }
768 if (momentum * Normal<0) {
771 G4cout <<
"G4UCNBoundaryProcess::MRDiffRefl: !" <<
G4endl;
774 return momentum.
unit();
804 GetMRMaxTransProbability(theta_i, Energy)/pDiffuseTrans <=
806 GetMRTransProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/
814 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
816 GetMRMaxTransProbability(theta_i, Energy)) > 1) {
817 G4cout <<
"MRMaxTrans Wahrscheinlichkeitsueberschreitung!" <<
G4endl;
819 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
821 GetMRMaxTransProbability(theta_i, Energy)) << G4endl;
823 SetMRMaxTransProbability(theta_i, Energy,
825 GetMRTransProbability(theta_i, Energy,
826 FermiPot, theta_o, phi_o));
829 if(++count > 10000) { accepted =
true; }
846 if (momentum*Normal<0) {
849 G4cout <<
"G4UCNBoundaryProcess::MRDiffTrans: !" <<
G4endl;
852 return momentum.
unit();
857 return (Energy - FermiPot);
869 if (momentum*Normal < 0) {
871 G4cout <<
"G4UCNBoundaryProcess::LDiffRefl: !" <<
G4endl;
874 return momentum.
unit();
920 G4cout <<
" *** No G4UCNMaterialPropertiesTable *** " <<
G4endl;
922 G4cout <<
" *** No MicroRoughness Table *** " <<
G4endl;
924 G4cout <<
" *** MicroRoughness Condition not satisfied *** " <<
G4endl;
936 G4cout <<
" *** MR Model Diffuse Reflection *** " <<
G4endl;
940 G4cout <<
" *** MR Model Diffuse Transmission *** " <<
G4endl;
949 G4cout <<
"Sum NoMRCondition: "
951 G4cout <<
"Sum No. E < V Loss: "
953 G4cout <<
"Sum No. E > V Ezero: "
955 G4cout <<
"Sum No. E < V SpinFlip: "
957 G4cout <<
"Sum No. E > V Specular Reflection: "
959 G4cout <<
"Sum No. E < V Specular Reflection: "
961 G4cout <<
"Sum No. E < V Lambertian Reflection: "
963 G4cout <<
"Sum No. E > V MR DiffuseReflection: "
965 G4cout <<
"Sum No. E < V MR DiffuseReflection: "
967 G4cout <<
"Sum No. E > V SnellTransmit: "
969 G4cout <<
"Sum No. E > V MR SnellTransmit: "
971 G4cout <<
"Sum No. E > V DiffuseTransmit: "
983 if (sd)
return sd->
Hit(&aStep);