Geant4  10.01.p03
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 = cost1 = cost2 = sint1 = sint2 = 0.;
141 
142  idx = idy = 0;
143  DichroicVector = NULL;
144 }
145 
146 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
147 // {
148 // }
149 
151  // Destructors
153 
155 
157  // Methods
159 
160 // PostStepDoIt
161 // ------------
162 //
163 
166 {
168 
169  aParticleChange.Initialize(aTrack);
171 
172  // Get hyperStep from G4ParallelWorldProcess
173  // NOTE: PostSetpDoIt of this process should be
174  // invoked after G4ParallelWorldProcess!
175 
176  const G4Step* pStep = &aStep;
177 
179 
180  if (hStep) pStep = hStep;
181 
182  G4bool isOnBoundary =
184 
185  if (isOnBoundary) {
186  Material1 = pStep->GetPreStepPoint()->GetMaterial();
187  Material2 = pStep->GetPostStepPoint()->GetMaterial();
188  } else {
191  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
192  }
193 
194  G4VPhysicalVolume* thePrePV =
195  pStep->GetPreStepPoint() ->GetPhysicalVolume();
196  G4VPhysicalVolume* thePostPV =
198 
199  if ( verboseLevel > 0 ) {
200  G4cout << " Photon at Boundary! " << G4endl;
201  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
202  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
203  }
204 
205  if (aTrack.GetStepLength()<=kCarTolerance/2){
208  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
209  }
210 
211  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
212 
213  thePhotonMomentum = aParticle->GetTotalMomentum();
214  OldMomentum = aParticle->GetMomentumDirection();
215  OldPolarization = aParticle->GetPolarization();
216 
217  if ( verboseLevel > 0 ) {
218  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
219  G4cout << " Old Polarization: " << OldPolarization << G4endl;
220  }
221 
222  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
223 
224  G4bool valid;
225  // Use the new method for Exit Normal in global coordinates,
226  // which provides the normal more reliably.
227 
228  // ID of Navigator which limits step
229 
231  std::vector<G4Navigator*>::iterator iNav =
233  GetActiveNavigatorsIterator();
235  (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
236 
237  if (valid) {
239  }
240  else
241  {
243  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
244  << " The Navigator reports that it returned an invalid normal"
245  << G4endl;
246  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
248  "Invalid Surface Normal - Geometry must return valid surface normal");
249  }
250 
251  if (OldMomentum * theGlobalNormal > 0.0) {
252 #ifdef G4OPTICAL_DEBUG
254  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
255  << " theGlobalNormal points in a wrong direction. "
256  << G4endl;
257  ed << " The momentum of the photon arriving at interface (oldMomentum)"
258  << " must exit the volume cross in the step. " << G4endl;
259  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
260  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
261  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
262  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
263  ed << G4endl;
264  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
265  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
266  ed,
267  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
268 #else
270 #endif
271  }
272 
273  G4MaterialPropertiesTable* aMaterialPropertiesTable;
274  G4MaterialPropertyVector* Rindex;
275 
276  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
277  if (aMaterialPropertiesTable) {
278  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
279  }
280  else {
285  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
286  }
287 
288  if (Rindex) {
289  Rindex1 = Rindex->Value(thePhotonMomentum);
290  }
291  else {
296  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
297  }
298 
299  theReflectivity = 1.;
300  theEfficiency = 0.;
301  theTransmittance = 0.;
302 
303  theSurfaceRoughness = 0.;
304 
305  theModel = glisur;
307 
309 
310  Rindex = NULL;
311  OpticalSurface = NULL;
312 
313  G4LogicalSurface* Surface = NULL;
314 
315  Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
316 
317  if (Surface == NULL){
318  G4bool enteredDaughter= (thePostPV->GetMotherLogical() ==
319  thePrePV ->GetLogicalVolume());
320  if(enteredDaughter){
321  Surface =
323  if(Surface == NULL)
324  Surface =
326  }
327  else {
328  Surface =
330  if(Surface == NULL)
331  Surface =
333  }
334  }
335 
336  if (Surface) OpticalSurface =
337  dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
338 
339  if (OpticalSurface) {
340 
341  type = OpticalSurface->GetType();
344 
345  aMaterialPropertiesTable = OpticalSurface->
346  GetMaterialPropertiesTable();
347 
348  if (aMaterialPropertiesTable) {
349 
352  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
353  if (Rindex) {
354  Rindex2 = Rindex->Value(thePhotonMomentum);
355  }
356  else {
361  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
362  }
363  }
364 
366  aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
368  aMaterialPropertiesTable->GetProperty("REALRINDEX");
370  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
371 
372  iTE = 1;
373  iTM = 1;
374 
375  if (PropertyPointer) {
376 
379 
380  } else if (PropertyPointer1 && PropertyPointer2) {
381 
383 
384  }
385 
387  aMaterialPropertiesTable->GetProperty("EFFICIENCY");
388  if (PropertyPointer) {
389  theEfficiency =
391  }
392 
394  aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
395  if (PropertyPointer) {
398  }
399 
400  if (aMaterialPropertiesTable->
401  ConstPropertyExists("SURFACEROUGHNESS"))
402  theSurfaceRoughness = aMaterialPropertiesTable->
403  GetConstProperty("SURFACEROUGHNESS");
404 
405  if ( theModel == unified ) {
407  aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
408  if (PropertyPointer) {
409  prob_sl =
411  } else {
412  prob_sl = 0.0;
413  }
414 
416  aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
417  if (PropertyPointer) {
418  prob_ss =
420  } else {
421  prob_ss = 0.0;
422  }
423 
425  aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
426  if (PropertyPointer) {
427  prob_bs =
429  } else {
430  prob_bs = 0.0;
431  }
432  }
433  }
434  else if (theFinish == polishedbackpainted ||
438  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
439  }
440  }
441 
442  if (type == dielectric_dielectric ) {
443  if (theFinish == polished || theFinish == ground ) {
444 
445  if (Material1 == Material2){
448  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
449  }
450  aMaterialPropertiesTable =
452  if (aMaterialPropertiesTable)
453  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
454  if (Rindex) {
455  Rindex2 = Rindex->Value(thePhotonMomentum);
456  }
457  else {
462  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
463  }
464  }
465  }
466 
467  if (type == dielectric_metal) {
468 
469  DielectricMetal();
470 
471  }
472  else if (type == dielectric_LUT) {
473 
474  DielectricLUT();
475 
476  }
477  else if (type == dielectric_dichroic) {
478 
480 
481  }
482  else if (type == dielectric_dielectric) {
483 
484  if ( theFinish == polishedbackpainted ||
487  }
488  else {
489  G4double rand = G4UniformRand();
490  if ( rand > theReflectivity ) {
491  if (rand > theReflectivity + theTransmittance) {
492  DoAbsorption();
493  } else {
497  }
498  }
499  else {
500  if ( theFinish == polishedfrontpainted ) {
501  DoReflection();
502  }
503  else if ( theFinish == groundfrontpainted ) {
505  DoReflection();
506  }
507  else {
509  }
510  }
511  }
512  }
513  else {
514 
515  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
516  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
517 
518  }
519 
520  NewMomentum = NewMomentum.unit();
522 
523  if ( verboseLevel > 0) {
524  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
525  G4cout << " New Polarization: " << NewPolarization << G4endl;
527  }
528 
531 
533  G4MaterialPropertyVector* groupvel =
535  G4double finalVelocity = groupvel->Value(thePhotonMomentum);
536  aParticleChange.ProposeVelocity(finalVelocity);
537  }
538 
539  if ( theStatus == Detection ) InvokeSD(pStep);
540 
541  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
542 }
543 
545 {
546  if ( theStatus == Undefined )
547  G4cout << " *** Undefined *** " << G4endl;
548  if ( theStatus == Transmission )
549  G4cout << " *** Transmission *** " << G4endl;
550  if ( theStatus == FresnelRefraction )
551  G4cout << " *** FresnelRefraction *** " << G4endl;
552  if ( theStatus == FresnelReflection )
553  G4cout << " *** FresnelReflection *** " << G4endl;
555  G4cout << " *** TotalInternalReflection *** " << G4endl;
557  G4cout << " *** LambertianReflection *** " << G4endl;
558  if ( theStatus == LobeReflection )
559  G4cout << " *** LobeReflection *** " << G4endl;
560  if ( theStatus == SpikeReflection )
561  G4cout << " *** SpikeReflection *** " << G4endl;
562  if ( theStatus == BackScattering )
563  G4cout << " *** BackScattering *** " << G4endl;
565  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
567  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
569  G4cout << " *** PolishedAirReflection *** " << G4endl;
571  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
573  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
575  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
577  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
579  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
581  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
583  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
585  G4cout << " *** EtchedAirReflection *** " << G4endl;
587  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
589  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
591  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
593  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
595  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
597  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
599  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
601  G4cout << " *** GroundAirReflection *** " << G4endl;
603  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
605  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
607  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
609  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
611  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
612  if ( theStatus == Absorption )
613  G4cout << " *** Absorption *** " << G4endl;
614  if ( theStatus == Detection )
615  G4cout << " *** Detection *** " << G4endl;
616  if ( theStatus == NotAtBoundary )
617  G4cout << " *** NotAtBoundary *** " << G4endl;
618  if ( theStatus == SameMaterial )
619  G4cout << " *** SameMaterial *** " << G4endl;
620  if ( theStatus == StepTooSmall )
621  G4cout << " *** StepTooSmall *** " << G4endl;
622  if ( theStatus == NoRINDEX )
623  G4cout << " *** NoRINDEX *** " << G4endl;
624  if ( theStatus == Dichroic )
625  G4cout << " *** Dichroic Transmission *** " << G4endl;
626 }
627 
630  const G4ThreeVector& Normal ) const
631 {
632  G4ThreeVector FacetNormal;
633 
634  if (theModel == unified || theModel == LUT) {
635 
636  /* This function code alpha to a random value taken from the
637  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
638  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
639  is a gaussian distribution with mean 0 and standard deviation
640  sigma_alpha. */
641 
642  G4double alpha;
643 
644  G4double sigma_alpha = 0.0;
645  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
646 
647  if (sigma_alpha == 0.0) return FacetNormal = Normal;
648 
649  G4double f_max = std::min(1.0,4.*sigma_alpha);
650 
651  do {
652  do {
653  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
654  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
655 
656  G4double phi = G4UniformRand()*twopi;
657 
658  G4double SinAlpha = std::sin(alpha);
659  G4double CosAlpha = std::cos(alpha);
660  G4double SinPhi = std::sin(phi);
661  G4double CosPhi = std::cos(phi);
662 
663  G4double unit_x = SinAlpha * CosPhi;
664  G4double unit_y = SinAlpha * SinPhi;
665  G4double unit_z = CosAlpha;
666 
667  FacetNormal.setX(unit_x);
668  FacetNormal.setY(unit_y);
669  FacetNormal.setZ(unit_z);
670 
671  G4ThreeVector tmpNormal = Normal;
672 
673  FacetNormal.rotateUz(tmpNormal);
674  } while (Momentum * FacetNormal >= 0.0);
675  }
676  else {
677 
678  G4double polish = 1.0;
679  if (OpticalSurface) polish = OpticalSurface->GetPolish();
680 
681  if (polish < 1.0) {
682  do {
683  G4ThreeVector smear;
684  do {
685  smear.setX(2.*G4UniformRand()-1.0);
686  smear.setY(2.*G4UniformRand()-1.0);
687  smear.setZ(2.*G4UniformRand()-1.0);
688  } while (smear.mag()>1.0);
689  smear = (1.-polish) * smear;
690  FacetNormal = Normal + smear;
691  } while (Momentum * FacetNormal >= 0.0);
692  FacetNormal = FacetNormal.unit();
693  }
694  else {
695  FacetNormal = Normal;
696  }
697  }
698  return FacetNormal;
699 }
700 
702 {
703  G4int n = 0;
704 
705  do {
706 
707  n++;
708 
709  G4double rand = G4UniformRand();
710  if ( rand > theReflectivity && n == 1 ) {
711  if (rand > theReflectivity + theTransmittance) {
712  DoAbsorption();
713  } else {
717  }
718  break;
719  }
720  else {
721 
723  if ( n > 1 ) {
725  if ( !G4BooleanRand(theReflectivity) ) {
726  DoAbsorption();
727  break;
728  }
729  }
730  }
731 
732  if ( theModel == glisur || theFinish == polished ) {
733 
734  DoReflection();
735 
736  } else {
737 
738  if ( n == 1 ) ChooseReflection();
739 
740  if ( theStatus == LambertianReflection ) {
741  DoReflection();
742  }
743  else if ( theStatus == BackScattering ) {
746  }
747  else {
748 
751  } else {
754  }
755  }
756 
758  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
760 
761  G4ThreeVector A_trans, A_paral;
762 
763  if (sint1 > 0.0 ) {
764  A_trans = OldMomentum.cross(theFacetNormal);
765  A_trans = A_trans.unit();
766  } else {
767  A_trans = OldPolarization;
768  }
769  A_paral = NewMomentum.cross(A_trans);
770  A_paral = A_paral.unit();
771 
772  if(iTE>0&&iTM>0) {
773  NewPolarization =
774  -OldPolarization + (2.*EdotN)*theFacetNormal;
775  } else if (iTE>0) {
776  NewPolarization = -A_trans;
777  } else if (iTM>0) {
778  NewPolarization = -A_paral;
779  }
780 
781  }
782 
783  }
784 
787 
788  }
789 
790  } while (NewMomentum * theGlobalNormal < 0.0);
791 }
792 
794 {
795  G4int thetaIndex, phiIndex;
796  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
797  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
798 
801 
802  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
803  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
804 
805  do {
806  G4double rand = G4UniformRand();
807  if ( rand > theReflectivity ) {
808  if (rand > theReflectivity + theTransmittance) {
809  DoAbsorption();
810  } else {
814  }
815  break;
816  }
817  else {
818  // Calculate Angle between Normal and Photon Momentum
819  G4double anglePhotonToNormal =
821  // Round it to closest integer
822  G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
823 
824  // Take random angles THETA and PHI,
825  // and see if below Probability - if not - Redo
826  do {
827  thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
828  phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
829  // Find probability with the new indeces from LUT
830  AngularDistributionValue = OpticalSurface ->
831  GetAngularDistributionValue(angleIncident,
832  thetaIndex,
833  phiIndex);
834  } while ( !G4BooleanRand(AngularDistributionValue) );
835 
836  thetaRad = (-90 + 4*thetaIndex)*pi/180;
837  phiRad = (-90 + 5*phiIndex)*pi/180;
838  // Rotate Photon Momentum in Theta, then in Phi
840  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
841  if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
842  PerpendicularVectorPhi =
843  PerpendicularVectorTheta.cross(NewMomentum);
844  }
845  else {
846  PerpendicularVectorTheta = NewMomentum.orthogonal();
847  PerpendicularVectorPhi =
848  PerpendicularVectorTheta.cross(NewMomentum);
849  }
850  NewMomentum =
851  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
852  PerpendicularVectorTheta);
853  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
854  // Rotate Polarization too:
857  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
858  }
859  } while (NewMomentum * theGlobalNormal <= 0.0);
860 }
861 
863 {
864  // Calculate Angle between Normal and Photon Momentum
865  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
866 
867  // Round it to closest integer
868  G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
869 
870  if (!DichroicVector) {
872  }
873 
874 
875  if (DichroicVector) {
876  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
878  DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
879 // G4cout << "wavelength: " << std::floor(wavelength/nm)
880 // << "nm" << G4endl;
881 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
882 // G4cout << "Transmittance: "
883 // << std::floor(theTransmittance/perCent) << "%" << G4endl;
884  } else {
886  ed << " G4OpBoundaryProcess/DielectricDichroic(): "
887  << " The dichroic surface has no G4Physics2DVector"
888  << G4endl;
889  G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
890  FatalException,ed,
891  "A dichroic surface must have an associated G4Physics2DVector");
892  }
893 
894  if ( !G4BooleanRand(theTransmittance) ) { // Not transmitted, so reflect
895 
896  if ( theModel == glisur || theFinish == polished ) {
897  DoReflection();
898  } else {
900  if ( theStatus == LambertianReflection ) {
901  DoReflection();
902  } else if ( theStatus == BackScattering ) {
905  } else {
906  do {
910  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
911  } while (NewMomentum * theGlobalNormal <= 0.0);
913  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
914  }
915  }
916 
917  } else {
918 
922 
923  }
924 }
925 
927 {
928  G4bool Inside = false;
929  G4bool Swap = false;
930 
931  G4bool SurfaceRoughnessCriterionPass = 1;
932  if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) {
933  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
934  G4double SurfaceRoughnessCriterion =
935  std::exp(-std::pow((4*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2));
936  SurfaceRoughnessCriterionPass =
937  G4BooleanRand(SurfaceRoughnessCriterion);
938  }
939 
940  leap:
941 
942  G4bool Through = false;
943  G4bool Done = false;
944 
945  do {
946 
947  if (Through) {
948  Swap = !Swap;
949  Through = false;
953  }
954 
955  if ( theFinish == polished ) {
957  }
958  else {
961  }
962 
965 
966  cost1 = - PdotN;
967  if (std::abs(cost1) < 1.0-kCarTolerance){
968  sint1 = std::sqrt(1.-cost1*cost1);
969  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
970  }
971  else {
972  sint1 = 0.0;
973  sint2 = 0.0;
974  }
975 
976  if (sint2 >= 1.0) {
977 
978  // Simulate total internal reflection
979 
980  if (Swap) Swap = !Swap;
981 
983 
984  if ( !SurfaceRoughnessCriterionPass ) theStatus =
986 
987  if ( theModel == unified && theFinish != polished )
989 
990  if ( theStatus == LambertianReflection ) {
991  DoReflection();
992  }
993  else if ( theStatus == BackScattering ) {
996  }
997  else {
998 
999  PdotN = OldMomentum * theFacetNormal;
1000  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1001  EdotN = OldPolarization * theFacetNormal;
1002  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1003 
1004  }
1005  }
1006  else if (sint2 < 1.0) {
1007 
1008  // Calculate amplitude for transmission (Q = P x N)
1009 
1010  if (cost1 > 0.0) {
1011  cost2 = std::sqrt(1.-sint2*sint2);
1012  }
1013  else {
1014  cost2 = -std::sqrt(1.-sint2*sint2);
1015  }
1016 
1017  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1018  G4double E1_perp, E1_parl;
1019 
1020  if (sint1 > 0.0) {
1021  A_trans = OldMomentum.cross(theFacetNormal);
1022  A_trans = A_trans.unit();
1023  E1_perp = OldPolarization * A_trans;
1024  E1pp = E1_perp * A_trans;
1025  E1pl = OldPolarization - E1pp;
1026  E1_parl = E1pl.mag();
1027  }
1028  else {
1029  A_trans = OldPolarization;
1030  // Here we Follow Jackson's conventions and we set the
1031  // parallel component = 1 in case of a ray perpendicular
1032  // to the surface
1033  E1_perp = 0.0;
1034  E1_parl = 1.0;
1035  }
1036 
1037  G4double s1 = Rindex1*cost1;
1038  G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1039  G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1040  G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1041  G4double s2 = Rindex2*cost2*E2_total;
1042 
1043  G4double TransCoeff;
1044 
1045  if (theTransmittance > 0) TransCoeff = theTransmittance;
1046  else if (cost1 != 0.0) TransCoeff = s2/s1;
1047  else TransCoeff = 0.0;
1048 
1049  G4double E2_abs, C_parl, C_perp;
1050 
1051  if ( !G4BooleanRand(TransCoeff) ) {
1052 
1053  // Simulate reflection
1054 
1055  if (Swap) Swap = !Swap;
1056 
1058 
1059  if ( !SurfaceRoughnessCriterionPass ) theStatus =
1061 
1062  if ( theModel == unified && theFinish != polished )
1063  ChooseReflection();
1064 
1065  if ( theStatus == LambertianReflection ) {
1066  DoReflection();
1067  }
1068  else if ( theStatus == BackScattering ) {
1071  }
1072  else {
1073 
1074  PdotN = OldMomentum * theFacetNormal;
1075  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1076 
1077  if (sint1 > 0.0) { // incident ray oblique
1078 
1079  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
1080  E2_perp = E2_perp - E1_perp;
1081  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1082  A_paral = NewMomentum.cross(A_trans);
1083  A_paral = A_paral.unit();
1084  E2_abs = std::sqrt(E2_total);
1085  C_parl = E2_parl/E2_abs;
1086  C_perp = E2_perp/E2_abs;
1087 
1088  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1089 
1090  }
1091 
1092  else { // incident ray perpendicular
1093 
1094  if (Rindex2 > Rindex1) {
1096  }
1097  else {
1099  }
1100 
1101  }
1102  }
1103  }
1104  else { // photon gets transmitted
1105 
1106  // Simulate transmission/refraction
1107 
1108  Inside = !Inside;
1109  Through = true;
1111 
1112  if (sint1 > 0.0) { // incident ray oblique
1113 
1114  G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
1116  NewMomentum = NewMomentum.unit();
1117  PdotN = -cost2;
1118  A_paral = NewMomentum.cross(A_trans);
1119  A_paral = A_paral.unit();
1120  E2_abs = std::sqrt(E2_total);
1121  C_parl = E2_parl/E2_abs;
1122  C_perp = E2_perp/E2_abs;
1123 
1124  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1125 
1126  }
1127  else { // incident ray perpendicular
1128 
1131 
1132  }
1133  }
1134  }
1135 
1136  OldMomentum = NewMomentum.unit();
1138 
1139  if (theStatus == FresnelRefraction) {
1140  Done = (NewMomentum * theGlobalNormal <= 0.0);
1141  }
1142  else {
1143  Done = (NewMomentum * theGlobalNormal >= 0.0);
1144  }
1145 
1146  } while (!Done);
1147 
1148  if (Inside && !Swap) {
1149  if( theFinish == polishedbackpainted ||
1151 
1152  G4double rand = G4UniformRand();
1153  if ( rand > theReflectivity ) {
1154  if (rand > theReflectivity + theTransmittance) {
1155  DoAbsorption();
1156  } else {
1160  }
1161  }
1162  else {
1163  if (theStatus != FresnelRefraction ) {
1165  }
1166  else {
1167  Swap = !Swap;
1170  }
1171  if ( theFinish == groundbackpainted )
1173 
1174  DoReflection();
1175 
1178 
1179  goto leap;
1180  }
1181  }
1182  }
1183 }
1184 
1185 // GetMeanFreePath
1186 // ---------------
1187 //
1189  G4double ,
1191 {
1192  *condition = Forced;
1193 
1194  return DBL_MAX;
1195 }
1196 
1198 {
1200  G4double magP= OldMomentum.mag();
1201  G4double magN= theFacetNormal.mag();
1202  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1203 
1204  return incidentangle;
1205 }
1206 
1208  G4double E1_parl,
1209  G4double incidentangle,
1210  G4double RealRindex,
1211  G4double ImaginaryRindex)
1212 {
1213 
1214  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1215  G4complex N(RealRindex, ImaginaryRindex);
1216  G4complex CosPhi;
1217 
1218  G4complex u(1,0); //unit number 1
1219 
1220  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1221  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1222  G4complex denominatorTE, denominatorTM;
1223  G4complex rTM, rTE;
1224 
1225  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1226  // Optics" written by Fowles
1227 
1228  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
1229 
1230  numeratorTE = std::cos(incidentangle) - N*CosPhi;
1231  denominatorTE = std::cos(incidentangle) + N*CosPhi;
1232  rTE = numeratorTE/denominatorTE;
1233 
1234  numeratorTM = N*std::cos(incidentangle) - CosPhi;
1235  denominatorTM = N*std::cos(incidentangle) + CosPhi;
1236  rTM = numeratorTM/denominatorTM;
1237 
1238  // This is my calculaton for reflectivity on a metalic surface
1239  // depending on the fraction of TE and TM polarization
1240  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1241  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1242 
1243  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1244  / (E1_perp*E1_perp + E1_parl*E1_parl);
1245  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1246  / (E1_perp*E1_perp + E1_parl*E1_parl);
1247  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1248 
1249  do {
1250  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1251  {iTE = -1;}else{iTE = 1;}
1252  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1253  {iTM = -1;}else{iTM = 1;}
1254  } while(iTE<0&&iTM<0);
1255 
1256  return real(Reflectivity);
1257 
1258 }
1259 
1261 {
1262  G4double RealRindex =
1264  G4double ImaginaryRindex =
1266 
1267  // calculate FacetNormal
1268  if ( theFinish == ground ) {
1269  theFacetNormal =
1271  } else {
1273  }
1274 
1276  cost1 = -PdotN;
1277 
1278  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1279  sint1 = std::sqrt(1. - cost1*cost1);
1280  } else {
1281  sint1 = 0.0;
1282  }
1283 
1284  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1285  G4double E1_perp, E1_parl;
1286 
1287  if (sint1 > 0.0 ) {
1288  A_trans = OldMomentum.cross(theFacetNormal);
1289  A_trans = A_trans.unit();
1290  E1_perp = OldPolarization * A_trans;
1291  E1pp = E1_perp * A_trans;
1292  E1pl = OldPolarization - E1pp;
1293  E1_parl = E1pl.mag();
1294  }
1295  else {
1296  A_trans = OldPolarization;
1297  // Here we Follow Jackson's conventions and we set the
1298  // parallel component = 1 in case of a ray perpendicular
1299  // to the surface
1300  E1_perp = 0.0;
1301  E1_parl = 1.0;
1302  }
1303 
1304  //calculate incident angle
1305  G4double incidentangle = GetIncidentAngle();
1306 
1307  //calculate the reflectivity depending on incident angle,
1308  //polarization and complex refractive
1309 
1310  theReflectivity =
1311  GetReflectivity(E1_perp, E1_parl, incidentangle,
1312  RealRindex, ImaginaryRindex);
1313 }
1314 
1316 {
1317  G4Step aStep = *pStep;
1318 
1320 
1322  if (sd) return sd->Hit(&aStep);
1323  else return false;
1324 }
G4double condition(const G4ErrorSymMatrix &m)
ThreeVector shoot(const G4int Ap, const G4int Af)
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4MaterialPropertyVector * GetProperty(const char *key)
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
G4StepStatus GetStepStatus() const
G4Material * GetMaterial() const
G4OpticalSurface * OpticalSurface
const G4double pi
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 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:93
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)
static const double perCent
Definition: G4SIunits.hh:296
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
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
G4VSensitiveDetector * GetSensitiveDetector() const
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 &)
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)