ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4OpBoundaryProcess.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4OpBoundaryProcess.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
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 // 2017-02-24 - add capability of simulating surface reflections
69 // with Look-Up-Tables (LUT) developed in DAVIS
70 //
71 // Author: Peter Gumplinger
72 // adopted from work by Werner Keil - April 2/96
73 // mail: gum@triumf.ca
74 //
76 
77 #include "G4ios.hh"
78 #include "G4PhysicalConstants.hh"
79 #include "G4OpProcessSubType.hh"
80 
81 #include "G4OpBoundaryProcess.hh"
82 #include "G4GeometryTolerance.hh"
83 
84 #include "G4VSensitiveDetector.hh"
85 #include "G4ParallelWorldProcess.hh"
86 
87 #include "G4SystemOfUnits.hh"
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
92  G4ProcessType type)
93  : G4VDiscreteProcess(processName, type)
94 {
95  if ( verboseLevel > 0) {
96  G4cout << GetProcessName() << " is created " << G4endl;
97  }
98 
99  SetProcessSubType(fOpBoundary);
100 
101  theStatus = Undefined;
102  theModel = glisur;
103  theFinish = polished;
104  theReflectivity = 1.;
105  theEfficiency = 0.;
106  theTransmittance = 0.;
107 
108  theSurfaceRoughness = 0.;
109 
110  prob_sl = 0.;
111  prob_ss = 0.;
112  prob_bs = 0.;
113 
114  //PropertyPointer = nullptr;
115  //PropertyPointer1 = nullptr;
116  //PropertyPointer2 = nullptr;
117 
118  fRealRIndexMPV = nullptr;
119  fImagRIndexMPV = nullptr;
120 
121  Material1 = nullptr;
122  Material2 = nullptr;
123 
124  OpticalSurface = nullptr;
125 
126  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
127 
128  iTE = iTM = 0;
129  thePhotonMomentum = 0.;
130  Rindex1 = Rindex2 = 1.;
131  cost1 = cost2 = sint1 = sint2 = 0.;
132 
133  idx = idy = 0;
134  DichroicVector = nullptr;
135 
136  fInvokeSD = true;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
141 G4OpBoundaryProcess::~G4OpBoundaryProcess()
142 {}
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 
146 G4VParticleChange*
147 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
148 {
149  theStatus = Undefined;
150 
151  aParticleChange.Initialize(aTrack);
152  aParticleChange.ProposeVelocity(aTrack.GetVelocity());
153 
154  // Get hyperStep from G4ParallelWorldProcess
155  // NOTE: PostSetpDoIt of this process should be
156  // invoked after G4ParallelWorldProcess!
157 
158  const G4Step* pStep = &aStep;
159 
160  const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
161 
162  if (hStep) pStep = hStep;
163 
164  G4bool isOnBoundary = (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
165 
166  if (isOnBoundary) {
167  Material1 = pStep->GetPreStepPoint()->GetMaterial();
168  Material2 = pStep->GetPostStepPoint()->GetMaterial();
169  } else {
170  theStatus = NotAtBoundary;
171  if (verboseLevel > 0) BoundaryProcessVerbose();
172  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
173  }
174 
175  G4VPhysicalVolume* thePrePV = pStep->GetPreStepPoint() ->GetPhysicalVolume();
176  G4VPhysicalVolume* thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume();
177 
178  if (verboseLevel > 0) {
179  G4cout << " Photon at Boundary! " << G4endl;
180  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
181  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
182  }
183 
184  if (aTrack.GetStepLength() <= kCarTolerance/2.) {
185  theStatus = StepTooSmall;
186  if (verboseLevel > 0) BoundaryProcessVerbose();
187  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
188  }
189 
190  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
191 
192  thePhotonMomentum = aParticle->GetTotalMomentum();
193  OldMomentum = aParticle->GetMomentumDirection();
194  OldPolarization = aParticle->GetPolarization();
195 
196  if ( verboseLevel > 0 ) {
197  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
198  G4cout << " Old Polarization: " << OldPolarization << G4endl;
199  }
200 
201  G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
202 
203  G4bool valid;
204  // Use the new method for Exit Normal in global coordinates,
205  // which provides the normal more reliably.
206 
207  // ID of Navigator which limits step
208  G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
209  std::vector<G4Navigator*>::iterator iNav =
210  G4TransportationManager::GetTransportationManager()->GetActiveNavigatorsIterator();
211  theGlobalNormal = (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint, &valid);
212 
213  if (valid) {
214  theGlobalNormal = -theGlobalNormal;
215  }
216  else
217  {
218  G4ExceptionDescription ed;
219  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
220  << " The Navigator reports that it returned an invalid normal"
221  << G4endl;
222  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
223  EventMustBeAborted,ed,
224  "Invalid Surface Normal - Geometry must return valid surface normal");
225  }
226 
227  if (OldMomentum * theGlobalNormal > 0.0) {
228 #ifdef G4OPTICAL_DEBUG
229  G4ExceptionDescription ed;
230  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
231  << " theGlobalNormal points in a wrong direction. "
232  << G4endl;
233  ed << " The momentum of the photon arriving at interface (oldMomentum)"
234  << " must exit the volume cross in the step. " << G4endl;
235  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
236  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
237  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
238  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
239  ed << G4endl;
240  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
241  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
242  ed,
243  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
244 #else
245  theGlobalNormal = -theGlobalNormal;
246 #endif
247  }
248 
249  G4MaterialPropertiesTable* aMaterialPropertiesTable;
250  G4MaterialPropertyVector* RindexMPV = nullptr;
251 
252  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
253  if (aMaterialPropertiesTable) {
254  RindexMPV = aMaterialPropertiesTable->GetProperty(kRINDEX);
255  }
256  else {
257  theStatus = NoRINDEX;
258  if (verboseLevel > 0) BoundaryProcessVerbose();
259  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
260  aParticleChange.ProposeTrackStatus(fStopAndKill);
261  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
262  }
263 
264  if (RindexMPV) {
265  Rindex1 = RindexMPV->Value(thePhotonMomentum);
266  }
267  else {
268  theStatus = NoRINDEX;
269  if (verboseLevel > 0) BoundaryProcessVerbose();
270  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
271  aParticleChange.ProposeTrackStatus(fStopAndKill);
272  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
273  }
274 
275  theReflectivity = 1.;
276  theEfficiency = 0.;
277  theTransmittance = 0.;
278  theSurfaceRoughness = 0.;
279 
280  theModel = glisur;
281  theFinish = polished;
282 
283  G4SurfaceType type = dielectric_dielectric;
284 
285  RindexMPV = nullptr;
286  OpticalSurface = nullptr;
287 
288  G4LogicalSurface* Surface = nullptr;
289 
290  Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
291 
292  if (Surface == nullptr) {
293  G4bool enteredDaughter =
294  (thePostPV->GetMotherLogical() == thePrePV ->GetLogicalVolume());
295  if (enteredDaughter) {
296  Surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
297  if (Surface == nullptr) {
298  Surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
299  }
300  }
301  else {
302  Surface = G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
303  if (Surface == nullptr) {
304  Surface = G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
305  }
306  }
307  }
308 
309  if (Surface) {
310  OpticalSurface = dynamic_cast<G4OpticalSurface*> (Surface->GetSurfaceProperty());
311  }
312 
313  if (OpticalSurface) {
314  type = OpticalSurface->GetType();
315  theModel = OpticalSurface->GetModel();
316  theFinish = OpticalSurface->GetFinish();
317 
318  aMaterialPropertiesTable = OpticalSurface->GetMaterialPropertiesTable();
319 
320  if (aMaterialPropertiesTable) {
321 
322  if (theFinish == polishedbackpainted || theFinish == groundbackpainted) {
323  RindexMPV = aMaterialPropertiesTable->GetProperty(kRINDEX);
324  if (RindexMPV) {
325  Rindex2 = RindexMPV->Value(thePhotonMomentum);
326  }
327  else {
328  theStatus = NoRINDEX;
329  if (verboseLevel > 0) BoundaryProcessVerbose();
330  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
331  aParticleChange.ProposeTrackStatus(fStopAndKill);
332  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
333  }
334  }
335 
336  fRealRIndexMPV = aMaterialPropertiesTable->GetProperty(kREALRINDEX);
337  fImagRIndexMPV = aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX);
338 
339  iTE = 1;
340  iTM = 1;
341 
342  G4MaterialPropertyVector* PropertyPointer = aMaterialPropertiesTable->GetProperty(kREFLECTIVITY);
343  if (PropertyPointer) {
344  theReflectivity = PropertyPointer->Value(thePhotonMomentum);
345  } else if (fRealRIndexMPV && fImagRIndexMPV) {
346  CalculateReflectivity();
347  }
348 
349  PropertyPointer = aMaterialPropertiesTable->GetProperty(kEFFICIENCY);
350  if (PropertyPointer) {
351  theEfficiency = PropertyPointer->Value(thePhotonMomentum);
352  }
353 
354  PropertyPointer = aMaterialPropertiesTable->GetProperty(kTRANSMITTANCE);
355  if (PropertyPointer) {
356  theTransmittance = PropertyPointer->Value(thePhotonMomentum);
357  }
358 
359  if (aMaterialPropertiesTable->ConstPropertyExists("SURFACEROUGHNESS")) {
360  theSurfaceRoughness = aMaterialPropertiesTable-> GetConstProperty(kSURFACEROUGHNESS);
361  }
362 
363  if (theModel == unified) {
364  PropertyPointer = aMaterialPropertiesTable->GetProperty(kSPECULARLOBECONSTANT);
365  if (PropertyPointer) {
366  prob_sl = PropertyPointer->Value(thePhotonMomentum);
367  } else {
368  prob_sl = 0.0;
369  }
370 
371  PropertyPointer = aMaterialPropertiesTable->GetProperty(kSPECULARSPIKECONSTANT);
372  if (PropertyPointer) {
373  prob_ss = PropertyPointer->Value(thePhotonMomentum);
374  } else {
375  prob_ss = 0.0;
376  }
377 
378  PropertyPointer = aMaterialPropertiesTable->GetProperty(kBACKSCATTERCONSTANT);
379  if (PropertyPointer) {
380  prob_bs = PropertyPointer->Value(thePhotonMomentum);
381  } else {
382  prob_bs = 0.0;
383  }
384  }
385  } // end of if(aMaterialPropertiesTable)
386  else if (theFinish == polishedbackpainted || theFinish == groundbackpainted ) {
387  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
388  aParticleChange.ProposeTrackStatus(fStopAndKill);
389  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
390  }
391  } // end of if(OpticalSurface)
392 
393  // DIELECTRIC-DIELECTRIC
394  if (type == dielectric_dielectric) {
395  if (theFinish == polished || theFinish == ground) {
396  if (Material1 == Material2) {
397  theStatus = SameMaterial;
398  if (verboseLevel > 0) BoundaryProcessVerbose();
399  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
400  }
401  aMaterialPropertiesTable = Material2->GetMaterialPropertiesTable();
402  if (aMaterialPropertiesTable) RindexMPV = aMaterialPropertiesTable->GetProperty(kRINDEX);
403  if (RindexMPV) {
404  Rindex2 = RindexMPV->Value(thePhotonMomentum);
405  }
406  else {
407  theStatus = NoRINDEX;
408  if (verboseLevel > 0) BoundaryProcessVerbose();
409  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
410  aParticleChange.ProposeTrackStatus(fStopAndKill);
411  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
412  }
413  }
414  if (theFinish == polishedbackpainted || theFinish == groundbackpainted) {
415  DielectricDielectric();
416  }
417  else {
418  G4double rand = G4UniformRand();
419  if (rand > theReflectivity) {
420  if (rand > theReflectivity + theTransmittance) {
421  DoAbsorption();
422  } else {
423  theStatus = Transmission;
424  NewMomentum = OldMomentum;
425  NewPolarization = OldPolarization;
426  }
427  }
428  else {
429  if (theFinish == polishedfrontpainted) { DoReflection(); }
430  else if (theFinish == groundfrontpainted) {
431  theStatus = LambertianReflection;
432  DoReflection();
433  }
434  else { DielectricDielectric(); }
435  }
436  }
437  }
438 
439  else if (type == dielectric_metal) {
440  DielectricMetal();
441  }
442 
443  else if (type == dielectric_LUT) {
444  DielectricLUT();
445  }
446 
447  else if (type == dielectric_LUTDAVIS) {
448  DielectricLUTDAVIS();
449  }
450 
451  else if (type == dielectric_dichroic) {
452  DielectricDichroic();
453  }
454 
455  else {
456  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
457  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
458  }
459 
460  NewMomentum = NewMomentum.unit();
461  NewPolarization = NewPolarization.unit();
462 
463  if (verboseLevel > 0) {
464  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
465  G4cout << " New Polarization: " << NewPolarization << G4endl;
466  BoundaryProcessVerbose();
467  }
468 
469  aParticleChange.ProposeMomentumDirection(NewMomentum);
470  aParticleChange.ProposePolarization(NewPolarization);
471 
472  if (theStatus == FresnelRefraction || theStatus == Transmission ) {
473  G4MaterialPropertyVector* groupvel =
474  Material2->GetMaterialPropertiesTable()->GetProperty(kGROUPVEL);
475  if (groupvel) {
476  aParticleChange.ProposeVelocity(groupvel->Value(thePhotonMomentum));
477  }
478  }
479 
480  if (theStatus == Detection && fInvokeSD) InvokeSD(pStep);
481 
482  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
483 }
484 
485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
486 
487 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
488 {
489  if ( theStatus == Undefined )
490  G4cout << " *** Undefined *** " << G4endl;
491  if ( theStatus == Transmission )
492  G4cout << " *** Transmission *** " << G4endl;
493  if ( theStatus == FresnelRefraction )
494  G4cout << " *** FresnelRefraction *** " << G4endl;
495  if ( theStatus == FresnelReflection )
496  G4cout << " *** FresnelReflection *** " << G4endl;
497  if ( theStatus == TotalInternalReflection )
498  G4cout << " *** TotalInternalReflection *** " << G4endl;
499  if ( theStatus == LambertianReflection )
500  G4cout << " *** LambertianReflection *** " << G4endl;
501  if ( theStatus == LobeReflection )
502  G4cout << " *** LobeReflection *** " << G4endl;
503  if ( theStatus == SpikeReflection )
504  G4cout << " *** SpikeReflection *** " << G4endl;
505  if ( theStatus == BackScattering )
506  G4cout << " *** BackScattering *** " << G4endl;
507  if ( theStatus == PolishedLumirrorAirReflection )
508  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
509  if ( theStatus == PolishedLumirrorGlueReflection )
510  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
511  if ( theStatus == PolishedAirReflection )
512  G4cout << " *** PolishedAirReflection *** " << G4endl;
513  if ( theStatus == PolishedTeflonAirReflection )
514  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
515  if ( theStatus == PolishedTiOAirReflection )
516  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
517  if ( theStatus == PolishedTyvekAirReflection )
518  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
519  if ( theStatus == PolishedVM2000AirReflection )
520  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
521  if ( theStatus == PolishedVM2000GlueReflection )
522  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
523  if ( theStatus == EtchedLumirrorAirReflection )
524  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
525  if ( theStatus == EtchedLumirrorGlueReflection )
526  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
527  if ( theStatus == EtchedAirReflection )
528  G4cout << " *** EtchedAirReflection *** " << G4endl;
529  if ( theStatus == EtchedTeflonAirReflection )
530  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
531  if ( theStatus == EtchedTiOAirReflection )
532  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
533  if ( theStatus == EtchedTyvekAirReflection )
534  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
535  if ( theStatus == EtchedVM2000AirReflection )
536  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
537  if ( theStatus == EtchedVM2000GlueReflection )
538  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
539  if ( theStatus == GroundLumirrorAirReflection )
540  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
541  if ( theStatus == GroundLumirrorGlueReflection )
542  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
543  if ( theStatus == GroundAirReflection )
544  G4cout << " *** GroundAirReflection *** " << G4endl;
545  if ( theStatus == GroundTeflonAirReflection )
546  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
547  if ( theStatus == GroundTiOAirReflection )
548  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
549  if ( theStatus == GroundTyvekAirReflection )
550  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
551  if ( theStatus == GroundVM2000AirReflection )
552  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
553  if ( theStatus == GroundVM2000GlueReflection )
554  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
555  if ( theStatus == Absorption )
556  G4cout << " *** Absorption *** " << G4endl;
557  if ( theStatus == Detection )
558  G4cout << " *** Detection *** " << G4endl;
559  if ( theStatus == NotAtBoundary )
560  G4cout << " *** NotAtBoundary *** " << G4endl;
561  if ( theStatus == SameMaterial )
562  G4cout << " *** SameMaterial *** " << G4endl;
563  if ( theStatus == StepTooSmall )
564  G4cout << " *** StepTooSmall *** " << G4endl;
565  if ( theStatus == NoRINDEX )
566  G4cout << " *** NoRINDEX *** " << G4endl;
567  if ( theStatus == Dichroic )
568  G4cout << " *** Dichroic Transmission *** " << G4endl;
569 }
570 
571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
572 
573 G4ThreeVector
574 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
575  const G4ThreeVector& Normal ) const
576 {
577  G4ThreeVector FacetNormal;
578 
579  if (theModel == unified || theModel == LUT || theModel== DAVIS) {
580 
581  /* This function code alpha to a random value taken from the
582  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
583  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
584  is a gaussian distribution with mean 0 and standard deviation
585  sigma_alpha. */
586 
587  G4double alpha;
588 
589  G4double sigma_alpha = 0.0;
590  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
591 
592  if (sigma_alpha == 0.0) return FacetNormal = Normal;
593 
594  G4double f_max = std::min(1.0, 4.*sigma_alpha);
595 
596  G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
597  G4ThreeVector tmpNormal;
598 
599  do {
600  do {
601  alpha = G4RandGauss::shoot(0.0, sigma_alpha);
602  // Loop checking, 13-Aug-2015, Peter Gumplinger
603  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi);
604 
605  phi = G4UniformRand()*twopi;
606 
607  SinAlpha = std::sin(alpha);
608  CosAlpha = std::cos(alpha);
609  SinPhi = std::sin(phi);
610  CosPhi = std::cos(phi);
611 
612  unit_x = SinAlpha * CosPhi;
613  unit_y = SinAlpha * SinPhi;
614  unit_z = CosAlpha;
615 
616  FacetNormal.setX(unit_x);
617  FacetNormal.setY(unit_y);
618  FacetNormal.setZ(unit_z);
619 
620  tmpNormal = Normal;
621 
622  FacetNormal.rotateUz(tmpNormal);
623  // Loop checking, 13-Aug-2015, Peter Gumplinger
624  } while (Momentum * FacetNormal >= 0.0);
625  }
626  else {
627  G4double polish = 1.0;
628  if (OpticalSurface) polish = OpticalSurface->GetPolish();
629 
630  if (polish < 1.0) {
631  do {
632  G4ThreeVector smear;
633  do {
634  smear.setX(2.*G4UniformRand()-1.0);
635  smear.setY(2.*G4UniformRand()-1.0);
636  smear.setZ(2.*G4UniformRand()-1.0);
637  // Loop checking, 13-Aug-2015, Peter Gumplinger
638  } while (smear.mag() > 1.0);
639  smear = (1.-polish) * smear;
640  FacetNormal = Normal + smear;
641  // Loop checking, 13-Aug-2015, Peter Gumplinger
642  } while (Momentum * FacetNormal >= 0.0);
643  FacetNormal = FacetNormal.unit();
644  }
645  else {
646  FacetNormal = Normal;
647  }
648  }
649  return FacetNormal;
650 }
651 
652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
653 
654 void G4OpBoundaryProcess::DielectricMetal()
655 {
656  G4int n = 0;
657  G4double rand, PdotN, EdotN;
658  G4ThreeVector A_trans, A_paral;
659 
660  do {
661  n++;
662 
663  rand = G4UniformRand();
664  if (rand > theReflectivity && n == 1) {
665  if (rand > theReflectivity + theTransmittance) {
666  DoAbsorption();
667  } else {
668  theStatus = Transmission;
669  NewMomentum = OldMomentum;
670  NewPolarization = OldPolarization;
671  }
672  break;
673  }
674  else {
675  if (fRealRIndexMPV && fImagRIndexMPV) {
676  if (n > 1) {
677  CalculateReflectivity();
678  if (!G4BooleanRand(theReflectivity)) {
679  DoAbsorption();
680  break;
681  }
682  }
683  }
684 
685  if (theModel == glisur || theFinish == polished) {
686  DoReflection();
687  } else {
688  if (n == 1) ChooseReflection();
689  if (theStatus == LambertianReflection) {
690  DoReflection();
691  }
692  else if (theStatus == BackScattering) {
693  NewMomentum = -OldMomentum;
694  NewPolarization = -OldPolarization;
695  }
696  else {
697  if (theStatus == LobeReflection) {
698  if (fRealRIndexMPV && fImagRIndexMPV){
699  //
700  } else {
701  theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
702  }
703  }
704 
705  PdotN = OldMomentum * theFacetNormal;
706  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
707  EdotN = OldPolarization * theFacetNormal;
708 
709  if (sint1 > 0.0) {
710  A_trans = OldMomentum.cross(theFacetNormal);
711  A_trans = A_trans.unit();
712  } else {
713  A_trans = OldPolarization;
714  }
715  A_paral = NewMomentum.cross(A_trans);
716  A_paral = A_paral.unit();
717 
718  if (iTE > 0 && iTM > 0) {
719  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
720  } else if (iTE > 0) {
721  NewPolarization = -A_trans;
722  } else if (iTM >0 ) {
723  NewPolarization = -A_paral;
724  }
725  }
726  }
727 
728  OldMomentum = NewMomentum;
729  OldPolarization = NewPolarization;
730  }
731 
732  // Loop checking, 13-Aug-2015, Peter Gumplinger
733  } while (NewMomentum * theGlobalNormal < 0.0);
734 }
735 
736 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
737 
738 void G4OpBoundaryProcess::DielectricLUT()
739 {
740  G4int thetaIndex, phiIndex;
741  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
742  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
743 
744  theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
745  (G4int(NoRINDEX)-G4int(groundbackpainted)));
746 
747  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
748  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
749 
750  G4double rand;
751 
752  do {
753  rand = G4UniformRand();
754  if (rand > theReflectivity) {
755  if (rand > theReflectivity + theTransmittance) {
756  DoAbsorption();
757  } else {
758  theStatus = Transmission;
759  NewMomentum = OldMomentum;
760  NewPolarization = OldPolarization;
761  }
762  break;
763  }
764  else {
765  // Calculate Angle between Normal and Photon Momentum
766  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
767  // Round it to closest integer
768  G4int angleIncident = G4int(std::floor(180./pi*anglePhotonToNormal+0.5));
769 
770  // Take random angles THETA and PHI,
771  // and see if below Probability - if not - Redo
772  do {
773  thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
774  phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
775  // Find probability with the new indeces from LUT
776  AngularDistributionValue = OpticalSurface->
777  GetAngularDistributionValue(angleIncident, thetaIndex, phiIndex);
778  // Loop checking, 13-Aug-2015, Peter Gumplinger
779  } while (!G4BooleanRand(AngularDistributionValue));
780 
781  thetaRad = (-90 + 4*thetaIndex)*pi/180.;
782  phiRad = (-90 + 5*phiIndex)*pi/180.;
783  // Rotate Photon Momentum in Theta, then in Phi
784  NewMomentum = -OldMomentum;
785 
786  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
787  if (PerpendicularVectorTheta.mag() < kCarTolerance ) {
788  PerpendicularVectorTheta = NewMomentum.orthogonal();
789  }
790  NewMomentum = NewMomentum.rotate(anglePhotonToNormal-thetaRad, PerpendicularVectorTheta);
791  PerpendicularVectorPhi = PerpendicularVectorTheta.cross(NewMomentum);
792  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
793 
794  // Rotate Polarization too:
795  theFacetNormal = (NewMomentum - OldMomentum).unit();
796  EdotN = OldPolarization * theFacetNormal;
797  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
798  }
799  // Loop checking, 13-Aug-2015, Peter Gumplinger
800  } while (NewMomentum * theGlobalNormal <= 0.0);
801 }
802 
803 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
804 
805 void G4OpBoundaryProcess::DielectricLUTDAVIS()
806 {
807  G4int angindex, random, angleIncident;
808  G4double ReflectivityValue, elevation, azimuth, EdotN;
809  G4double anglePhotonToNormal;
810 
811  G4int LUTbin = OpticalSurface->GetLUTbins();
812  G4double rand = G4UniformRand();
813 
814  do {
815  anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
816  angleIncident = G4int(std::floor(180./pi*anglePhotonToNormal+0.5));
817 
818  ReflectivityValue = OpticalSurface->GetReflectivityLUTValue(angleIncident);
819 
820  if (rand > ReflectivityValue) {
821  if (theEfficiency > 0.) {
822  DoAbsorption();
823  break;
824  }
825  else {
826  theStatus = Transmission;
827 
828  if (angleIncident <= 0.01) {
829  NewMomentum = OldMomentum;
830  break;
831  }
832 
833  do {
834  random = G4RandFlat::shootInt(1, LUTbin+1);
835  angindex = (((random*2)-1)) + angleIncident*LUTbin*2 + 3640000;
836 
837  azimuth = OpticalSurface->GetAngularDistributionValueLUT(angindex-1);
838  elevation= OpticalSurface->GetAngularDistributionValueLUT(angindex);
839  } while (elevation == 0. && azimuth == 0.);
840 
841  NewMomentum = -OldMomentum;
842 
843  G4ThreeVector v = theGlobalNormal.cross(-NewMomentum);
844  G4ThreeVector vNorm = v/v.mag();
845  G4ThreeVector u = vNorm.cross(theGlobalNormal);
846 
847  u = u *= (std::sin(elevation) * std::cos(azimuth));
848  v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
849  G4ThreeVector w = theGlobalNormal *= (std::cos(elevation));
850  NewMomentum = G4ThreeVector(u+v+w);
851 
852  // Rotate Polarization too:
853  theFacetNormal = (NewMomentum - OldMomentum).unit();
854  EdotN = OldPolarization * theFacetNormal;
855  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
856  }
857  }
858  else {
859  theStatus = LobeReflection;
860 
861  if (angleIncident == 0) {
862  NewMomentum = -OldMomentum;
863  break;
864  }
865 
866  do {
867  random = G4RandFlat::shootInt(1, LUTbin+1);
868  angindex = (((random*2)-1)) + (angleIncident-1) * LUTbin*2;
869 
870  azimuth = OpticalSurface->GetAngularDistributionValueLUT(angindex-1);
871  elevation = OpticalSurface->GetAngularDistributionValueLUT(angindex);
872  } while (elevation == 0. && azimuth == 0.);
873 
874  NewMomentum = -OldMomentum;
875 
876  G4ThreeVector v = theGlobalNormal.cross(-NewMomentum);
877  G4ThreeVector vNorm = v/v.mag();
878  G4ThreeVector u = vNorm.cross(theGlobalNormal);
879 
880  u = u *= (std::sin(elevation) * std::cos(azimuth));
881  v = vNorm *= (std::sin(elevation) * std::sin(azimuth));
882  G4ThreeVector w = theGlobalNormal*=(std::cos(elevation));
883 
884  NewMomentum = G4ThreeVector(u+v+w);
885 
886  // Rotate Polarization too: (needs revision)
887  NewPolarization = OldPolarization;
888  }
889  } while (NewMomentum * theGlobalNormal <= 0.0);
890 }
891 
892 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
893 
894 void G4OpBoundaryProcess::DielectricDichroic()
895 {
896  // Calculate Angle between Normal and Photon Momentum
897  G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
898 
899  // Round it to closest integer
900  G4double angleIncident = std::floor(180./pi*anglePhotonToNormal+0.5);
901 
902  if (!DichroicVector) {
903  if (OpticalSurface) DichroicVector = OpticalSurface->GetDichroicVector();
904  }
905 
906  if (DichroicVector) {
907  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
908  theTransmittance = DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
909  // G4cout << "wavelength: " << std::floor(wavelength/nm)
910  // << "nm" << G4endl;
911  // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
912  // G4cout << "Transmittance: "
913  // << std::floor(theTransmittance/perCent) << "%" << G4endl;
914  } else {
915  G4ExceptionDescription ed;
916  ed << " G4OpBoundaryProcess/DielectricDichroic(): "
917  << " The dichroic surface has no G4Physics2DVector"
918  << G4endl;
919  G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
920  FatalException,ed,
921  "A dichroic surface must have an associated G4Physics2DVector");
922  }
923 
924  if (!G4BooleanRand(theTransmittance)) { // Not transmitted, so reflect
925  if (theModel == glisur || theFinish == polished) {
926  DoReflection();
927  } else {
928  ChooseReflection();
929  if (theStatus == LambertianReflection) {
930  DoReflection();
931  } else if (theStatus == BackScattering) {
932  NewMomentum = -OldMomentum;
933  NewPolarization = -OldPolarization;
934  } else {
935  G4double PdotN, EdotN;
936  do {
937  if (theStatus == LobeReflection) {
938  theFacetNormal = GetFacetNormal(OldMomentum, theGlobalNormal);
939  }
940  PdotN = OldMomentum * theFacetNormal;
941  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
942  // Loop checking, 13-Aug-2015, Peter Gumplinger
943  } while (NewMomentum * theGlobalNormal <= 0.0);
944 
945  EdotN = OldPolarization * theFacetNormal;
946  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
947  }
948  }
949  } else {
950  theStatus = Dichroic;
951  NewMomentum = OldMomentum;
952  NewPolarization = OldPolarization;
953  }
954 }
955 
956 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
957 
958 void G4OpBoundaryProcess::DielectricDielectric()
959 {
960  G4bool Inside = false;
961  G4bool Swap = false;
962 
963  G4bool SurfaceRoughnessCriterionPass = true;
964  if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) {
965  G4double wavelength = h_Planck*c_light/thePhotonMomentum;
966  G4double SurfaceRoughnessCriterion =
967  std::exp(-std::pow((4.*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2));
968  SurfaceRoughnessCriterionPass = G4BooleanRand(SurfaceRoughnessCriterion);
969  }
970 
971  leap:
972 
973  G4bool Through = false;
974  G4bool Done = false;
975 
976  G4double PdotN, EdotN;
977 
978  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
979  G4double E1_perp, E1_parl;
980  G4double s1, s2, E2_perp, E2_parl, E2_total, TransCoeff;
981  G4double E2_abs, C_parl, C_perp;
982  G4double alpha;
983 
984  do {
985  if (Through) {
986  Swap = !Swap;
987  Through = false;
988  theGlobalNormal = -theGlobalNormal;
989  G4SwapPtr(Material1, Material2);
990  G4SwapObj(&Rindex1, &Rindex2);
991  }
992 
993  if (theFinish == polished) {
994  theFacetNormal = theGlobalNormal;
995  }
996  else {
997  theFacetNormal = GetFacetNormal(OldMomentum, theGlobalNormal);
998  }
999 
1000  PdotN = OldMomentum * theFacetNormal;
1001  EdotN = OldPolarization * theFacetNormal;
1002 
1003  cost1 = -PdotN;
1004  if (std::abs(cost1) < 1.0-kCarTolerance) {
1005  sint1 = std::sqrt(1.-cost1*cost1);
1006  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
1007  // this isn't a sine as we might expect from the name; can be > 1
1008  }
1009  else {
1010  sint1 = 0.0;
1011  sint2 = 0.0;
1012  }
1013 
1014  // TOTAL INTERNAL REFLECTION
1015  if (sint2 >= 1.0) {
1016  // Simulate total internal reflection
1017 
1018  if (Swap) Swap = !Swap;
1019 
1020  theStatus = TotalInternalReflection;
1021 
1022  if (!SurfaceRoughnessCriterionPass) theStatus = LambertianReflection;
1023 
1024  if (theModel == unified && theFinish != polished) ChooseReflection();
1025 
1026  if (theStatus == LambertianReflection) {
1027  DoReflection();
1028  }
1029  else if (theStatus == BackScattering) {
1030  NewMomentum = -OldMomentum;
1031  NewPolarization = -OldPolarization;
1032  }
1033  else {
1034  PdotN = OldMomentum * theFacetNormal;
1035  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1036  EdotN = OldPolarization * theFacetNormal;
1037  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1038  }
1039  }
1040  // NOT TIR
1041  else if (sint2 < 1.0) {
1042  // Calculate amplitude for transmission (Q = P x N)
1043 
1044  if (cost1 > 0.0) {
1045  cost2 = std::sqrt(1.-sint2*sint2);
1046  }
1047  else {
1048  cost2 = -std::sqrt(1.-sint2*sint2);
1049  }
1050 
1051  if (sint1 > 0.0) {
1052  A_trans = OldMomentum.cross(theFacetNormal);
1053  A_trans = A_trans.unit();
1054  E1_perp = OldPolarization * A_trans;
1055  E1pp = E1_perp * A_trans;
1056  E1pl = OldPolarization - E1pp;
1057  E1_parl = E1pl.mag();
1058  }
1059  else {
1060  A_trans = OldPolarization;
1061  // Here we Follow Jackson's conventions and we set the
1062  // parallel component = 1 in case of a ray perpendicular
1063  // to the surface
1064  E1_perp = 0.0;
1065  E1_parl = 1.0;
1066  }
1067 
1068  s1 = Rindex1*cost1;
1069  E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
1070  E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
1071  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1072  s2 = Rindex2*cost2*E2_total;
1073 
1074  if (theTransmittance > 0.) TransCoeff = theTransmittance;
1075  else if (cost1 != 0.0) TransCoeff = s2/s1;
1076  else TransCoeff = 0.0;
1077 
1078  // NOT TIR: REFLECTION
1079  if (!G4BooleanRand(TransCoeff)) {
1080  // Simulate reflection
1081 
1082  if (Swap) Swap = !Swap;
1083 
1084  theStatus = FresnelReflection;
1085 
1086  if (!SurfaceRoughnessCriterionPass) theStatus = LambertianReflection;
1087 
1088  if (theModel == unified && theFinish != polished) ChooseReflection();
1089 
1090  if (theStatus == LambertianReflection) {
1091  DoReflection();
1092  }
1093  else if (theStatus == BackScattering) {
1094  NewMomentum = -OldMomentum;
1095  NewPolarization = -OldPolarization;
1096  }
1097  else {
1098  PdotN = OldMomentum * theFacetNormal;
1099  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1100 
1101  if (sint1 > 0.0) { // incident ray oblique
1102  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
1103  E2_perp = E2_perp - E1_perp;
1104  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
1105  A_paral = NewMomentum.cross(A_trans);
1106  A_paral = A_paral.unit();
1107  E2_abs = std::sqrt(E2_total);
1108  C_parl = E2_parl/E2_abs;
1109  C_perp = E2_perp/E2_abs;
1110 
1111  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1112  }
1113  else { // incident ray perpendicular
1114  if (Rindex2 > Rindex1) {
1115  NewPolarization = - OldPolarization;
1116  }
1117  else {
1118  NewPolarization = OldPolarization;
1119  }
1120  }
1121  }
1122  }
1123  // NOT TIR: TRANSMISSION
1124  else { // photon gets transmitted
1125  // Simulate transmission/refraction
1126 
1127  Inside = !Inside;
1128  Through = true;
1129  theStatus = FresnelRefraction;
1130 
1131  if (sint1 > 0.0) { // incident ray oblique
1132  alpha = cost1 - cost2*(Rindex2/Rindex1);
1133  NewMomentum = OldMomentum + alpha*theFacetNormal;
1134  NewMomentum = NewMomentum.unit();
1135  A_paral = NewMomentum.cross(A_trans);
1136  A_paral = A_paral.unit();
1137  E2_abs = std::sqrt(E2_total);
1138  C_parl = E2_parl/E2_abs;
1139  C_perp = E2_perp/E2_abs;
1140 
1141  NewPolarization = C_parl*A_paral + C_perp*A_trans;
1142  }
1143  else { // incident ray perpendicular
1144  NewMomentum = OldMomentum;
1145  NewPolarization = OldPolarization;
1146  }
1147  }
1148  }
1149 
1150  OldMomentum = NewMomentum.unit();
1151  OldPolarization = NewPolarization.unit();
1152 
1153  if (theStatus == FresnelRefraction) {
1154  Done = (NewMomentum * theGlobalNormal <= 0.0);
1155  }
1156  else {
1157  Done = (NewMomentum * theGlobalNormal >= -kCarTolerance);
1158  }
1159 
1160  // Loop checking, 13-Aug-2015, Peter Gumplinger
1161  } while (!Done);
1162 
1163  if (Inside && !Swap) {
1164  if (theFinish == polishedbackpainted || theFinish == groundbackpainted) {
1165  G4double rand = G4UniformRand();
1166  if (rand > theReflectivity) {
1167  if (rand > theReflectivity + theTransmittance) {
1168  DoAbsorption();
1169  } else {
1170  theStatus = Transmission;
1171  NewMomentum = OldMomentum;
1172  NewPolarization = OldPolarization;
1173  }
1174  }
1175  else {
1176  if (theStatus != FresnelRefraction ) {
1177  theGlobalNormal = -theGlobalNormal;
1178  }
1179  else {
1180  Swap = !Swap;
1181  G4SwapPtr(Material1, Material2);
1182  G4SwapObj(&Rindex1, &Rindex2);
1183  }
1184  if (theFinish == groundbackpainted) theStatus = LambertianReflection;
1185 
1186  DoReflection();
1187 
1188  theGlobalNormal = -theGlobalNormal;
1189  OldMomentum = NewMomentum;
1190 
1191  goto leap;
1192  }
1193  }
1194  }
1195 }
1196 
1197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1198 
1199 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
1200  G4double ,
1201  G4ForceCondition* condition)
1202 {
1203  *condition = Forced;
1204  return DBL_MAX;
1205 }
1206 
1207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1208 
1209 G4double G4OpBoundaryProcess::GetIncidentAngle()
1210 {
1211  G4double PdotN = OldMomentum * theFacetNormal;
1212  G4double magP = OldMomentum.mag();
1213  G4double magN = theFacetNormal.mag();
1214  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1215 
1216  return incidentangle;
1217 }
1218 
1219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1220 
1221 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1222  G4double E1_parl,
1223  G4double incidentangle,
1224  G4double RealRindex,
1225  G4double ImaginaryRindex)
1226 {
1227  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1228  G4complex N1(Rindex1, 0.), N2(RealRindex, ImaginaryRindex);
1229  G4complex CosPhi;
1230 
1231  G4complex u(1.,0.); //unit number 1
1232 
1233  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1234  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1235  G4complex denominatorTE, denominatorTM;
1236  G4complex rTM, rTE;
1237 
1238  G4MaterialPropertiesTable* aMaterialPropertiesTable =
1239  Material1->GetMaterialPropertiesTable();
1240  G4MaterialPropertyVector* aPropertyPointerR =
1241  aMaterialPropertiesTable->GetProperty(kREALRINDEX);
1242  G4MaterialPropertyVector* aPropertyPointerI =
1243  aMaterialPropertiesTable->GetProperty(kIMAGINARYRINDEX);
1244  if (aPropertyPointerR && aPropertyPointerI) {
1245  G4double RRindex = aPropertyPointerR->Value(thePhotonMomentum);
1246  G4double IRindex = aPropertyPointerI->Value(thePhotonMomentum);
1247  N1 = G4complex(RRindex,IRindex);
1248  }
1249 
1250  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1251  // Optics" written by Fowles
1252 
1253  CosPhi = std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))*(N1*N1)/(N2*N2)));
1254 
1255  numeratorTE = N1*std::cos(incidentangle) - N2*CosPhi;
1256  denominatorTE = N1*std::cos(incidentangle) + N2*CosPhi;
1257  rTE = numeratorTE/denominatorTE;
1258 
1259  numeratorTM = N2*std::cos(incidentangle) - N1*CosPhi;
1260  denominatorTM = N2*std::cos(incidentangle) + N1*CosPhi;
1261  rTM = numeratorTM/denominatorTM;
1262 
1263  // This is my calculaton for reflectivity on a metalic surface
1264  // depending on the fraction of TE and TM polarization
1265  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1266  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1267 
1268  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1269  / (E1_perp*E1_perp + E1_parl*E1_parl);
1270  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1271  / (E1_perp*E1_perp + E1_parl*E1_parl);
1272  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1273 
1274  do {
1275  if (G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE)) {
1276  iTE = -1;
1277  } else {
1278  iTE = 1;
1279  }
1280  if (G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM)) {
1281  iTM = -1;
1282  } else {
1283  iTM = 1;
1284  }
1285  // Loop checking, 13-Aug-2015, Peter Gumplinger
1286  } while (iTE < 0 && iTM < 0);
1287 
1288  return real(Reflectivity);
1289 }
1290 
1291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1292 
1293 void G4OpBoundaryProcess::CalculateReflectivity()
1294 {
1295  G4double RealRindex = fRealRIndexMPV->Value(thePhotonMomentum);
1296  G4double ImaginaryRindex = fImagRIndexMPV->Value(thePhotonMomentum);
1297 
1298  // calculate FacetNormal
1299  if (theFinish == ground) {
1300  theFacetNormal = GetFacetNormal(OldMomentum, theGlobalNormal);
1301  } else {
1302  theFacetNormal = theGlobalNormal;
1303  }
1304 
1305  G4double PdotN = OldMomentum * theFacetNormal;
1306  cost1 = -PdotN;
1307 
1308  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1309  sint1 = std::sqrt(1. - cost1*cost1);
1310  } else {
1311  sint1 = 0.0;
1312  }
1313 
1314  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1315  G4double E1_perp, E1_parl;
1316 
1317  if (sint1 > 0.0) {
1318  A_trans = OldMomentum.cross(theFacetNormal);
1319  A_trans = A_trans.unit();
1320  E1_perp = OldPolarization * A_trans;
1321  E1pp = E1_perp * A_trans;
1322  E1pl = OldPolarization - E1pp;
1323  E1_parl = E1pl.mag();
1324  }
1325  else {
1326  A_trans = OldPolarization;
1327  // Here we Follow Jackson's conventions and we set the
1328  // parallel component = 1 in case of a ray perpendicular
1329  // to the surface
1330  E1_perp = 0.0;
1331  E1_parl = 1.0;
1332  }
1333 
1334  //calculate incident angle
1335  G4double incidentangle = GetIncidentAngle();
1336 
1337  //calculate the reflectivity depending on incident angle,
1338  //polarization and complex refractive
1339  theReflectivity = GetReflectivity(E1_perp, E1_parl, incidentangle,
1340  RealRindex, ImaginaryRindex);
1341 }
1342 
1343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1344 
1345 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1346 {
1347  G4Step aStep = *pStep;
1348  aStep.AddTotalEnergyDeposit(thePhotonMomentum);
1349 
1350  G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1351  if (sd) return sd->Hit(&aStep);
1352  else return false;
1353 }