Geant4  10.03
G4OpBoundaryProcess.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 //
27 // Optical Photon Boundary Process Class Implementation
29 //
30 // File: G4OpBoundaryProcess.cc
31 // Description: Discrete Process -- reflection/refraction at
32 // optical interfaces
33 // Version: 1.1
34 // Created: 1997-06-18
35 // Modified: 1998-05-25 - Correct parallel component of polarization
36 // (thanks to: Stefano Magni + Giovanni Pieri)
37 // 1998-05-28 - NULL Rindex pointer before reuse
38 // (thanks to: Stefano Magni)
39 // 1998-06-11 - delete *sint1 in oblique reflection
40 // (thanks to: Giovanni Pieri)
41 // 1998-06-19 - move from GetLocalExitNormal() to the new
42 // method: GetLocalExitNormal(&valid) to get
43 // the surface normal in all cases
44 // 1998-11-07 - NULL OpticalSurface pointer before use
45 // comparison not sharp for: std::abs(cost1) < 1.0
46 // remove sin1, sin2 in lines 556,567
47 // (thanks to Stefano Magni)
48 // 1999-10-10 - Accommodate changes done in DoAbsorption by
49 // changing logic in DielectricMetal
50 // 2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
51 // might be used uninitialized in this function
52 // moved E2_perp, E2_parl and E2_total out of 'if'
53 // 2003-11-27 - Modified line 168-9 to reflect changes made to
54 // G4OpticalSurface class ( by Fan Lei)
55 // 2004-02-02 - Set theStatus = Undefined at start of DoIt
56 // 2005-07-28 - add G4ProcessType to constructor
57 // 2006-11-04 - add capability of calculating the reflectivity
58 // off a metal surface by way of a complex index
59 // of refraction - Thanks to Sehwook Lee and John
60 // Hauptman (Dept. of Physics - Iowa State Univ.)
61 // 2009-11-10 - add capability of simulating surface reflections
62 // with Look-Up-Tables (LUT) containing measured
63 // optical reflectance for a variety of surface
64 // treatments - Thanks to Martin Janecek and
65 // William Moses (Lawrence Berkeley National Lab.)
66 // 2013-06-01 - add the capability of simulating the transmission
67 // of a dichronic filter
68 //
69 // Author: Peter Gumplinger
70 // adopted from work by Werner Keil - April 2/96
71 // mail: gum@triumf.ca
72 //
74 
75 #include "G4ios.hh"
76 #include "G4PhysicalConstants.hh"
77 #include "G4OpProcessSubType.hh"
78 
79 #include "G4OpBoundaryProcess.hh"
80 #include "G4GeometryTolerance.hh"
81 
82 #include "G4VSensitiveDetector.hh"
84 
85 #include "G4SystemOfUnits.hh"
86 
88 // Class Implementation
90 
92  // Operators
94 
95 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
96 // {
97 // }
98 
100  // Constructors
102 
104  G4ProcessType type)
105  : G4VDiscreteProcess(processName, type)
106 {
107  if ( verboseLevel > 0) {
108  G4cout << GetProcessName() << " is created " << G4endl;
109  }
110 
112 
114  theModel = glisur;
116  theReflectivity = 1.;
117  theEfficiency = 0.;
118  theTransmittance = 0.;
119 
120  theSurfaceRoughness = 0.;
121 
122  prob_sl = 0.;
123  prob_ss = 0.;
124  prob_bs = 0.;
125 
126  PropertyPointer = NULL;
127  PropertyPointer1 = NULL;
128  PropertyPointer2 = NULL;
129 
130  Material1 = NULL;
131  Material2 = NULL;
132 
133  OpticalSurface = NULL;
134 
137 
138  iTE = iTM = 0;
139  thePhotonMomentum = 0.;
140  Rindex1 = Rindex2 = 1.;
141  cost1 = cost2 = sint1 = sint2 = 0.;
142 
143  idx = idy = 0;
144  DichroicVector = NULL;
145 
146  fInvokeSD = true;
147 }
148 
149 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
150 // {
151 // }
152 
154  // Destructors
156 
158 
160  // Methods
162 
163 // PostStepDoIt
164 // ------------
165 //
166 
169 {
171 
172  aParticleChange.Initialize(aTrack);
174 
175  // Get hyperStep from G4ParallelWorldProcess
176  // NOTE: PostSetpDoIt of this process should be
177  // invoked after G4ParallelWorldProcess!
178 
179  const G4Step* pStep = &aStep;
180 
182 
183  if (hStep) pStep = hStep;
184 
185  G4bool isOnBoundary =
187 
188  if (isOnBoundary) {
189  Material1 = pStep->GetPreStepPoint()->GetMaterial();
190  Material2 = pStep->GetPostStepPoint()->GetMaterial();
191  } else {
194  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
195  }
196 
197  G4VPhysicalVolume* thePrePV =
198  pStep->GetPreStepPoint() ->GetPhysicalVolume();
199  G4VPhysicalVolume* thePostPV =
201 
202  if ( verboseLevel > 0 ) {
203  G4cout << " Photon at Boundary! " << G4endl;
204  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
205  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
206  }
207 
208  if (aTrack.GetStepLength()<=kCarTolerance/2){
211  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
212  }
213 
214  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
215 
216  thePhotonMomentum = aParticle->GetTotalMomentum();
217  OldMomentum = aParticle->GetMomentumDirection();
218  OldPolarization = aParticle->GetPolarization();
219 
220  if ( verboseLevel > 0 ) {
221  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
222  G4cout << " Old Polarization: " << OldPolarization << G4endl;
223  }
224 
225  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
226 
227  G4bool valid;
228  // Use the new method for Exit Normal in global coordinates,
229  // which provides the normal more reliably.
230 
231  // ID of Navigator which limits step
232 
234  std::vector<G4Navigator*>::iterator iNav =
236  GetActiveNavigatorsIterator();
238  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
239 
240  if (valid) {
242  }
243  else
244  {
246  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
247  << " The Navigator reports that it returned an invalid normal"
248  << G4endl;
249  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
251  "Invalid Surface Normal - Geometry must return valid surface normal");
252  }
253 
254  if (OldMomentum * theGlobalNormal > 0.0) {
255 #ifdef G4OPTICAL_DEBUG
257  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
258  << " theGlobalNormal points in a wrong direction. "
259  << G4endl;
260  ed << " The momentum of the photon arriving at interface (oldMomentum)"
261  << " must exit the volume cross in the step. " << G4endl;
262  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
263  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
264  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
265  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
266  ed << G4endl;
267  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
268  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
269  ed,
270  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
271 #else
273 #endif
274  }
275 
276  G4MaterialPropertiesTable* aMaterialPropertiesTable;
277  G4MaterialPropertyVector* Rindex;
278 
279  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
280  if (aMaterialPropertiesTable) {
281  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
282  }
283  else {
288  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
289  }
290 
291  if (Rindex) {
292  Rindex1 = Rindex->Value(thePhotonMomentum);
293  }
294  else {
299  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
300  }
301 
302  theReflectivity = 1.;
303  theEfficiency = 0.;
304  theTransmittance = 0.;
305 
306  theSurfaceRoughness = 0.;
307 
308  theModel = glisur;
310 
312 
313  Rindex = NULL;
314  OpticalSurface = NULL;
315 
316  G4LogicalSurface* Surface = NULL;
317 
318  Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
319 
320  if (Surface == NULL){
321  G4bool enteredDaughter= (thePostPV->GetMotherLogical() ==
322  thePrePV ->GetLogicalVolume());
323  if(enteredDaughter){
324  Surface =
326  if(Surface == NULL)
327  Surface =
329  }
330  else {
331  Surface =
333  if(Surface == NULL)
334  Surface =
336  }
337  }
338 
339  if (Surface) OpticalSurface =
340  dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
341 
342  if (OpticalSurface) {
343 
344  type = OpticalSurface->GetType();
347 
348  aMaterialPropertiesTable = OpticalSurface->
349  GetMaterialPropertiesTable();
350 
351  if (aMaterialPropertiesTable) {
352 
355  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
356  if (Rindex) {
357  Rindex2 = Rindex->Value(thePhotonMomentum);
358  }
359  else {
364  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
365  }
366  }
367 
369  aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
371  aMaterialPropertiesTable->GetProperty("REALRINDEX");
373  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
374 
375  iTE = 1;
376  iTM = 1;
377 
378  if (PropertyPointer) {
379 
382 
383  } else if (PropertyPointer1 && PropertyPointer2) {
384 
386 
387  }
388 
390  aMaterialPropertiesTable->GetProperty("EFFICIENCY");
391  if (PropertyPointer) {
392  theEfficiency =
394  }
395 
397  aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
398  if (PropertyPointer) {
401  }
402 
403  if (aMaterialPropertiesTable->
404  ConstPropertyExists("SURFACEROUGHNESS"))
405  theSurfaceRoughness = aMaterialPropertiesTable->
406  GetConstProperty("SURFACEROUGHNESS");
407 
408  if ( theModel == unified ) {
410  aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
411  if (PropertyPointer) {
412  prob_sl =
414  } else {
415  prob_sl = 0.0;
416  }
417 
419  aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
420  if (PropertyPointer) {
421  prob_ss =
423  } else {
424  prob_ss = 0.0;
425  }
426 
428  aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
429  if (PropertyPointer) {
430  prob_bs =
432  } else {
433  prob_bs = 0.0;
434  }
435  }
436  }
437  else if (theFinish == polishedbackpainted ||
441  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
442  }
443  }
444 
445  if (type == dielectric_dielectric ) {
446  if (theFinish == polished || theFinish == ground ) {
447 
448  if (Material1 == Material2){
451  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
452  }
453  aMaterialPropertiesTable =
455  if (aMaterialPropertiesTable)
456  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
457  if (Rindex) {
458  Rindex2 = Rindex->Value(thePhotonMomentum);
459  }
460  else {
465  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
466  }
467  }
468  }
469 
470  if (type == dielectric_metal) {
471 
472  DielectricMetal();
473 
474  }
475  else if (type == dielectric_LUT) {
476 
477  DielectricLUT();
478 
479  }
480  else if (type == dielectric_dichroic) {
481 
483 
484  }
485  else if (type == dielectric_dielectric) {
486 
487  if ( theFinish == polishedbackpainted ||
490  }
491  else {
492  G4double rand = G4UniformRand();
493  if ( rand > theReflectivity ) {
494  if (rand > theReflectivity + theTransmittance) {
495  DoAbsorption();
496  } else {
500  }
501  }
502  else {
503  if ( theFinish == polishedfrontpainted ) {
504  DoReflection();
505  }
506  else if ( theFinish == groundfrontpainted ) {
508  DoReflection();
509  }
510  else {
512  }
513  }
514  }
515  }
516  else {
517 
518  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
519  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
520 
521  }
522 
523  NewMomentum = NewMomentum.unit();
525 
526  if ( verboseLevel > 0) {
527  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
528  G4cout << " New Polarization: " << NewPolarization << G4endl;
530  }
531 
534 
536  G4MaterialPropertyVector* groupvel =
538  G4double finalVelocity = groupvel->Value(thePhotonMomentum);
539  aParticleChange.ProposeVelocity(finalVelocity);
540  }
541 
542  if ( theStatus == Detection && fInvokeSD ) InvokeSD(pStep);
543 
544  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
545 }
546 
548 {
549  if ( theStatus == Undefined )
550  G4cout << " *** Undefined *** " << G4endl;
551  if ( theStatus == Transmission )
552  G4cout << " *** Transmission *** " << G4endl;
553  if ( theStatus == FresnelRefraction )
554  G4cout << " *** FresnelRefraction *** " << G4endl;
555  if ( theStatus == FresnelReflection )
556  G4cout << " *** FresnelReflection *** " << G4endl;
558  G4cout << " *** TotalInternalReflection *** " << G4endl;
560  G4cout << " *** LambertianReflection *** " << G4endl;
561  if ( theStatus == LobeReflection )
562  G4cout << " *** LobeReflection *** " << G4endl;
563  if ( theStatus == SpikeReflection )
564  G4cout << " *** SpikeReflection *** " << G4endl;
565  if ( theStatus == BackScattering )
566  G4cout << " *** BackScattering *** " << G4endl;
568  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
570  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
572  G4cout << " *** PolishedAirReflection *** " << G4endl;
574  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
576  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
578  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
580  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
582  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
584  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
586  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
588  G4cout << " *** EtchedAirReflection *** " << G4endl;
590  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
592  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
594  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
596  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
598  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
600  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
602  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
604  G4cout << " *** GroundAirReflection *** " << G4endl;
606  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
608  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
610  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
612  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
614  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
615  if ( theStatus == Absorption )
616  G4cout << " *** Absorption *** " << G4endl;
617  if ( theStatus == Detection )
618  G4cout << " *** Detection *** " << G4endl;
619  if ( theStatus == NotAtBoundary )
620  G4cout << " *** NotAtBoundary *** " << G4endl;
621  if ( theStatus == SameMaterial )
622  G4cout << " *** SameMaterial *** " << G4endl;
623  if ( theStatus == StepTooSmall )
624  G4cout << " *** StepTooSmall *** " << G4endl;
625  if ( theStatus == NoRINDEX )
626  G4cout << " *** NoRINDEX *** " << G4endl;
627  if ( theStatus == Dichroic )
628  G4cout << " *** Dichroic Transmission *** " << G4endl;
629 }
630 
633  const G4ThreeVector& Normal ) const
634 {
635  G4ThreeVector FacetNormal;
636 
637  if (theModel == unified || theModel == LUT) {
638 
639  /* This function code alpha to a random value taken from the
640  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
641  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
642  is a gaussian distribution with mean 0 and standard deviation
643  sigma_alpha. */
644 
645  G4double alpha;
646 
647  G4double sigma_alpha = 0.0;
648  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
649 
650  if (sigma_alpha == 0.0) return FacetNormal = Normal;
651 
652  G4double f_max = std::min(1.0,4.*sigma_alpha);
653 
654  G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
655  G4ThreeVector tmpNormal;
656 
657  do {
658  do {
659  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
660  // Loop checking, 13-Aug-2015, Peter Gumplinger
661  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
662 
663  phi = G4UniformRand()*twopi;
664 
665  SinAlpha = std::sin(alpha);
666  CosAlpha = std::cos(alpha);
667  SinPhi = std::sin(phi);
668  CosPhi = std::cos(phi);
669 
670  unit_x = SinAlpha * CosPhi;
671  unit_y = SinAlpha * SinPhi;
672  unit_z = CosAlpha;
673 
674  FacetNormal.setX(unit_x);
675  FacetNormal.setY(unit_y);
676  FacetNormal.setZ(unit_z);
677 
678  tmpNormal = Normal;
679 
680  FacetNormal.rotateUz(tmpNormal);
681  // Loop checking, 13-Aug-2015, Peter Gumplinger
682  } while (Momentum * FacetNormal >= 0.0);
683  }
684  else {
685 
686  G4double polish = 1.0;
687  if (OpticalSurface) polish = OpticalSurface->GetPolish();
688 
689  if (polish < 1.0) {
690  do {
691  G4ThreeVector smear;
692  do {
693  smear.setX(2.*G4UniformRand()-1.0);
694  smear.setY(2.*G4UniformRand()-1.0);
695  smear.setZ(2.*G4UniformRand()-1.0);
696  // Loop checking, 13-Aug-2015, Peter Gumplinger
697  } while (smear.mag()>1.0);
698  smear = (1.-polish) * smear;
699  FacetNormal = Normal + smear;
700  // Loop checking, 13-Aug-2015, Peter Gumplinger
701  } while (Momentum * FacetNormal >= 0.0);
702  FacetNormal = FacetNormal.unit();
703  }
704  else {
705  FacetNormal = Normal;
706  }
707  }
708  return FacetNormal;
709 }
710 
712 {
713  G4int n = 0;
714  G4double rand, PdotN, EdotN;
715  G4ThreeVector A_trans, A_paral;
716 
717  do {
718 
719  n++;
720 
721  rand = G4UniformRand();
722  if ( rand > theReflectivity && n == 1 ) {
723  if (rand > theReflectivity + theTransmittance) {
724  DoAbsorption();
725  } else {
729  }
730  break;
731  }
732  else {
733 
735  if ( n > 1 ) {
737  if ( !G4BooleanRand(theReflectivity) ) {
738  DoAbsorption();
739  break;
740  }
741  }
742  }
743 
744  if ( theModel == glisur || theFinish == polished ) {
745 
746  DoReflection();
747 
748  } else {
749 
750  if ( n == 1 ) ChooseReflection();
751 
752  if ( theStatus == LambertianReflection ) {
753  DoReflection();
754  }
755  else if ( theStatus == BackScattering ) {
758  }
759  else {
760 
763  } else {
766  }
767  }
768 
769  PdotN = OldMomentum * theFacetNormal;
770  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
772 
773  if (sint1 > 0.0 ) {
774  A_trans = OldMomentum.cross(theFacetNormal);
775  A_trans = A_trans.unit();
776  } else {
777  A_trans = OldPolarization;
778  }
779  A_paral = NewMomentum.cross(A_trans);
780  A_paral = A_paral.unit();
781 
782  if(iTE>0&&iTM>0) {
783  NewPolarization =
784  -OldPolarization + (2.*EdotN)*theFacetNormal;
785  } else if (iTE>0) {
786  NewPolarization = -A_trans;
787  } else if (iTM>0) {
788  NewPolarization = -A_paral;
789  }
790 
791  }
792 
793  }
794 
797 
798  }
799 
800  // Loop checking, 13-Aug-2015, Peter Gumplinger
801  } while (NewMomentum * theGlobalNormal < 0.0);
802 }
803 
805 {
806  G4int thetaIndex, phiIndex;
807  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
808  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
809 
812 
813  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
814  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
815 
816  G4double rand;
817 
818  do {
819  rand = G4UniformRand();
820  if ( rand > theReflectivity ) {
821  if (rand > theReflectivity + theTransmittance) {
822  DoAbsorption();
823  } else {
827  }
828  break;
829  }
830  else {
831  // Calculate Angle between Normal and Photon Momentum
832  G4double anglePhotonToNormal =
834  // Round it to closest integer
835  G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
836 
837  // Take random angles THETA and PHI,
838  // and see if below Probability - if not - Redo
839  do {
840  thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
841  phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
842  // Find probability with the new indeces from LUT
843  AngularDistributionValue = OpticalSurface ->
844  GetAngularDistributionValue(angleIncident,
845  thetaIndex,
846  phiIndex);
847  // Loop checking, 13-Aug-2015, Peter Gumplinger
848  } while ( !G4BooleanRand(AngularDistributionValue) );
849 
850  thetaRad = (-90 + 4*thetaIndex)*pi/180;
851  phiRad = (-90 + 5*phiIndex)*pi/180;
852  // Rotate Photon Momentum in Theta, then in Phi
854 
855  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
856  if (PerpendicularVectorTheta.mag() < kCarTolerance )
857  PerpendicularVectorTheta = NewMomentum.orthogonal();
858  NewMomentum =
859  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
860  PerpendicularVectorTheta);
861  PerpendicularVectorPhi =
862  PerpendicularVectorTheta.cross(NewMomentum);
863  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
864 
865  // Rotate Polarization too:
868  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
869  }
870  // Loop checking, 13-Aug-2015, Peter Gumplinger
871  } while (NewMomentum * theGlobalNormal <= 0.0);
872 }
873 
875 {
876  // Calculate Angle between Normal and Photon Momentum
877  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
878 
879  // Round it to closest integer
880  G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
881 
882  if (!DichroicVector) {
884  }
885 
886 
887  if (DichroicVector) {
888  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
890  DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
891 // G4cout << "wavelength: " << std::floor(wavelength/nm)
892 // << "nm" << G4endl;
893 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
894 // G4cout << "Transmittance: "
895 // << std::floor(theTransmittance/perCent) << "%" << G4endl;
896  } else {
898  ed << " G4OpBoundaryProcess/DielectricDichroic(): "
899  << " The dichroic surface has no G4Physics2DVector"
900  << G4endl;
901  G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
902  FatalException,ed,
903  "A dichroic surface must have an associated G4Physics2DVector");
904  }
905 
906  if ( !G4BooleanRand(theTransmittance) ) { // Not transmitted, so reflect
907 
908  if ( theModel == glisur || theFinish == polished ) {
909  DoReflection();
910  } else {
912  if ( theStatus == LambertianReflection ) {
913  DoReflection();
914  } else if ( theStatus == BackScattering ) {
917  } else {
918  G4double PdotN, EdotN;
919  do {
922  PdotN = OldMomentum * theFacetNormal;
923  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
924  // Loop checking, 13-Aug-2015, Peter Gumplinger
925  } while (NewMomentum * theGlobalNormal <= 0.0);
927  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
928  }
929  }
930 
931  } else {
932 
936 
937  }
938 }
939 
941 {
942  G4bool Inside = false;
943  G4bool Swap = false;
944 
945  G4bool SurfaceRoughnessCriterionPass = 1;
946  if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) {
947  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
948  G4double SurfaceRoughnessCriterion =
949  std::exp(-std::pow((4*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2));
950  SurfaceRoughnessCriterionPass =
951  G4BooleanRand(SurfaceRoughnessCriterion);
952  }
953 
954  leap:
955 
956  G4bool Through = false;
957  G4bool Done = false;
958 
959  G4double PdotN, EdotN;
960 
961  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
962  G4double E1_perp, E1_parl;
963  G4double s1, s2, E2_perp, E2_parl, E2_total, TransCoeff;
964  G4double E2_abs, C_parl, C_perp;
965  G4double alpha;
966 
967  do {
968 
969  if (Through) {
970  Swap = !Swap;
971  Through = false;
975  }
976 
977  if ( theFinish == polished ) {
979  }
980  else {
983  }
984 
985  PdotN = OldMomentum * theFacetNormal;
987 
988  cost1 = - PdotN;
989  if (std::abs(cost1) < 1.0-kCarTolerance){
990  sint1 = std::sqrt(1.-cost1*cost1);
991  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
992  }
993  else {
994  sint1 = 0.0;
995  sint2 = 0.0;
996  }
997 
998  if (sint2 >= 1.0) {
999 
1000  // Simulate total internal reflection
1001 
1002  if (Swap) Swap = !Swap;
1003 
1005 
1006  if ( !SurfaceRoughnessCriterionPass ) theStatus =
1008 
1009  if ( theModel == unified && theFinish != polished )
1010  ChooseReflection();
1011 
1012  if ( theStatus == LambertianReflection ) {
1013  DoReflection();
1014  }
1015  else if ( theStatus == BackScattering ) {
1018  }
1019  else {
1020 
1021  PdotN = OldMomentum * theFacetNormal;
1022  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1023  EdotN = OldPolarization * theFacetNormal;
1024  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1025 
1026  }
1027  }
1028  else if (sint2 < 1.0) {
1029 
1030  // Calculate amplitude for transmission (Q = P x N)
1031 
1032  if (cost1 > 0.0) {
1033  cost2 = std::sqrt(1.-sint2*sint2);
1034  }
1035  else {
1036  cost2 = -std::sqrt(1.-sint2*sint2);
1037  }
1038 
1039  if (sint1 > 0.0) {
1040  A_trans = OldMomentum.cross(theFacetNormal);
1041  A_trans = A_trans.unit();
1042  E1_perp = OldPolarization * A_trans;
1043  E1pp = E1_perp * A_trans;
1044  E1pl = OldPolarization - E1pp;
1045  E1_parl = E1pl.mag();
1046  }
1047  else {
1048  A_trans = OldPolarization;
1049  // Here we Follow Jackson's conventions and we set the
1050  // parallel component = 1 in case of a ray perpendicular
1051  // to the surface
1052  E1_perp = 0.0;
1053  E1_parl = 1.0;
1054  }
1055 
1056  s1 = Rindex1*cost1;
1057  E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1058  E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1059  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1060  s2 = Rindex2*cost2*E2_total;
1061 
1062  if (theTransmittance > 0) TransCoeff = theTransmittance;
1063  else if (cost1 != 0.0) TransCoeff = s2/s1;
1064  else TransCoeff = 0.0;
1065 
1066  if ( !G4BooleanRand(TransCoeff) ) {
1067 
1068  // Simulate reflection
1069 
1070  if (Swap) Swap = !Swap;
1071 
1073 
1074  if ( !SurfaceRoughnessCriterionPass ) theStatus =
1076 
1077  if ( theModel == unified && theFinish != polished )
1078  ChooseReflection();
1079 
1080  if ( theStatus == LambertianReflection ) {
1081  DoReflection();
1082  }
1083  else if ( theStatus == BackScattering ) {
1086  }
1087  else {
1088 
1089  PdotN = OldMomentum * theFacetNormal;
1090  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1091 
1092  if (sint1 > 0.0) { // incident ray oblique
1093 
1094  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
1095  E2_perp = E2_perp - E1_perp;
1096  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1097  A_paral = NewMomentum.cross(A_trans);
1098  A_paral = A_paral.unit();
1099  E2_abs = std::sqrt(E2_total);
1100  C_parl = E2_parl/E2_abs;
1101  C_perp = E2_perp/E2_abs;
1102 
1103  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1104 
1105  }
1106 
1107  else { // incident ray perpendicular
1108 
1109  if (Rindex2 > Rindex1) {
1111  }
1112  else {
1114  }
1115 
1116  }
1117  }
1118  }
1119  else { // photon gets transmitted
1120 
1121  // Simulate transmission/refraction
1122 
1123  Inside = !Inside;
1124  Through = true;
1126 
1127  if (sint1 > 0.0) { // incident ray oblique
1128 
1129  alpha = cost1 - cost2*(Rindex2/Rindex1);
1131  NewMomentum = NewMomentum.unit();
1132 // PdotN = -cost2;
1133  A_paral = NewMomentum.cross(A_trans);
1134  A_paral = A_paral.unit();
1135  E2_abs = std::sqrt(E2_total);
1136  C_parl = E2_parl/E2_abs;
1137  C_perp = E2_perp/E2_abs;
1138 
1139  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1140 
1141  }
1142  else { // incident ray perpendicular
1143 
1146 
1147  }
1148  }
1149  }
1150 
1151  OldMomentum = NewMomentum.unit();
1153 
1154  if (theStatus == FresnelRefraction) {
1155  Done = (NewMomentum * theGlobalNormal <= 0.0);
1156  }
1157  else {
1158  Done = (NewMomentum * theGlobalNormal >= -kCarTolerance);
1159  }
1160 
1161  // Loop checking, 13-Aug-2015, Peter Gumplinger
1162  } while (!Done);
1163 
1164  if (Inside && !Swap) {
1165  if( theFinish == polishedbackpainted ||
1167 
1168  G4double rand = G4UniformRand();
1169  if ( rand > theReflectivity ) {
1170  if (rand > theReflectivity + theTransmittance) {
1171  DoAbsorption();
1172  } else {
1176  }
1177  }
1178  else {
1179  if (theStatus != FresnelRefraction ) {
1181  }
1182  else {
1183  Swap = !Swap;
1186  }
1187  if ( theFinish == groundbackpainted )
1189 
1190  DoReflection();
1191 
1194 
1195  goto leap;
1196  }
1197  }
1198  }
1199 }
1200 
1201 // GetMeanFreePath
1202 // ---------------
1203 //
1205  G4double ,
1207 {
1208  *condition = Forced;
1209 
1210  return DBL_MAX;
1211 }
1212 
1214 {
1216  G4double magP= OldMomentum.mag();
1217  G4double magN= theFacetNormal.mag();
1218  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1219 
1220  return incidentangle;
1221 }
1222 
1224  G4double E1_parl,
1225  G4double incidentangle,
1226  G4double RealRindex,
1227  G4double ImaginaryRindex)
1228 {
1229  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1230  G4complex N1(Rindex1, 0), N2(RealRindex, ImaginaryRindex);
1231  G4complex CosPhi;
1232 
1233  G4complex u(1,0); //unit number 1
1234 
1235  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1236  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1237  G4complex denominatorTE, denominatorTM;
1238  G4complex rTM, rTE;
1239 
1240  G4MaterialPropertiesTable* aMaterialPropertiesTable =
1242  G4MaterialPropertyVector* aPropertyPointerR =
1243  aMaterialPropertiesTable->GetProperty("REALRINDEX");
1244  G4MaterialPropertyVector* aPropertyPointerI =
1245  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
1246  if (aPropertyPointerR && aPropertyPointerI) {
1247  G4double RRindex = aPropertyPointerR->Value(thePhotonMomentum);
1248  G4double IRindex = aPropertyPointerI->Value(thePhotonMomentum);
1249  N1 = G4complex(RRindex,IRindex);
1250  }
1251 
1252  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1253  // Optics" written by Fowles
1254 
1255  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))*(N1*N1)/(N2*N2)));
1256 
1257  numeratorTE = N1*std::cos(incidentangle) - N2*CosPhi;
1258  denominatorTE = N1*std::cos(incidentangle) + N2*CosPhi;
1259  rTE = numeratorTE/denominatorTE;
1260 
1261  numeratorTM = N2*std::cos(incidentangle) - N1*CosPhi;
1262  denominatorTM = N2*std::cos(incidentangle) + N1*CosPhi;
1263  rTM = numeratorTM/denominatorTM;
1264 
1265  // This is my calculaton for reflectivity on a metalic surface
1266  // depending on the fraction of TE and TM polarization
1267  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1268  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1269 
1270  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1271  / (E1_perp*E1_perp + E1_parl*E1_parl);
1272  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1273  / (E1_perp*E1_perp + E1_parl*E1_parl);
1274  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1275 
1276  do {
1277  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1278  {iTE = -1;}else{iTE = 1;}
1279  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1280  {iTM = -1;}else{iTM = 1;}
1281  // Loop checking, 13-Aug-2015, Peter Gumplinger
1282  } while(iTE<0&&iTM<0);
1283 
1284  return real(Reflectivity);
1285 
1286 }
1287 
1289 {
1290  G4double RealRindex =
1292  G4double ImaginaryRindex =
1294 
1295  // calculate FacetNormal
1296  if ( theFinish == ground ) {
1297  theFacetNormal =
1299  } else {
1301  }
1302 
1304  cost1 = -PdotN;
1305 
1306  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1307  sint1 = std::sqrt(1. - cost1*cost1);
1308  } else {
1309  sint1 = 0.0;
1310  }
1311 
1312  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1313  G4double E1_perp, E1_parl;
1314 
1315  if (sint1 > 0.0 ) {
1316  A_trans = OldMomentum.cross(theFacetNormal);
1317  A_trans = A_trans.unit();
1318  E1_perp = OldPolarization * A_trans;
1319  E1pp = E1_perp * A_trans;
1320  E1pl = OldPolarization - E1pp;
1321  E1_parl = E1pl.mag();
1322  }
1323  else {
1324  A_trans = OldPolarization;
1325  // Here we Follow Jackson's conventions and we set the
1326  // parallel component = 1 in case of a ray perpendicular
1327  // to the surface
1328  E1_perp = 0.0;
1329  E1_parl = 1.0;
1330  }
1331 
1332  //calculate incident angle
1333  G4double incidentangle = GetIncidentAngle();
1334 
1335  //calculate the reflectivity depending on incident angle,
1336  //polarization and complex refractive
1337 
1338  theReflectivity =
1339  GetReflectivity(E1_perp, E1_parl, incidentangle,
1340  RealRindex, ImaginaryRindex);
1341 }
1342 
1344 {
1345  G4Step aStep = *pStep;
1346 
1348 
1350  if (sd) return sd->Hit(&aStep);
1351  else return false;
1352 }
G4double condition(const G4ErrorSymMatrix &m)
ThreeVector shoot(const G4int Ap, const G4int Af)
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4OpticalSurfaceModel GetModel() const
G4int verboseLevel
Definition: G4VProcess.hh:368
G4double GetReflectivity(G4double E1_perp, G4double E1_parl, G4double incidentangle, G4double RealRindex, G4double ImaginaryRindex)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
G4Physics2DVector * DichroicVector
G4double GetVelocity() const
const G4DynamicParticle * GetDynamicParticle() const
static constexpr double perCent
Definition: G4SIunits.hh:332
G4StepStatus GetStepStatus() const
G4Material * GetMaterial() const
G4OpticalSurface * OpticalSurface
G4double GetSurfaceTolerance() const
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
const G4SurfaceType & GetType() const
G4OpBoundaryProcessStatus
G4OpBoundaryProcess(const G4String &processName="OpBoundary", G4ProcessType type=fOptical)
G4SurfaceType
G4MaterialPropertyVector * PropertyPointer1
int G4int
Definition: G4Types.hh:78
G4OpticalSurfaceFinish theFinish
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
static G4LogicalBorderSurface * GetSurface(const G4VPhysicalVolume *vol1, const G4VPhysicalVolume *vol2)
G4double GetTotalMomentum() const
G4StepPoint * GetPreStepPoint() const
G4OpBoundaryProcessStatus theStatus
std::complex< G4double > G4complex
Definition: G4Types.hh:81
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double Value(G4double x, G4double y, size_t &lastidx, size_t &lastidy) const
G4double GetPolish() const
G4VPhysicalVolume * GetPhysicalVolume() const
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
G4int GetPhiIndexMax(void) const
G4bool G4BooleanRand(const G4double prob) const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4bool Hit(G4Step *aStep)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *condition)
void BoundaryProcessVerbose(void) const
Definition: G4Step.hh:76
const G4int n
G4LogicalVolume * GetMotherLogical() const
G4double Value(G4double theEnergy, size_t &lastidx) const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4OpticalSurfaceFinish GetFinish() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
G4OpticalSurfaceModel theModel
G4bool InvokeSD(const G4Step *step)
virtual void Initialize(const G4Track &)
G4LogicalVolume * GetLogicalVolume() const
G4MaterialPropertyVector * PropertyPointer
G4Physics2DVector * GetDichroicVector()
void G4SwapPtr(T *&a, T *&b)
Definition: templates.hh:121
static constexpr double nm
Definition: G4SIunits.hh:112
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
G4StepPoint * GetPostStepPoint() const
void AddTotalEnergyDeposit(G4double value)
const G4ThreeVector & GetPolarization() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
G4SurfaceProperty * GetSurfaceProperty() const
static constexpr double pi
Definition: G4SIunits.hh:75
G4VSensitiveDetector * GetSensitiveDetector() const
static constexpr double halfpi
Definition: G4SIunits.hh:77
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4MaterialPropertyVector * PropertyPointer2
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4double GetSigmaAlpha() const
G4ForceCondition
static const G4double alpha
void ProposeVelocity(G4double finalVelocity)
#define DBL_MAX
Definition: templates.hh:83
G4int GetThetaIndexMax(void) const
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4MaterialPropertyVector * GetProperty(const char *key)
G4double GetStepLength() const
static G4GeometryTolerance * GetInstance()
static const G4Step * GetHyperStep()
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GetFacetNormal(const G4ThreeVector &Momentum, const G4ThreeVector &Normal) const
G4ProcessType
static G4LogicalSkinSurface * GetSurface(const G4LogicalVolume *vol)