Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UCNBoundaryProcess.cc
Go to the documentation of this file.
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 // $Id: G4UCNBoundaryPrcess.cc 69576 2013-05-08 13:48:13Z gcosmo $
27 //
29 // UCN BoundaryProcess Class Implementation
31 //
32 // File: G4UCNBoundaryProcess.cc
33 // Description: Discrete Process -- Boundary Process of UCN
34 // Version: 1.0
35 // Created: 2014-05-12
36 // Author: Peter Gumplinger
37 // adopted from Geant4UCN by Peter Fierlinger (9.9.04) and
38 // Marcin Kuzniak (21.4.06)
39 // 1/v energy dependent absorption cross section
40 // inside materials
41 // Updated: 2007 Extensions for the microroughness model by Stefan Heule
42 //
43 // mail: gum@triumf.ca
44 //
46 
47 #include "G4UCNProcessSubType.hh"
48 
49 
50 #include "G4UCNBoundaryProcess.hh"
52 
53 #include "G4GeometryTolerance.hh"
54 
55 #include "G4StepPoint.hh"
56 #include "G4ParticleDefinition.hh"
57 
59 
62 
63 #include "G4VSensitiveDetector.hh"
64 
65 #include "G4SystemOfUnits.hh"
66 #include "G4PhysicalConstants.hh"
67 
69  G4ProcessType type)
70  : G4VDiscreteProcess(processName, type)
71 {
72 
73  if (verboseLevel > 0) G4cout << GetProcessName() << " is created " << G4endl;
74 
76 
77  theStatus = Undefined;
78 
79  fMessenger = new G4UCNBoundaryProcessMessenger(this);
80 
81  neV = 1.0e-9*eV;
82 
83  kCarTolerance = G4GeometryTolerance::GetInstance()
85 
86  Material1 = NULL;
87  Material2 = NULL;
88 
89  aMaterialPropertiesTable1 = NULL;
90  aMaterialPropertiesTable2 = NULL;
91 
92  UseMicroRoughnessReflection = false;
93  DoMicroRoughnessReflection = false;
94 
95  nNoMPT = nNoMRT = nNoMRCondition = 0;
96  nAbsorption = nEzero = nFlip = 0;
97  aSpecularReflection = bSpecularReflection = 0;
98  bLambertianReflection = 0;
99  aMRDiffuseReflection = bMRDiffuseReflection = 0;
100  nSnellTransmit = mSnellTransmit = 0;
101  aMRDiffuseTransmit = 0;
102  ftheta_o = fphi_o = 0;
103 }
104 
106 {
107  delete fMessenger;
108 }
109 
112 {
113  aParticleChange.Initialize(aTrack);
115 
116  // Get hyperStep from G4ParallelWorldProcess
117  // NOTE: PostSetpDoIt of this process should be
118  // invoked after G4ParallelWorldProcess!
119 
120  const G4Step* pStep = &aStep;
121 
123 
124  if (hStep) pStep = hStep;
125 
126  G4bool isOnBoundary =
128 
129  if (isOnBoundary) {
130  Material1 = pStep->GetPreStepPoint()->GetMaterial();
131  Material2 = pStep->GetPostStepPoint()->GetMaterial();
132  } else {
133  theStatus = NotAtBoundary;
134  if ( verboseLevel > 1 ) BoundaryProcessVerbose();
135  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
136  }
137 
138  if (aTrack.GetStepLength()<=kCarTolerance/2) {
139  theStatus = StepTooSmall;
140  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
141  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
142  }
143 
144  aMaterialPropertiesTable1 = (G4UCNMaterialPropertiesTable*)Material1->
145  GetMaterialPropertiesTable();
146  aMaterialPropertiesTable2 = (G4UCNMaterialPropertiesTable*)Material2->
147  GetMaterialPropertiesTable();
148 
149  G4String volnam1 = pStep->GetPreStepPoint() ->GetPhysicalVolume()->GetName();
150  G4String volnam2 = pStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
151 
152  if (verboseLevel > 0) {
153  G4cout << " UCN at Boundary! " << G4endl;
154  G4cout << " vol1: " << volnam1 << ", vol2: " << volnam2 << G4endl;
155  G4cout << " Ekin: " << aTrack.GetKineticEnergy()/neV <<"neV"<< G4endl;
156  G4cout << " MomDir: " << aTrack.GetMomentumDirection() << G4endl;
157  }
158 
159  if (Material1 == Material2) {
160  theStatus = SameMaterial;
161  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
162  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
163  }
164 
165  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
166 
167  G4bool valid;
168  // Use the new method for Exit Normal in global coordinates,
169  // which provides the normal more reliably.
170 
171  // ID of Navigator which limits step
172 
174  std::vector<G4Navigator*>::iterator iNav =
176  GetActiveNavigatorsIterator();
177 
178  G4ThreeVector theGlobalNormal =
179  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
180 
181  if (valid) {
182  theGlobalNormal = -theGlobalNormal;
183  }
184  else
185  {
187  ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
188  << " The Navigator reports that it returned an invalid normal"
189  << G4endl;
190  G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun01",
192  "Invalid Surface Normal - Geometry must return valid surface normal");
193  }
194 
195  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
196 
197  G4ThreeVector OldMomentum = aParticle->GetMomentumDirection();
198 
199  if (OldMomentum * theGlobalNormal > 0.0) {
200 #ifdef G4OPTICAL_DEBUG
202  ed << " G4UCNBoundaryProcess/PostStepDoIt(): "
203  << " theGlobalNormal points in a wrong direction. "
204  << G4endl;
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;
211  ed << G4endl;
212  G4Exception("G4UCNBoundaryProcess::PostStepDoIt", "UCNBoun02",
213  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
214  ed,
215  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
216 #else
217  theGlobalNormal = -theGlobalNormal;
218 #endif
219  }
220 
221  G4ThreeVector theNeutronMomentum = aTrack.GetMomentum();
222 
223  G4double theMomentumNormal = theNeutronMomentum*theGlobalNormal;
224  G4double theVelocityNormal = aTrack.GetVelocity() *
225  (OldMomentum * theGlobalNormal);
226 
227  G4double Enormal = theMomentumNormal*theMomentumNormal/2./neutron_mass_c2;
228  G4double Energy = aTrack.GetKineticEnergy();
229 
230  G4double FermiPot2 = 0.;
231  G4double pDiffuse = 0.;
232  G4double pSpinFlip = 0.;
233  G4double pUpScatter = 0.;
234 
235  if (aMaterialPropertiesTable2) {
236  FermiPot2 = aMaterialPropertiesTable2->GetConstProperty("FERMIPOT")*neV;
237  pDiffuse = aMaterialPropertiesTable2->GetConstProperty("DIFFUSION");
238  pSpinFlip = aMaterialPropertiesTable2->GetConstProperty("SPINFLIP");
239  pUpScatter = aMaterialPropertiesTable2->GetConstProperty("LOSS");
240  }
241 
242  G4double FermiPot1 = 0.;
243  if (aMaterialPropertiesTable1)
244  FermiPot1 = aMaterialPropertiesTable1->GetConstProperty("FERMIPOT")*neV;
245 
246  G4double FermiPotDiff = FermiPot2 - FermiPot1;
247 
248  if ( verboseLevel > 1 )
249  G4cout << "UCNBoundaryProcess: new FermiPot: " << FermiPot2/neV
250  << "neV, old FermiPot:" << FermiPot1/neV << "neV" << G4endl;
251 
252  // Use microroughness diffuse reflection behavior if activated
253 
254  DoMicroRoughnessReflection = UseMicroRoughnessReflection;
255 
256  G4double theta_i = 0.;
257 
258  if (!aMaterialPropertiesTable2) {
259 
260  nNoMPT++;
261  theStatus = NoMPT;
262  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
263  DoMicroRoughnessReflection = false;
264 
265  } else {
266 
267  if (!aMaterialPropertiesTable2->GetMicroRoughnessTable()) {
268 
269  nNoMRT++;
270  theStatus = NoMRT;
271  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
272 
273  DoMicroRoughnessReflection = false;
274  }
275 
276  // Angle theta_in between surface and momentum direction,
277  // Phi_in is defined to be 0
278 
279  theta_i = OldMomentum.angle(-theGlobalNormal);
280 
281  // Checks the MR-conditions
282 
283  if (!aMaterialPropertiesTable2->
284  ConditionsValid(Energy, FermiPotDiff, theta_i)) {
285 
286  nNoMRCondition++;
287  theStatus = NoMRCondition;
288  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
289 
290  DoMicroRoughnessReflection = false;
291  }
292  }
293 
294  G4double MRpDiffuse = 0.;
295  G4double MRpDiffuseTrans = 0.;
296 
297  // If microroughness is available and active for material in question
298 
299  if (DoMicroRoughnessReflection) {
300 
301  // Integral probability for non-specular reflection with microroughness
302 
303  MRpDiffuse = aMaterialPropertiesTable2->
304  GetMRIntProbability(theta_i, Energy);
305 
306  // Integral probability for non-specular transmission with microroughness
307 
308  MRpDiffuseTrans = aMaterialPropertiesTable2->
309  GetMRIntTransProbability(theta_i, Energy);
310 
311  if ( verboseLevel > 1 ) {
312  G4cout << "theta: " << theta_i/degree << "degree" << G4endl;
313  G4cout << "Energy: " << Energy/neV << "neV"
314  << ", Enormal: " << Enormal/neV << "neV" << G4endl;
315  G4cout << "FermiPotDiff: " << FermiPotDiff/neV << "neV " << G4endl;
316  G4cout << "Reflection_prob: " << MRpDiffuse
317  << ", Transmission_prob: " << MRpDiffuseTrans << G4endl;
318  }
319  }
320 
321  if (!High(Enormal, FermiPotDiff)) {
322 
323  // Below critical velocity
324 
325  if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> BELOW critical velocity" << G4endl;
326 
327  // Loss on reflection
328 
329  if (Loss(pUpScatter, theVelocityNormal, FermiPotDiff)) {
330 
331  // kill it.
334 
335  nAbsorption++;
336  theStatus = Absorption;
337  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
338 
339  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
340  }
341 
342  // spinflips
343 
344  if (SpinFlip(pSpinFlip)) {
345  nFlip++;
346  theStatus = Flip;
347  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
348 
349  G4ThreeVector NewPolarization = -1. * aParticle->GetPolarization();
350  aParticleChange.ProposePolarization(NewPolarization);
351  }
352 
353  // Reflect from surface
354 
355  G4ThreeVector NewMomentum;
356 
357  // If microroughness is available and active - do non-specular reflection
358 
359  if (DoMicroRoughnessReflection)
360  NewMomentum = MRreflect(MRpDiffuse, OldMomentum, theGlobalNormal,
361  Energy, FermiPotDiff);
362  else
363 
364  // Else do it with the Lambert model as implemented by Peter Fierlinger
365 
366  NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
367 
369 
370  } else {
371 
372  // Above critical velocity
373 
374  if (verboseLevel > 0) G4cout << "G4UCNBoundaryProcess -> ABOVE critical velocity" << G4endl;
375 
376  // If it is faster than the criticial velocity,
377  // there is a probability to be still reflected.
378  // This formula is (only) valid for low loss materials
379 
380  // If microroughness available and active, do reflection with it
381 
382  G4ThreeVector NewMomentum;
383 
384  if (DoMicroRoughnessReflection) {
385 
386  G4double Enew;
387 
388  NewMomentum =
389  MRreflectHigh(MRpDiffuse, MRpDiffuseTrans, 0., OldMomentum,
390  theGlobalNormal, Energy, FermiPotDiff, Enew);
391 
392  if (Enew == 0.) {
395  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
396  } else {
401  }
402 
403  } else {
404 
405  G4double reflectivity = Reflectivity(FermiPotDiff, Enormal);
406 
407  if ( verboseLevel > 1 ) G4cout << "UCNBoundaryProcess: reflectivity "
408  << reflectivity << G4endl;
409 
410  if (G4UniformRand() < reflectivity) {
411 
412  // Reflect from surface
413 
414  NewMomentum = Reflect(pDiffuse, OldMomentum, theGlobalNormal);
416 
417  } else {
418 
419  // --- Transmission because it is faster than the critical velocity
420 
421  G4double Enew = Transmit(FermiPotDiff, Energy);
422 
423  // --- Change of the normal momentum component
424  // p = sqrt(2*m*Ekin)
425 
426  G4double mass = -std::sqrt(theMomentumNormal*theMomentumNormal -
427  neutron_mass_c2*2.*FermiPotDiff);
428 
429  // --- Momentum direction in new media
430 
431  NewMomentum =
432  theNeutronMomentum - (theMomentumNormal-mass)*theGlobalNormal;
433 
434  nSnellTransmit++;
435  theStatus = SnellTransmit;
436  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
437 
442 
443  if (verboseLevel > 1) {
444  G4cout << "Energy: " << Energy/neV << "neV, Enormal: "
445  << Enormal/neV << "neV, fpdiff " << FermiPotDiff/neV
446  << "neV, Enew " << Enew/neV << "neV" << G4endl;
447  G4cout << "UCNBoundaryProcess: Transmit and set the new Energy "
448  << aParticleChange.GetEnergy()/neV << "neV" << G4endl;
449  }
450  }
451  }
452  }
453 
454  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
455 }
456 
458  G4double ,
460 {
461  *condition = Forced;
462 
463  return DBL_MAX;
464 }
465 
466 G4bool G4UCNBoundaryProcess::Loss(G4double pUpScatter,
467  G4double theVelocityNormal,
468  G4double FermiPot)
469 {
470  // The surface roughness is not taken into account here.
471  // One could use e.g. ultracold neutrons, R. Golub, p.35,
472  // where mu is increased by roughness parameters sigma and
473  // omega, which are related to the height and the distance of
474  // "disturbances" on the surface
475 
476  G4double vBound = std::sqrt(2.*FermiPot/neutron_mass_c2*c_squared);
477  G4double vRatio = theVelocityNormal/vBound;
478 
479  G4double pLoss = (2*pUpScatter*vRatio)/(std::sqrt(1-(vRatio*vRatio)));
480 
481  // Check, if enhancement for surface roughness should be done
482 
483  if (DoMicroRoughnessReflection) {
484  if (aMaterialPropertiesTable2) {
486  G4double b = aMaterialPropertiesTable2->GetRMS();
487  G4double w = aMaterialPropertiesTable2->GetCorrLen();
488 
489  // cf. Golub's book p. 35, eq. 2.103
490 
491  pLoss *= std::sqrt(1+2*b*b*vBound*vBound/
492  (hdm*hdm+0.85*hdm*vBound*w+2*vBound*vBound*w*w));
493  }
494  }
495 
496  return (G4UniformRand() <= std::fabs(pLoss));
497 }
498 
499 G4bool G4UCNBoundaryProcess::SpinFlip(G4double pSpinFlip)
500 {
501  return (G4UniformRand() <= pSpinFlip);
502 }
503 
504 G4double G4UCNBoundaryProcess::Reflectivity(G4double FermiPot, G4double Enormal)
505 {
506  G4double r = (std::sqrt(Enormal) - std::sqrt(Enormal - FermiPot)) /
507  (std::sqrt(Enormal) + std::sqrt(Enormal - FermiPot));
508 
509  return r*r;
510 }
511 
512 G4ThreeVector G4UCNBoundaryProcess::Reflect(G4double pDiffuse,
513  G4ThreeVector OldMomentum,
514  G4ThreeVector Normal)
515 {
516  G4double PdotN = OldMomentum * Normal;
517 
518  G4ThreeVector NewMomentum = OldMomentum - (2.*PdotN)*Normal;
519  NewMomentum.unit();
520 
521  // Reflect diffuse
522 
523  if (NewMomentum == OldMomentum || G4UniformRand() < pDiffuse) {
524 
525  NewMomentum = LDiffRefl(Normal);
526 
527  bLambertianReflection++;
528  theStatus = LambertianReflection;
529  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
530 
531  return NewMomentum;
532  }
533 
534  // Reflect specular
535 
536  bSpecularReflection++;
537  theStatus = SpecularReflection;
538  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
539 
540  return NewMomentum;
541 }
542 
544  G4ThreeVector OldMomentum,
545  G4ThreeVector Normal,
546  G4double Energy,
547  G4double FermiPot)
548 {
549  // Only for Enormal <= VFermi
550 
551  G4ThreeVector NewMomentum;
552 
553  // Checks if the reflection should be non-specular
554 
555  if (G4UniformRand() <= pDiffuse) {
556 
557  // Reflect diffuse
558 
559  // Determines the angles of the non-specular reflection
560 
561  NewMomentum =
562  MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
563 
564  bMRDiffuseReflection++;
565  theStatus = MRDiffuseReflection;
566  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
567 
568  return NewMomentum;
569 
570  } else {
571 
572  // Reflect specular
573 
574  G4double PdotN = OldMomentum * Normal;
575 
576  NewMomentum = OldMomentum - (2.*PdotN)*Normal;
577  NewMomentum.unit();
578 
579  bSpecularReflection++;
580  theStatus = SpecularReflection;
581  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
582 
583  return NewMomentum;
584  }
585 }
586 
588  G4double pDiffuseTrans,
589  G4double pLoss,
590  G4ThreeVector OldMomentum,
591  G4ThreeVector Normal,
592  G4double Energy,
593  G4double FermiPot,
594  G4double &Enew)
595 {
596  // Only for Enormal > VFermi
597 
598  G4double costheta = OldMomentum*Normal;
599 
600  G4double Enormal = Energy * (costheta*costheta);
601 
602 // G4double pSpecular = Reflectivity(Enormal,FermiPot)*
603  G4double pSpecular = Reflectivity(FermiPot,Enormal)*
604  (1.-pDiffuse-pDiffuseTrans-pLoss);
605 
606  G4ThreeVector NewMomentum;
607 
608  G4double decide = G4UniformRand();
609 
610  if (decide < pSpecular) {
611 
612  // Reflect specularly
613 
614  G4double PdotN = OldMomentum * Normal;
615  NewMomentum = OldMomentum - (2.*PdotN)*Normal;
616  NewMomentum.unit();
617 
618  Enew = Energy;
619 
620  aSpecularReflection++;
621  theStatus = SpecularReflection;
622  if ( verboseLevel ) BoundaryProcessVerbose();
623 
624  return NewMomentum;
625  }
626 
627  if (decide < pSpecular + pDiffuse) {
628 
629  // Reflect diffusely
630 
631  // Determines the angles of the non-specular reflection
632 
633  NewMomentum =
634  MRDiffRefl(Normal, Energy, FermiPot, OldMomentum, pDiffuse);
635 
636  if (verboseLevel > 0) G4cout << "Diffuse normal " << Normal
637  << ", " << NewMomentum << G4endl;
638  Enew = Energy;
639 
640  aMRDiffuseReflection++;
641  theStatus = MRDiffuseReflection;
642  if ( verboseLevel ) BoundaryProcessVerbose();
643 
644  return NewMomentum;
645  }
646 
647  if (decide < pSpecular + pDiffuse + pDiffuseTrans) {
648 
649  // Transmit diffusely
650 
651  // Determines the angles of the non-specular transmission
652 
653  NewMomentum =
654  MRDiffTrans(Normal, Energy, FermiPot, OldMomentum, pDiffuseTrans);
655 
656  Enew = Energy - FermiPot;
657 
658  aMRDiffuseTransmit++;
659  theStatus = MRDiffuseTransmit;
660  if ( verboseLevel ) BoundaryProcessVerbose();
661 
662  return NewMomentum;
663  }
664 
665  if (decide < pSpecular + pDiffuse + pDiffuseTrans + pLoss) {
666 
667  // Loss
668 
669  Enew = 0.;
670 
671  nEzero++;
672  theStatus = Ezero;
673  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
674 
675  return NewMomentum;
676  }
677 
678  // last case: Refractive transmission
679 
680  Enew = Energy - FermiPot;
681 
682  G4double palt = std::sqrt(2*neutron_mass_c2/c_squared*Energy);
683  G4double produ = OldMomentum * Normal;
684 
685  NewMomentum = palt*OldMomentum-
686  (palt*produ+std::sqrt(palt*palt*produ*produ-2*neutron_mass_c2/
687  c_squared*FermiPot))*Normal;
688 
689  mSnellTransmit++;
690  theStatus = SnellTransmit;
691  if ( verboseLevel > 0 ) BoundaryProcessVerbose();
692 
693  return NewMomentum.unit();
694 }
695 
696 G4ThreeVector G4UCNBoundaryProcess::MRDiffRefl(G4ThreeVector Normal,
697  G4double Energy,
698  G4double FermiPot,
699  G4ThreeVector OldMomentum,
700  G4double pDiffuse)
701 {
702  G4bool accepted = false;
703 
704  G4double theta_o, phi_o;
705 
706  // Polar angle of incidence
707 
708  G4double theta_i = OldMomentum.polarAngle(-Normal);
709 
710  // Azimuthal angle of incidence
711 
712  // G4double phi_i = -OldMomentum.azimAngle(-Normal);
713 
714  // accept-reject method for MR-reflection
715 
716  G4int count = 0;
717  while (!accepted) {
718  theta_o = G4UniformRand()*pi/2;
719  phi_o = G4UniformRand()*pi*2-pi;
720  // Box over distribution is increased by 50% to ensure no value is above
721  if (1.5*G4UniformRand()*
722  aMaterialPropertiesTable2->
723  GetMRMaxProbability(theta_i, Energy)/pDiffuse <=
724  aMaterialPropertiesTable2->
725  GetMRProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/pDiffuse)
726 
727  accepted = true;
728 
729  // For the case that the box is nevertheless exceeded
730 
731  if (aMaterialPropertiesTable2->
732  GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
733  (1.5*aMaterialPropertiesTable2->
734  GetMRMaxProbability(theta_i, Energy)) > 1) {
735  G4cout << "MRMax Wahrscheinlichkeitsueberschreitung!" << G4endl;
736  G4cout << aMaterialPropertiesTable2->
737  GetMRProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
738  (1.5*aMaterialPropertiesTable2->
739  GetMRMaxProbability(theta_i, Energy)) << G4endl;
740  aMaterialPropertiesTable2->
741  SetMRMaxProbability(theta_i, Energy,
742  aMaterialPropertiesTable2->
743  GetMRProbability(theta_i, Energy,
744  FermiPot, theta_o, phi_o));
745  }
746  // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
747  if(++count > 10000) { accepted = true; }
748  }
749 
750  // Creates vector in the local coordinate system of the reflection
751 
752  G4ThreeVector localmomentum;
753  localmomentum.setRThetaPhi(1., theta_o, phi_o);
754 
755  ftheta_o = theta_o;
756  fphi_o = phi_o;
757 
758  // Get coordinate transform matrix
759 
760  G4RotationMatrix TransCoord =
761  GetCoordinateTransformMatrix(Normal, OldMomentum);
762 
763  //transfer to the coordinates of the global system
764 
765  G4ThreeVector momentum = TransCoord*localmomentum;
766 
767  //momentum.rotateUz(Normal);
768 
769  if (momentum * Normal<0) {
770  momentum*=-1;
771  // something has gone wrong...
772  G4cout << "G4UCNBoundaryProcess::MRDiffRefl: !" << G4endl;
773  }
774 
775  return momentum.unit();
776 }
777 
778 G4ThreeVector G4UCNBoundaryProcess::MRDiffTrans(G4ThreeVector Normal,
779  G4double Energy,
780  G4double FermiPot,
781  G4ThreeVector OldMomentum,
782  G4double pDiffuseTrans)
783 {
784  G4bool accepted = false;
785 
786  G4double theta_o, phi_o;
787 
788  // Polar angle of incidence
789 
790  G4double theta_i = OldMomentum.polarAngle(-Normal);
791 
792  // azimuthal angle of incidence
793 
794  // G4double phi_i = -OldMomentum.azimAngle(-Normal);
795 
796  G4int count = 0;
797  while (!accepted) {
798  theta_o = G4UniformRand()*pi/2;
799  phi_o = G4UniformRand()*pi*2-pi;
800 
801  // Box over distribution is increased by 50% to ensure no value is above
802 
803  if (1.5*G4UniformRand()*
804  aMaterialPropertiesTable2->
805  GetMRMaxTransProbability(theta_i, Energy)/pDiffuseTrans <=
806  aMaterialPropertiesTable2->
807  GetMRTransProbability(theta_i,Energy,FermiPot,theta_o,phi_o)/
808  pDiffuseTrans)
809 
810  accepted=true;
811 
812  // For the case that the box is nevertheless exceeded
813 
814  if(aMaterialPropertiesTable2->
815  GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
816  (1.5*aMaterialPropertiesTable2->
817  GetMRMaxTransProbability(theta_i, Energy)) > 1) {
818  G4cout << "MRMaxTrans Wahrscheinlichkeitsueberschreitung!" << G4endl;
819  G4cout << aMaterialPropertiesTable2->
820  GetMRTransProbability(theta_i, Energy, FermiPot, theta_o, phi_o)/
821  (1.5*aMaterialPropertiesTable2->
822  GetMRMaxTransProbability(theta_i, Energy)) << G4endl;
823  aMaterialPropertiesTable2->
824  SetMRMaxTransProbability(theta_i, Energy,
825  aMaterialPropertiesTable2->
826  GetMRTransProbability(theta_i, Energy,
827  FermiPot, theta_o, phi_o));
828  }
829  // Loop checking, 31-Aug-2015, Vladimir Ivanchenko
830  if(++count > 10000) { accepted = true; }
831  }
832 
833  // Creates vector in the local coordinate system of the reflection
834 
835  G4ThreeVector localmomentum;
836  localmomentum.setRThetaPhi(1., pi-theta_o, phi_o);
837 
838  // Get coordinate transform matrix
839 
840  G4RotationMatrix TransCoord =
841  GetCoordinateTransformMatrix(Normal, OldMomentum);
842 
843  // Transfer to the coordinates of the global system
844 
845  G4ThreeVector momentum = TransCoord*localmomentum;
846 
847  if (momentum*Normal<0) {
848  // something has gone wrong...
849  momentum*=-1;
850  G4cout << "G4UCNBoundaryProcess::MRDiffTrans: !" << G4endl;
851  }
852 
853  return momentum.unit();
854 }
855 
856 G4double G4UCNBoundaryProcess::Transmit(G4double FermiPot, G4double Energy)
857 {
858  return (Energy - FermiPot);
859 }
860 
861 G4ThreeVector G4UCNBoundaryProcess::LDiffRefl(G4ThreeVector Normal)
862 {
863  G4ThreeVector momentum;
864 
865  // cosine distribution - Lambert's law
866 
867  momentum.setRThetaPhi(1., std::acos(std::sqrt(G4UniformRand())), 2.*pi*G4UniformRand());
868  momentum.rotateUz(Normal);
869 
870  if (momentum*Normal < 0) {
871  momentum*=-1;
872  G4cout << "G4UCNBoundaryProcess::LDiffRefl: !" << G4endl;
873  }
874 
875  return momentum.unit();
876 }
877 
878 G4RotationMatrix G4UCNBoundaryProcess::
879  GetCoordinateTransformMatrix(G4ThreeVector Normal,
880  G4ThreeVector direction)
881 {
882  // Definition of the local coordinate system (c.s. of the reflection)
883 
884  G4ThreeVector es1, es2, es3;
885 
886  // z-Coordinate is the surface normal, has already length 1
887 
888  es3 = Normal;
889 
890  // Perpendicular part of incident direction w.r.t. normal
891  es1 = direction.perpPart(Normal);
892 
893  // Set to unit length: x-Coordinate
894 
895  es1.setMag(1.);
896  es2 = es1;
897 
898  // y-Coordinate is the pi/2-rotation of x-axis around z-axis
899 
900  es2.rotate(90.*degree, es3);
901 
902  // Transformation matrix consists just of the three coordinate vectors
903  // as the global coordinate system is the 'standard' coordinate system
904 
905  G4RotationMatrix matrix(es1, es2, es3);
906 
907  return matrix;
908 }
909 
910 void G4UCNBoundaryProcess::BoundaryProcessVerbose() const
911 {
912  if ( theStatus == Undefined )
913  G4cout << " *** Undefined *** " << G4endl;
914  if ( theStatus == NotAtBoundary )
915  G4cout << " *** NotAtBoundary *** " << G4endl;
916  if ( theStatus == SameMaterial )
917  G4cout << " *** SameMaterial *** " << G4endl;
918  if ( theStatus == StepTooSmall )
919  G4cout << " *** StepTooSmall *** " << G4endl;
920  if ( theStatus == NoMPT )
921  G4cout << " *** No G4UCNMaterialPropertiesTable *** " << G4endl;
922  if ( theStatus == NoMRT )
923  G4cout << " *** No MicroRoughness Table *** " << G4endl;
924  if ( theStatus == NoMRCondition )
925  G4cout << " *** MicroRoughness Condition not satisfied *** " << G4endl;
926  if ( theStatus == Absorption )
927  G4cout << " *** Loss on Surface *** " << G4endl;
928  if ( theStatus == Ezero )
929  G4cout << " *** Ezero on Surface *** " << G4endl;
930  if ( theStatus == Flip )
931  G4cout << " *** Spin Flip on Surface *** " << G4endl;
932  if ( theStatus == SpecularReflection )
933  G4cout << " *** Specular Reflection *** " << G4endl;
934  if ( theStatus == LambertianReflection )
935  G4cout << " *** LambertianR Reflection *** " << G4endl;
936  if ( theStatus == MRDiffuseReflection )
937  G4cout << " *** MR Model Diffuse Reflection *** " << G4endl;
938  if ( theStatus == SnellTransmit )
939  G4cout << " *** Snell Transmission *** " << G4endl;
940  if ( theStatus == MRDiffuseTransmit )
941  G4cout << " *** MR Model Diffuse Transmission *** " << G4endl;
942 }
943 
945 {
946  G4cout << "Sum NoMT: "
947  << nNoMPT << G4endl;
948  G4cout << "Sum NoMRT: "
949  << nNoMRT << G4endl;
950  G4cout << "Sum NoMRCondition: "
951  << nNoMRCondition << G4endl;
952  G4cout << "Sum No. E < V Loss: "
953  << nAbsorption << G4endl;
954  G4cout << "Sum No. E > V Ezero: "
955  << nEzero << G4endl;
956  G4cout << "Sum No. E < V SpinFlip: "
957  << nFlip << G4endl;
958  G4cout << "Sum No. E > V Specular Reflection: "
959  << aSpecularReflection << G4endl;
960  G4cout << "Sum No. E < V Specular Reflection: "
961  << bSpecularReflection << G4endl;
962  G4cout << "Sum No. E < V Lambertian Reflection: "
963  << bLambertianReflection << G4endl;
964  G4cout << "Sum No. E > V MR DiffuseReflection: "
965  << aMRDiffuseReflection << G4endl;
966  G4cout << "Sum No. E < V MR DiffuseReflection: "
967  << bMRDiffuseReflection << G4endl;
968  G4cout << "Sum No. E > V SnellTransmit: "
969  << nSnellTransmit << G4endl;
970  G4cout << "Sum No. E > V MR SnellTransmit: "
971  << mSnellTransmit << G4endl;
972  G4cout << "Sum No. E > V DiffuseTransmit: "
973  << aMRDiffuseTransmit << G4endl;
974  G4cout << " " << G4endl;
975 }
976 
977 G4bool G4UCNBoundaryProcess::InvokeSD(const G4Step* pStep)
978 {
979  G4Step aStep = *pStep;
980 
982 
984  if (sd) return sd->Hit(&aStep);
985  else return false;
986 }
G4double condition(const G4ErrorSymMatrix &m)
G4int verboseLevel
Definition: G4VProcess.hh:368
double angle(const Hep3Vector &) const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetVelocity() const
const G4DynamicParticle * GetDynamicParticle() const
G4StepStatus GetStepStatus() const
G4Material * GetMaterial() const
G4double GetSurfaceTolerance() const
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
int G4int
Definition: G4Types.hh:78
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4StepPoint * GetPreStepPoint() const
tuple b
Definition: test.py:12
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4VPhysicalVolume * GetPhysicalVolume() const
G4UCNBoundaryProcess(const G4String &processName="UCNBoundaryProcess", G4ProcessType type=fUCN)
const G4String & GetName() const
static constexpr double degree
Definition: G4SIunits.hh:144
const G4ThreeVector & GetPosition() const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
static constexpr double eV
Definition: G4SIunits.hh:215
G4bool Hit(G4Step *aStep)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Hep3Vector perpPart() const
Definition: G4Step.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void setRThetaPhi(double r, double theta, double phi)
void BoundaryProcessSummary() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
float neutron_mass_c2
Definition: hepunit.py:276
G4ThreeVector MRreflect(G4double, G4ThreeVector, G4ThreeVector, G4double, G4double)
G4ThreeVector GetMomentum() const
virtual void Initialize(const G4Track &)
const G4ThreeVector & GetMomentumDirection() const
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *condition)
G4ThreeVector MRreflectHigh(G4double, G4double, G4double, G4ThreeVector, G4ThreeVector, G4double, G4double, G4double &)
Hep3Vector unit() const
G4StepPoint * GetPostStepPoint() const
void AddTotalEnergyDeposit(G4double value)
const G4ThreeVector & GetPolarization() const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
G4double GetEnergy() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
G4VSensitiveDetector * GetSensitiveDetector() const
void setMag(double)
Definition: ThreeVector.cc:25
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4Track * GetTrack() const
void ProposeVelocity(G4double finalVelocity)
#define DBL_MAX
Definition: templates.hh:83
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
float c_light
Definition: hepunit.py:257
G4double GetStepLength() const
static G4GeometryTolerance * GetInstance()
static const G4Step * GetHyperStep()
G4double GetConstProperty(const char *key) const
G4ProcessType
double polarAngle(const Hep3Vector &v2) const
Definition: SpaceVectorD.cc:28