ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UCNBoundaryProcess.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UCNBoundaryProcess.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
28 // UCN BoundaryProcess Class Implementation
30 //
31 // File: G4UCNBoundaryProcess.cc
32 // Description: Discrete Process -- Boundary Process of UCN
33 // Version: 1.0
34 // Created: 2014-05-12
35 // Author: Peter Gumplinger
36 // adopted from Geant4UCN by Peter Fierlinger (9.9.04) and
37 // Marcin Kuzniak (21.4.06)
38 // 1/v energy dependent absorption cross section
39 // inside materials
40 // Updated: 2007 Extensions for the microroughness model by Stefan Heule
41 //
42 // mail: gum@triumf.ca
43 //
45 
46 #include "G4UCNProcessSubType.hh"
47 
48 
49 #include "G4UCNBoundaryProcess.hh"
51 
52 #include "G4GeometryTolerance.hh"
53 
54 #include "G4StepPoint.hh"
55 #include "G4ParticleDefinition.hh"
56 
58 
61 
62 #include "G4VSensitiveDetector.hh"
63 
64 #include "G4SystemOfUnits.hh"
65 #include "G4PhysicalConstants.hh"
66 
68  G4ProcessType type)
69  : G4VDiscreteProcess(processName, type)
70 {
71 
72  if (verboseLevel > 0) G4cout << GetProcessName() << " is created " << G4endl;
73 
75 
77 
79 
80  neV = 1.0e-9*eV;
81 
84 
85  Material1 = NULL;
86  Material2 = NULL;
87 
90 
93 
95  nAbsorption = nEzero = nFlip = 0;
100  aMRDiffuseTransmit = 0;
101  ftheta_o = fphi_o = 0;
102 }
103 
105 {
106  delete fMessenger;
107 }
108 
111 {
112  aParticleChange.Initialize(aTrack);
114 
115  // Get hyperStep from G4ParallelWorldProcess
116  // NOTE: PostSetpDoIt of this process should be
117  // invoked after G4ParallelWorldProcess!
118 
119  const G4Step* pStep = &aStep;
120 
122 
123  if (hStep) pStep = hStep;
124 
125  G4bool isOnBoundary =
127 
128  if (isOnBoundary) {
129  Material1 = pStep->GetPreStepPoint()->GetMaterial();
130  Material2 = pStep->GetPostStepPoint()->GetMaterial();
131  } else {
134  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
135  }
136 
137  if (aTrack.GetStepLength()<=kCarTolerance/2) {
140  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
141  }
142 
144  GetMaterialPropertiesTable();
146  GetMaterialPropertiesTable();
147 
148  G4String volnam1 = pStep->GetPreStepPoint() ->GetPhysicalVolume()->GetName();
149  G4String volnam2 = pStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
150 
151  if (verboseLevel > 0) {
152  G4cout << " UCN at Boundary! " << G4endl;
153  G4cout << " vol1: " << volnam1 << ", vol2: " << volnam2 << G4endl;
154  G4cout << " Ekin: " << aTrack.GetKineticEnergy()/neV <<"neV"<< G4endl;
155  G4cout << " MomDir: " << aTrack.GetMomentumDirection() << G4endl;
156  }
157 
158  if (Material1 == Material2) {
161  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
162  }
163 
164  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
165 
166  G4bool valid;
167  // Use the new method for Exit Normal in global coordinates,
168  // which provides the normal more reliably.
169 
170  // ID of Navigator which limits step
171 
173  std::vector<G4Navigator*>::iterator iNav =
175  GetActiveNavigatorsIterator();
176 
177  G4ThreeVector theGlobalNormal =
178  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
179 
180  if (valid) {
181  theGlobalNormal = -theGlobalNormal;
182  }
183  else
184  {
186  ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
187  << " The Navigator reports that it returned an invalid normal"
188  << G4endl;
189  G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun01",
191  "Invalid Surface Normal - Geometry must return valid surface normal");
192  }
193 
194  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
195 
196  G4ThreeVector OldMomentum = aParticle->GetMomentumDirection();
197 
198  if (OldMomentum * theGlobalNormal > 0.0) {
199 #ifdef G4OPTICAL_DEBUG
201  ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
202  << " theGlobalNormal points in a wrong direction. "
203  << G4endl;
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;
210  ed << G4endl;
211  G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun02",
212  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
213  ed,
214  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
215 #else
216  theGlobalNormal = -theGlobalNormal;
217 #endif
218  }
219 
220  G4ThreeVector theNeutronMomentum = aTrack.GetMomentum();
221 
222  G4double theMomentumNormal = theNeutronMomentum*theGlobalNormal;
223  G4double theVelocityNormal = aTrack.GetVelocity() *
224  (OldMomentum * theGlobalNormal);
225 
226  G4double Enormal = theMomentumNormal*theMomentumNormal/2./neutron_mass_c2;
227  G4double Energy = aTrack.GetKineticEnergy();
228 
229  G4double FermiPot2 = 0.;
230  G4double pDiffuse = 0.;
231  G4double pSpinFlip = 0.;
232  G4double pUpScatter = 0.;
233 
235  FermiPot2 = aMaterialPropertiesTable2->GetConstProperty("FERMIPOT")*neV;
236  pDiffuse = aMaterialPropertiesTable2->GetConstProperty("DIFFUSION");
237  pSpinFlip = aMaterialPropertiesTable2->GetConstProperty("SPINFLIP");
238  pUpScatter = aMaterialPropertiesTable2->GetConstProperty("LOSS");
239  }
240 
241  G4double FermiPot1 = 0.;
243  FermiPot1 = aMaterialPropertiesTable1->GetConstProperty("FERMIPOT")*neV;
244 
245  G4double FermiPotDiff = FermiPot2 - FermiPot1;
246 
247  if ( verboseLevel > 1 )
248  G4cout << "UCNBoundaryProcess: new FermiPot: " << FermiPot2/neV
249  << "neV, old FermiPot:" << FermiPot1/neV << "neV" << G4endl;
250 
251  // Use microroughness diffuse reflection behavior if activated
252 
254 
255  G4double theta_i = 0.;
256 
258 
259  nNoMPT++;
260  theStatus = NoMPT;
263 
264  } else {
265 
267 
268  nNoMRT++;
269  theStatus = NoMRT;
271 
273  }
274 
275  // Angle theta_in between surface and momentum direction,
276  // Phi_in is defined to be 0
277 
278  theta_i = OldMomentum.angle(-theGlobalNormal);
279 
280  // Checks the MR-conditions
281 
283  ConditionsValid(Energy, FermiPotDiff, theta_i)) {
284 
285  nNoMRCondition++;
288 
290  }
291  }
292 
293  G4double MRpDiffuse = 0.;
294  G4double MRpDiffuseTrans = 0.;
295 
296  // If microroughness is available and active for material in question
297 
299 
300  // Integral probability for non-specular reflection with microroughness
301 
302  MRpDiffuse = aMaterialPropertiesTable2->
303  GetMRIntProbability(theta_i, Energy);
304 
305  // Integral probability for non-specular transmission with microroughness
306 
307  MRpDiffuseTrans = aMaterialPropertiesTable2->
308  GetMRIntTransProbability(theta_i, Energy);
309 
310  if ( verboseLevel > 1 ) {
311  G4cout << "theta: " << theta_i/degree << "degree" << G4endl;
312  G4cout << "Energy: " << Energy/neV << "neV"
313  << ", Enormal: " << Enormal/neV << "neV" << G4endl;
314  G4cout << "FermiPotDiff: " << FermiPotDiff/neV << "neV " << G4endl;
315  G4cout << "Reflection_prob: " << MRpDiffuse
316  << ", Transmission_prob: " << MRpDiffuseTrans << G4endl;
317  }
318  }
319 
320  if (!High(Enormal, FermiPotDiff)) {
321 
322  // Below critical velocity
323 
324  if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> BELOW critical velocity" << G4endl;
325 
326  // Loss on reflection
327 
328  if (Loss(pUpScatter, theVelocityNormal, FermiPotDiff)) {
329 
330  // kill it.
333 
334  nAbsorption++;
337 
338  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
339  }
340 
341  // spinflips
342 
343  if (SpinFlip(pSpinFlip)) {
344  nFlip++;
345  theStatus = Flip;
347 
348  G4ThreeVector NewPolarization = -1. * aParticle->GetPolarization();
349  aParticleChange.ProposePolarization(NewPolarization);
350  }
351 
352  // Reflect from surface
353 
354  G4ThreeVector NewMomentum;
355 
356  // If microroughness is available and active - do non-specular reflection
357 
359  NewMomentum = MRreflect(MRpDiffuse, OldMomentum, theGlobalNormal,
360  Energy, FermiPotDiff);
361  else
362 
363  // Else do it with the Lambert model as implemented by Peter Fierlinger
364 
365  NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
366 
368 
369  } else {
370 
371  // Above critical velocity
372 
373  if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> ABOVE critical velocity" << G4endl;
374 
375  // If it is faster than the criticial velocity,
376  // there is a probability to be still reflected.
377  // This formula is (only) valid for low loss materials
378 
379  // If microroughness available and active, do reflection with it
380 
381  G4ThreeVector NewMomentum;
382 
384 
385  G4double Enew;
386 
387  NewMomentum =
388  MRreflectHigh(MRpDiffuse, MRpDiffuseTrans, 0., OldMomentum,
389  theGlobalNormal, Energy, FermiPotDiff, Enew);
390 
391  if (Enew == 0.) {
394  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
395  } else {
400  }
401 
402  } else {
403 
404  G4double reflectivity = Reflectivity(FermiPotDiff, Enormal);
405 
406  if ( verboseLevel > 1 ) G4cout << "UCNBoundaryProcess: reflectivity "
407  << reflectivity << G4endl;
408 
409  if (G4UniformRand() < reflectivity) {
410 
411  // Reflect from surface
412 
413  NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
415 
416  } else {
417 
418  // --- Transmission because it is faster than the critical velocity
419 
420  G4double Enew = Transmit(FermiPotDiff, Energy);
421 
422  // --- Change of the normal momentum component
423  // p = sqrt(2*m*Ekin)
424 
425  G4double mass = -std::sqrt(theMomentumNormal*theMomentumNormal -
426  neutron_mass_c2*2.*FermiPotDiff);
427 
428  // --- Momentum direction in new media
429 
430  NewMomentum =
431  theNeutronMomentum - (theMomentumNormal-mass)*theGlobalNormal;
432 
433  nSnellTransmit++;
436 
441 
442  if (verboseLevel > 1) {
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 "
447  << aParticleChange.GetEnergy()/neV << "neV" << G4endl;
448  }
449  }
450  }
451  }
452 
453  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
454 }
455 
457  G4double ,
459 {
460  *condition = Forced;
461 
462  return DBL_MAX;
463 }
464 
466  G4double theVelocityNormal,
467  G4double FermiPot)
468 {
469  // The surface roughness is not taken into account here.
470  // One could use e.g. ultracold neutrons, R. Golub, p.35,
471  // where mu is increased by roughness parameters sigma and
472  // omega, which are related to the height and the distance of
473  // "disturbances" on the surface
474 
475  G4double vBound = std::sqrt(2.*FermiPot/neutron_mass_c2*c_squared);
476  G4double vRatio = theVelocityNormal/vBound;
477 
478  G4double pLoss = (2*pUpScatter*vRatio)/(std::sqrt(1-(vRatio*vRatio)));
479 
480  // Check, if enhancement for surface roughness should be done
481 
484  const G4double hdm = hbar_Planck*c_squared/neutron_mass_c2;
487 
488  // cf. Golub's book p. 35, eq. 2.103
489 
490  pLoss *= std::sqrt(1+2*b*b*vBound*vBound/
491  (hdm*hdm+0.85*hdm*vBound*w+2*vBound*vBound*w*w));
492  }
493  }
494 
495  return (G4UniformRand() <= std::fabs(pLoss));
496 }
497 
499 {
500  return (G4UniformRand() <= pSpinFlip);
501 }
502 
504 {
505  G4double r = (std::sqrt(Enormal) - std::sqrt(Enormal - FermiPot)) /
506  (std::sqrt(Enormal) + std::sqrt(Enormal - FermiPot));
507 
508  return r*r;
509 }
510 
512  G4ThreeVector OldMomentum,
513  G4ThreeVector Normal)
514 {
515  G4double PdotN = OldMomentum * Normal;
516 
517  G4ThreeVector NewMomentum = OldMomentum - (2.*PdotN)*Normal;
518  NewMomentum.unit();
519 
520  // Reflect diffuse
521 
522  if (NewMomentum == OldMomentum || G4UniformRand() < pDiffuse) {
523 
524  NewMomentum = LDiffRefl(Normal);
525 
529 
530  return NewMomentum;
531  }
532 
533  // Reflect specular
534 
538 
539  return NewMomentum;
540 }
541 
543  G4ThreeVector OldMomentum,
544  G4ThreeVector Normal,
545  G4double Energy,
546  G4double FermiPot)
547 {
548  // Only for Enormal <= VFermi
549 
550  G4ThreeVector NewMomentum;
551 
552  // Checks if the reflection should be non-specular
553 
554  if (G4UniformRand() <= pDiffuse) {
555 
556  // Reflect diffuse
557 
558  // Determines the angles of the non-specular reflection
559 
560  NewMomentum =
561  MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
562 
566 
567  return NewMomentum;
568 
569  } else {
570 
571  // Reflect specular
572 
573  G4double PdotN = OldMomentum * Normal;
574 
575  NewMomentum = OldMomentum - (2.*PdotN)*Normal;
576  NewMomentum.unit();
577 
581 
582  return NewMomentum;
583  }
584 }
585 
587  G4double pDiffuseTrans,
588  G4double pLoss,
589  G4ThreeVector OldMomentum,
590  G4ThreeVector Normal,
591  G4double Energy,
592  G4double FermiPot,
593  G4double &Enew)
594 {
595  // Only for Enormal > VFermi
596 
597  G4double costheta = OldMomentum*Normal;
598 
599  G4double Enormal = Energy * (costheta*costheta);
600 
601 // G4double pSpecular = Reflectivity(Enormal,FermiPot)*
602  G4double pSpecular = Reflectivity(FermiPot,Enormal)*
603  (1.-pDiffuse-pDiffuseTrans-pLoss);
604 
605  G4ThreeVector NewMomentum;
606 
607  G4double decide = G4UniformRand();
608 
609  if (decide < pSpecular) {
610 
611  // Reflect specularly
612 
613  G4double PdotN = OldMomentum * Normal;
614  NewMomentum = OldMomentum - (2.*PdotN)*Normal;
615  NewMomentum.unit();
616 
617  Enew = Energy;
618 
622 
623  return NewMomentum;
624  }
625 
626  if (decide < pSpecular + pDiffuse) {
627 
628  // Reflect diffusely
629 
630  // Determines the angles of the non-specular reflection
631 
632  NewMomentum =
633  MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
634 
635  if (verboseLevel > 0) G4cout << "Diffuse normal " << Normal
636  << ", " << NewMomentum << G4endl;
637  Enew = Energy;
638 
642 
643  return NewMomentum;
644  }
645 
646  if (decide < pSpecular + pDiffuse + pDiffuseTrans) {
647 
648  // Transmit diffusely
649 
650  // Determines the angles of the non-specular transmission
651 
652  NewMomentum =
653  MRDiffTrans(Normal, Energy, FermiPot, OldMomentum, pDiffuseTrans);
654 
655  Enew = Energy - FermiPot;
656 
660 
661  return NewMomentum;
662  }
663 
664  if (decide < pSpecular + pDiffuse + pDiffuseTrans + pLoss) {
665 
666  // Loss
667 
668  Enew = 0.;
669 
670  nEzero++;
671  theStatus = Ezero;
673 
674  return NewMomentum;
675  }
676 
677  // last case: Refractive transmission
678 
679  Enew = Energy - FermiPot;
680 
681  G4double palt = std::sqrt(2*neutron_mass_c2/c_squared*Energy);
682  G4double produ = OldMomentum * Normal;
683 
684  NewMomentum = palt*OldMomentum-
685  (palt*produ+std::sqrt(palt*palt*produ*produ-2*neutron_mass_c2/
686  c_squared*FermiPot))*Normal;
687 
688  mSnellTransmit++;
691 
692  return NewMomentum.unit();
693 }
694 
696  G4double Energy,
697  G4double FermiPot,
698  G4ThreeVector OldMomentum,
699  G4double pDiffuse)
700 {
701  G4bool accepted = false;
702 
703  G4double theta_o, phi_o;
704 
705  // Polar angle of incidence
706 
707  G4double theta_i = OldMomentum.polarAngle(-Normal);
708 
709  // Azimuthal angle of incidence
710 
711  // G4double phi_i = -OldMomentum.azimAngle(-Normal);
712 
713  // accept-reject method for MR-reflection
714 
715  G4int count = 0;
716  while (!accepted) {
717  theta_o = G4UniformRand()*pi/2;
718  phi_o = G4UniformRand()*pi*2-pi;
719  // Box over distribution is increased by 50% to ensure no value is above
720  if (1.5*G4UniformRand()*
722  GetMRMaxProbability(theta_i, Energy)/pDiffuse <=
724  GetMRProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/pDiffuse)
725 
726  accepted = true;
727 
728  // For the case that the box is nevertheless exceeded
729 
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));
744  }
745  // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
746  if(++count > 10000) { accepted = true; }
747  }
748 
749  // Creates vector in the local coordinate system of the reflection
750 
751  G4ThreeVector localmomentum;
752  localmomentum.setRThetaPhi(1., theta_o, phi_o);
753 
754  ftheta_o = theta_o;
755  fphi_o = phi_o;
756 
757  // Get coordinate transform matrix
758 
759  G4RotationMatrix TransCoord =
760  GetCoordinateTransformMatrix(Normal, OldMomentum);
761 
762  //transfer to the coordinates of the global system
763 
764  G4ThreeVector momentum = TransCoord*localmomentum;
765 
766  //momentum.rotateUz(Normal);
767 
768  if (momentum * Normal<0) {
769  momentum*=-1;
770  // something has gone wrong...
771  G4cout << "G4UCNBoundaryProcess::MRDiffRefl: !" << G4endl;
772  }
773 
774  return momentum.unit();
775 }
776 
778  G4double Energy,
779  G4double FermiPot,
780  G4ThreeVector OldMomentum,
781  G4double pDiffuseTrans)
782 {
783  G4bool accepted = false;
784 
785  G4double theta_o, phi_o;
786 
787  // Polar angle of incidence
788 
789  G4double theta_i = OldMomentum.polarAngle(-Normal);
790 
791  // azimuthal angle of incidence
792 
793  // G4double phi_i = -OldMomentum.azimAngle(-Normal);
794 
795  G4int count = 0;
796  while (!accepted) {
797  theta_o = G4UniformRand()*pi/2;
798  phi_o = G4UniformRand()*pi*2-pi;
799 
800  // Box over distribution is increased by 50% to ensure no value is above
801 
802  if (1.5*G4UniformRand()*
804  GetMRMaxTransProbability(theta_i, Energy)/pDiffuseTrans <=
806  GetMRTransProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/
807  pDiffuseTrans)
808 
809  accepted=true;
810 
811  // For the case that the box is nevertheless exceeded
812 
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));
827  }
828  // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
829  if(++count > 10000) { accepted = true; }
830  }
831 
832  // Creates vector in the local coordinate system of the reflection
833 
834  G4ThreeVector localmomentum;
835  localmomentum.setRThetaPhi(1., pi-theta_o, phi_o);
836 
837  // Get coordinate transform matrix
838 
839  G4RotationMatrix TransCoord =
840  GetCoordinateTransformMatrix(Normal, OldMomentum);
841 
842  // Transfer to the coordinates of the global system
843 
844  G4ThreeVector momentum = TransCoord*localmomentum;
845 
846  if (momentum*Normal<0) {
847  // something has gone wrong...
848  momentum*=-1;
849  G4cout << "G4UCNBoundaryProcess::MRDiffTrans: !" << G4endl;
850  }
851 
852  return momentum.unit();
853 }
854 
856 {
857  return (Energy - FermiPot);
858 }
859 
861 {
863 
864  // cosine distribution - Lambert's law
865 
866  momentum.setRThetaPhi(1., std::acos(std::sqrt(G4UniformRand())), 2.*pi*G4UniformRand());
867  momentum.rotateUz(Normal);
868 
869  if (momentum*Normal < 0) {
870  momentum*=-1;
871  G4cout << "G4UCNBoundaryProcess::LDiffRefl: !" << G4endl;
872  }
873 
874  return momentum.unit();
875 }
876 
879  G4ThreeVector direction)
880 {
881  // Definition of the local coordinate system (c.s. of the reflection)
882 
883  G4ThreeVector es1, es2, es3;
884 
885  // z-Coordinate is the surface normal, has already length 1
886 
887  es3 = Normal;
888 
889  // Perpendicular part of incident direction w.r.t. normal
890  es1 = direction.perpPart(Normal);
891 
892  // Set to unit length: x-Coordinate
893 
894  es1.setMag(1.);
895  es2 = es1;
896 
897  // y-Coordinate is the pi/2-rotation of x-axis around z-axis
898 
899  es2.rotate(90.*degree, es3);
900 
901  // Transformation matrix consists just of the three coordinate vectors
902  // as the global coordinate system is the 'standard' coordinate system
903 
904  G4RotationMatrix matrix(es1, es2, es3);
905 
906  return matrix;
907 }
908 
910 {
911  if ( theStatus == Undefined )
912  G4cout << " *** Undefined *** " << G4endl;
913  if ( theStatus == NotAtBoundary )
914  G4cout << " *** NotAtBoundary *** " << G4endl;
915  if ( theStatus == SameMaterial )
916  G4cout << " *** SameMaterial *** " << G4endl;
917  if ( theStatus == StepTooSmall )
918  G4cout << " *** StepTooSmall *** " << G4endl;
919  if ( theStatus == NoMPT )
920  G4cout << " *** No G4UCNMaterialPropertiesTable *** " << G4endl;
921  if ( theStatus == NoMRT )
922  G4cout << " *** No MicroRoughness Table *** " << G4endl;
923  if ( theStatus == NoMRCondition )
924  G4cout << " *** MicroRoughness Condition not satisfied *** " << G4endl;
925  if ( theStatus == Absorption )
926  G4cout << " *** Loss on Surface *** " << G4endl;
927  if ( theStatus == Ezero )
928  G4cout << " *** Ezero on Surface *** " << G4endl;
929  if ( theStatus == Flip )
930  G4cout << " *** Spin Flip on Surface *** " << G4endl;
931  if ( theStatus == SpecularReflection )
932  G4cout << " *** Specular Reflection *** " << G4endl;
934  G4cout << " *** LambertianR Reflection *** " << G4endl;
936  G4cout << " *** MR Model Diffuse Reflection *** " << G4endl;
937  if ( theStatus == SnellTransmit )
938  G4cout << " *** Snell Transmission *** " << G4endl;
939  if ( theStatus == MRDiffuseTransmit )
940  G4cout << " *** MR Model Diffuse Transmission *** " << G4endl;
941 }
942 
944 {
945  G4cout << "Sum NoMT: "
946  << nNoMPT << G4endl;
947  G4cout << "Sum NoMRT: "
948  << nNoMRT << G4endl;
949  G4cout << "Sum NoMRCondition: "
950  << nNoMRCondition << G4endl;
951  G4cout << "Sum No. E < V Loss: "
952  << nAbsorption << G4endl;
953  G4cout << "Sum No. E > V Ezero: "
954  << nEzero << G4endl;
955  G4cout << "Sum No. E < V SpinFlip: "
956  << nFlip << G4endl;
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: "
968  << nSnellTransmit << G4endl;
969  G4cout << "Sum No. E > V MR SnellTransmit: "
970  << mSnellTransmit << G4endl;
971  G4cout << "Sum No. E > V DiffuseTransmit: "
973  G4cout << " " << G4endl;
974 }
975 
977 {
978  G4Step aStep = *pStep;
979 
981 
983  if (sd) return sd->Hit(&aStep);
984  else return false;
985 }