120 const G4Step* pStep = &aStep;
124 if (hStep) pStep = hStep;
145 GetMaterialPropertiesTable();
147 GetMaterialPropertiesTable();
154 G4cout <<
" vol1: " << volnam1 <<
", vol2: " << volnam2 <<
G4endl;
174 std::vector<G4Navigator*>::iterator iNav =
176 GetActiveNavigatorsIterator();
179 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
182 theGlobalNormal = -theGlobalNormal;
187 ed <<
" G4UCNBoundaryProcess/PostStepDoIt(): "
188 <<
" The Navigator reports that it returned an invalid normal"
190 G4Exception(
"G4UCNBoundaryProcess::PostStepDoIt",
"UCNBoun01",
192 "Invalid Surface Normal - Geometry must return valid surface normal");
199 if (OldMomentum * theGlobalNormal > 0.0) {
200 #ifdef G4OPTICAL_DEBUG
202 ed <<
" G4UCNBoundaryProcess/PostStepDoIt(): "
203 <<
" theGlobalNormal points in a wrong direction. "
205 ed <<
" The momentum of the photon arriving at interface (oldMomentum)"
206 <<
" must exit the volume cross in the step. " <<
G4endl;
207 ed <<
" So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." <<
G4endl;
208 ed <<
" >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal <<
G4endl;
209 ed <<
" Old Momentum (during step) = " << OldMomentum <<
G4endl;
210 ed <<
" Global Normal (Exiting New Vol) = " << theGlobalNormal <<
G4endl;
212 G4Exception(
"G4UCNBoundaryProcess::PostStepDoIt",
"UCNBoun02",
215 "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
217 theGlobalNormal = -theGlobalNormal;
223 G4double theMomentumNormal = theNeutronMomentum*theGlobalNormal;
225 (OldMomentum * theGlobalNormal);
227 G4double Enormal = theMomentumNormal*theMomentumNormal/2./neutron_mass_c2;
246 G4double FermiPotDiff = FermiPot2 - FermiPot1;
249 G4cout <<
"UCNBoundaryProcess: new FermiPot: " << FermiPot2/
neV
250 <<
"neV, old FermiPot:" << FermiPot1/
neV <<
"neV" <<
G4endl;
279 theta_i = OldMomentum.angle(-theGlobalNormal);
284 ConditionsValid(Energy, FermiPotDiff, theta_i)) {
304 GetMRIntProbability(theta_i, Energy);
309 GetMRIntTransProbability(theta_i, Energy);
313 G4cout <<
"Energy: " << Energy/
neV <<
"neV"
314 <<
", Enormal: " << Enormal/
neV <<
"neV" <<
G4endl;
316 G4cout <<
"Reflection_prob: " << MRpDiffuse
317 <<
", Transmission_prob: " << MRpDiffuseTrans <<
G4endl;
321 if (!
High(Enormal, FermiPotDiff)) {
329 if (
Loss(pUpScatter, theVelocityNormal, FermiPotDiff)) {
360 NewMomentum =
MRreflect(MRpDiffuse, OldMomentum, theGlobalNormal,
361 Energy, FermiPotDiff);
366 NewMomentum =
Reflect(pDiffuse, OldMomentum, theGlobalNormal);
390 theGlobalNormal, Energy, FermiPotDiff, Enew);
407 << reflectivity <<
G4endl;
413 NewMomentum =
Reflect(pDiffuse, OldMomentum, theGlobalNormal);
425 G4double mass = -std::sqrt(theMomentumNormal*theMomentumNormal -
426 neutron_mass_c2*2.*FermiPotDiff);
431 theNeutronMomentum - (theMomentumNormal-mass)*theGlobalNormal;
442 G4cout <<
"Energy: " << Energy/
neV <<
"neV, Enormal: "
443 << Enormal/
neV <<
"neV, fpdiff " << FermiPotDiff/
neV
444 <<
"neV, Enew " << Enew/
neV <<
"neV" <<
G4endl;
445 G4cout <<
"UCNBoundaryProcess: Transmit and set the new Energy "
474 G4double vBound = std::sqrt(2.*FermiPot/neutron_mass_c2*c_squared);
475 G4double vRatio = theVelocityNormal/vBound;
477 G4double pLoss = (2*pUpScatter*vRatio)/(std::sqrt(1-(vRatio*vRatio)));
483 const G4double hdm = hbar_Planck*c_squared/neutron_mass_c2;
489 pLoss *= 1+2*b*b*vBound*vBound/
490 (hdm*hdm+0.85*hdm*vBound*w+2*vBound*vBound*w*
w);
504 G4double r = (std::sqrt(Enormal) - std::sqrt(Enormal - FermiPot)) /
505 (std::sqrt(Enormal) + std::sqrt(Enormal - FermiPot));
514 G4double PdotN = OldMomentum * Normal;
521 if (NewMomentum == OldMomentum ||
G4UniformRand() < pDiffuse) {
560 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
572 G4double PdotN = OldMomentum * Normal;
574 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
596 G4double costheta = OldMomentum*Normal;
598 G4double Enormal = Energy * (costheta*costheta);
602 (1.-pDiffuse-pDiffuseTrans-pLoss);
608 if (decide < pSpecular) {
612 G4double PdotN = OldMomentum * Normal;
613 NewMomentum = OldMomentum - (2.*PdotN)*Normal;
625 if (decide < pSpecular + pDiffuse) {
632 MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
635 <<
", " << NewMomentum <<
G4endl;
645 if (decide < pSpecular + pDiffuse + pDiffuseTrans) {
652 MRDiffTrans(Normal, Energy, FermiPot, OldMomentum, pDiffuseTrans);
654 Enew = Energy - FermiPot;
663 if (decide < pSpecular + pDiffuse + pDiffuseTrans + pLoss) {
678 Enew = Energy - FermiPot;
680 G4double palt = std::sqrt(2*neutron_mass_c2/c_squared*Energy);
681 G4double produ = OldMomentum * Normal;
683 NewMomentum = palt*OldMomentum-
684 (palt*produ+std::sqrt(palt*palt*produ*produ-2*neutron_mass_c2/
685 c_squared*FermiPot))*Normal;
691 return NewMomentum.unit();
706 G4double theta_i = OldMomentum.polarAngle(-Normal);
721 GetMRMaxProbability(theta_i, Energy)/pDiffuse <=
723 GetMRProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/pDiffuse)
730 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
732 GetMRMaxProbability(theta_i, Energy)) > 1) {
733 G4cout <<
"MRMax Wahrscheinlichkeitsueberschreitung!" <<
G4endl;
735 GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
737 GetMRMaxProbability(theta_i, Energy)) << G4endl;
739 SetMRMaxProbability(theta_i, Energy,
741 GetMRProbability(theta_i, Energy,
742 FermiPot, theta_o, phi_o));
745 if(++count > 10000) { accepted =
true; }
751 localmomentum.setRThetaPhi(1., theta_o, phi_o);
767 if (momentum * Normal<0) {
770 G4cout <<
"G4UCNBoundaryProcess::MRDiffRefl: !" <<
G4endl;
773 return momentum.unit();
788 G4double theta_i = OldMomentum.polarAngle(-Normal);
803 GetMRMaxTransProbability(theta_i, Energy)/pDiffuseTrans <=
805 GetMRTransProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/
813 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
815 GetMRMaxTransProbability(theta_i, Energy)) > 1) {
816 G4cout <<
"MRMaxTrans Wahrscheinlichkeitsueberschreitung!" <<
G4endl;
818 GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
820 GetMRMaxTransProbability(theta_i, Energy)) << G4endl;
822 SetMRMaxTransProbability(theta_i, Energy,
824 GetMRTransProbability(theta_i, Energy,
825 FermiPot, theta_o, phi_o));
828 if(++count > 10000) { accepted =
true; }
834 localmomentum.setRThetaPhi(1.,
pi-theta_o, phi_o);
845 if (momentum*Normal<0) {
848 G4cout <<
"G4UCNBoundaryProcess::MRDiffTrans: !" <<
G4endl;
851 return momentum.unit();
856 return (Energy - FermiPot);
866 momentum.rotateUz(Normal);
868 if (momentum*Normal < 0) {
870 G4cout <<
"G4UCNBoundaryProcess::LDiffRefl: !" <<
G4endl;
873 return momentum.unit();
889 es1 = direction.perpPart(Normal);
898 es2.rotate(90.*
degree, es3);
919 G4cout <<
" *** No G4UCNMaterialPropertiesTable *** " <<
G4endl;
921 G4cout <<
" *** No MicroRoughness Table *** " <<
G4endl;
923 G4cout <<
" *** MicroRoughness Condition not satisfied *** " <<
G4endl;
935 G4cout <<
" *** MR Model Diffuse Reflection *** " <<
G4endl;
939 G4cout <<
" *** MR Model Diffuse Transmission *** " <<
G4endl;
948 G4cout <<
"Sum NoMRCondition: "
950 G4cout <<
"Sum No. E < V Loss: "
952 G4cout <<
"Sum No. E > V Ezero: "
954 G4cout <<
"Sum No. E < V SpinFlip: "
956 G4cout <<
"Sum No. E > V Specular Reflection: "
958 G4cout <<
"Sum No. E < V Specular Reflection: "
960 G4cout <<
"Sum No. E < V Lambertian Reflection: "
962 G4cout <<
"Sum No. E > V MR DiffuseReflection: "
964 G4cout <<
"Sum No. E < V MR DiffuseReflection: "
966 G4cout <<
"Sum No. E > V SnellTransmit: "
968 G4cout <<
"Sum No. E > V MR SnellTransmit: "
970 G4cout <<
"Sum No. E > V DiffuseTransmit: "
982 if (sd)
return sd->
Hit(&aStep);
G4int bLambertianReflection
G4double condition(const G4ErrorSymMatrix &m)
G4int aMRDiffuseReflection
G4ThreeVector MRDiffTrans(G4ThreeVector, G4double, G4double, G4ThreeVector, G4double)
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
G4double GetVelocity() const
CLHEP::HepRotation G4RotationMatrix
const G4DynamicParticle * GetDynamicParticle() const
G4StepStatus GetStepStatus() const
G4Material * GetMaterial() const
G4RotationMatrix GetCoordinateTransformMatrix(G4ThreeVector, G4ThreeVector)
G4double GetSurfaceTolerance() const
const G4double w[NPOINTSGL]
G4bool High(G4double, G4double)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4UCNBoundaryProcessMessenger * fMessenger
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4double Transmit(G4double, G4double)
G4ThreeVector MRDiffRefl(G4ThreeVector, G4double, G4double, G4ThreeVector, G4double)
G4int aSpecularReflection
void BoundaryProcessVerbose() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4StepPoint * GetPreStepPoint() const
static G4int GetHypNavigatorID()
G4double GetKineticEnergy() const
G4bool DoMicroRoughnessReflection
G4GLOB_DLL std::ostream G4cout
G4VPhysicalVolume * GetPhysicalVolume() const
G4UCNBoundaryProcess(const G4String &processName="UCNBoundaryProcess", G4ProcessType type=fUCN)
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
G4int bSpecularReflection
const G4ThreeVector & GetMomentumDirection() const
G4bool UseMicroRoughnessReflection
G4ThreeVector LDiffRefl(G4ThreeVector)
G4bool Hit(G4Step *aStep)
void SetProcessSubType(G4int)
G4int bMRDiffuseReflection
const G4String & GetProcessName() const
void BoundaryProcessSummary() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4TransportationManager * GetTransportationManager()
G4double GetCorrLen() const
G4bool Loss(G4double, G4double, G4double)
G4UCNMaterialPropertiesTable * aMaterialPropertiesTable1
G4ThreeVector MRreflect(G4double, G4ThreeVector, G4ThreeVector, G4double, G4double)
G4ThreeVector GetMomentum() const
virtual void Initialize(const G4Track &)
const G4ThreeVector & GetMomentumDirection() const
G4bool InvokeSD(const G4Step *step)
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *condition)
G4ThreeVector MRreflectHigh(G4double, G4double, G4double, G4ThreeVector, G4ThreeVector, G4double, G4double, G4double &)
G4double * GetMicroRoughnessTable()
G4StepPoint * GetPostStepPoint() const
void AddTotalEnergyDeposit(G4double value)
const G4ThreeVector & GetPolarization() const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
static const double degree
G4double GetEnergy() const
G4ThreeVector Reflect(G4double, G4ThreeVector, G4ThreeVector)
G4VSensitiveDetector * GetSensitiveDetector() const
virtual ~G4UCNBoundaryProcess()
G4UCNMaterialPropertiesTable * aMaterialPropertiesTable2
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeTrackStatus(G4TrackStatus status)
G4UCNBoundaryProcessStatus theStatus
G4Track * GetTrack() const
G4double GetConstProperty(const char *key)
void ProposeVelocity(G4double finalVelocity)
G4double Reflectivity(G4double, G4double)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4bool SpinFlip(G4double)
G4double GetStepLength() const
static G4GeometryTolerance * GetInstance()
static const G4Step * GetHyperStep()