Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 //
67 // Author: Peter Gumplinger
68 // adopted from work by Werner Keil - April 2/96
69 // mail: gum@triumf.ca
70 //
72 
73 #include "G4ios.hh"
74 #include "G4PhysicalConstants.hh"
75 #include "G4OpProcessSubType.hh"
76 
77 #include "G4OpBoundaryProcess.hh"
78 #include "G4GeometryTolerance.hh"
79 
81 // Class Implementation
83 
85  // Operators
87 
88 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
89 // {
90 // }
91 
93  // Constructors
95 
96 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
97  G4ProcessType type)
98  : G4VDiscreteProcess(processName, type)
99 {
100  if ( verboseLevel > 0) {
101  G4cout << GetProcessName() << " is created " << G4endl;
102  }
103 
104  SetProcessSubType(fOpBoundary);
105 
106  theStatus = Undefined;
107  theModel = glisur;
108  theFinish = polished;
109  theReflectivity = 1.;
110  theEfficiency = 0.;
111  theTransmittance = 0.;
112 
113  prob_sl = 0.;
114  prob_ss = 0.;
115  prob_bs = 0.;
116 
117  PropertyPointer = NULL;
118  PropertyPointer1 = NULL;
119  PropertyPointer2 = NULL;
120 
121  Material1 = NULL;
122  Material2 = NULL;
123 
124  OpticalSurface = NULL;
125 
126  kCarTolerance = G4GeometryTolerance::GetInstance()
127  ->GetSurfaceTolerance();
128 
129  iTE = iTM = 0;
130  thePhotonMomentum = 0.;
131  Rindex1 = Rindex2 = cost1 = cost2 = sint1 = sint2 = 0.;
132 }
133 
134 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
135 // {
136 // }
137 
139  // Destructors
141 
142 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
143 
145  // Methods
147 
148 // PostStepDoIt
149 // ------------
150 //
151 G4VParticleChange*
152 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
153 {
154  theStatus = Undefined;
155 
156  aParticleChange.Initialize(aTrack);
157  aParticleChange.ProposeVelocity(aTrack.GetVelocity());
158 
159  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
160  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
161 
162  if ( verboseLevel > 0 ) {
163  G4cout << " Photon at Boundary! " << G4endl;
164  G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume();
165  G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume();
166  if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
167  if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
168  }
169 
170  if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
171  theStatus = NotAtBoundary;
172  if ( verboseLevel > 0) BoundaryProcessVerbose();
173  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
174  }
175  if (aTrack.GetStepLength()<=kCarTolerance/2){
176  theStatus = StepTooSmall;
177  if ( verboseLevel > 0) BoundaryProcessVerbose();
178  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
179  }
180 
181  Material1 = pPreStepPoint -> GetMaterial();
182  Material2 = pPostStepPoint -> GetMaterial();
183 
184  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
185 
186  thePhotonMomentum = aParticle->GetTotalMomentum();
187  OldMomentum = aParticle->GetMomentumDirection();
188  OldPolarization = aParticle->GetPolarization();
189 
190  if ( verboseLevel > 0 ) {
191  G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
192  G4cout << " Old Polarization: " << OldPolarization << G4endl;
193  }
194 
195  G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
196 
197  G4Navigator* theNavigator =
198  G4TransportationManager::GetTransportationManager()->
199  GetNavigatorForTracking();
200 
201  G4bool valid;
202  // Use the new method for Exit Normal in global coordinates,
203  // which provides the normal more reliably.
204  theGlobalNormal =
205  theNavigator->GetGlobalExitNormal(theGlobalPoint,&valid);
206 
207  if (valid) {
208  theGlobalNormal = -theGlobalNormal;
209  }
210  else
211  {
212  G4ExceptionDescription ed;
213  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
214  << " The Navigator reports that it returned an invalid normal"
215  << G4endl;
216  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
217  EventMustBeAborted,ed,
218  "Invalid Surface Normal - Geometry must return valid surface normal");
219  }
220 
221  if (OldMomentum * theGlobalNormal > 0.0) {
222 #ifdef G4OPTICAL_DEBUG
223  G4ExceptionDescription ed;
224  ed << " G4OpBoundaryProcess/PostStepDoIt(): "
225  << " theGlobalNormal points in a wrong direction. "
226  << G4endl;
227  ed << " The momentum of the photon arriving at interface (oldMomentum)"
228  << " must exit the volume cross in the step. " << G4endl;
229  ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
230  ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
231  ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
232  ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
233  ed << G4endl;
234  G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
235  EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
236  ed,
237  "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
238 #else
239  theGlobalNormal = -theGlobalNormal;
240 #endif
241  }
242 
243  G4MaterialPropertiesTable* aMaterialPropertiesTable;
244  G4MaterialPropertyVector* Rindex;
245 
246  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
247  if (aMaterialPropertiesTable) {
248  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
249  }
250  else {
251  theStatus = NoRINDEX;
252  if ( verboseLevel > 0) BoundaryProcessVerbose();
253  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
254  aParticleChange.ProposeTrackStatus(fStopAndKill);
255  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
256  }
257 
258  if (Rindex) {
259  Rindex1 = Rindex->Value(thePhotonMomentum);
260  }
261  else {
262  theStatus = NoRINDEX;
263  if ( verboseLevel > 0) BoundaryProcessVerbose();
264  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
265  aParticleChange.ProposeTrackStatus(fStopAndKill);
266  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
267  }
268 
269  theReflectivity = 1.;
270  theEfficiency = 0.;
271  theTransmittance = 0.;
272 
273  theModel = glisur;
274  theFinish = polished;
275 
276  G4SurfaceType type = dielectric_dielectric;
277 
278  Rindex = NULL;
279  OpticalSurface = NULL;
280 
281  G4LogicalSurface* Surface = NULL;
282 
283  Surface = G4LogicalBorderSurface::GetSurface
284  (pPreStepPoint ->GetPhysicalVolume(),
285  pPostStepPoint->GetPhysicalVolume());
286 
287  if (Surface == NULL){
288  G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
289  ->GetMotherLogical() ==
290  pPreStepPoint->GetPhysicalVolume()
291  ->GetLogicalVolume());
292  if(enteredDaughter){
293  Surface = G4LogicalSkinSurface::GetSurface
294  (pPostStepPoint->GetPhysicalVolume()->
295  GetLogicalVolume());
296  if(Surface == NULL)
297  Surface = G4LogicalSkinSurface::GetSurface
298  (pPreStepPoint->GetPhysicalVolume()->
299  GetLogicalVolume());
300  }
301  else {
302  Surface = G4LogicalSkinSurface::GetSurface
303  (pPreStepPoint->GetPhysicalVolume()->
304  GetLogicalVolume());
305  if(Surface == NULL)
306  Surface = G4LogicalSkinSurface::GetSurface
307  (pPostStepPoint->GetPhysicalVolume()->
308  GetLogicalVolume());
309  }
310  }
311 
312  if (Surface) OpticalSurface =
313  dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
314 
315  if (OpticalSurface) {
316 
317  type = OpticalSurface->GetType();
318  theModel = OpticalSurface->GetModel();
319  theFinish = OpticalSurface->GetFinish();
320 
321  aMaterialPropertiesTable = OpticalSurface->
322  GetMaterialPropertiesTable();
323 
324  if (aMaterialPropertiesTable) {
325 
326  if (theFinish == polishedbackpainted ||
327  theFinish == groundbackpainted ) {
328  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
329  if (Rindex) {
330  Rindex2 = Rindex->Value(thePhotonMomentum);
331  }
332  else {
333  theStatus = NoRINDEX;
334  if ( verboseLevel > 0) BoundaryProcessVerbose();
335  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
336  aParticleChange.ProposeTrackStatus(fStopAndKill);
337  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
338  }
339  }
340 
341  PropertyPointer =
342  aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
343  PropertyPointer1 =
344  aMaterialPropertiesTable->GetProperty("REALRINDEX");
345  PropertyPointer2 =
346  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
347 
348  iTE = 1;
349  iTM = 1;
350 
351  if (PropertyPointer) {
352 
353  theReflectivity =
354  PropertyPointer->Value(thePhotonMomentum);
355 
356  } else if (PropertyPointer1 && PropertyPointer2) {
357 
358  CalculateReflectivity();
359 
360  }
361 
362  PropertyPointer =
363  aMaterialPropertiesTable->GetProperty("EFFICIENCY");
364  if (PropertyPointer) {
365  theEfficiency =
366  PropertyPointer->Value(thePhotonMomentum);
367  }
368 
369  PropertyPointer =
370  aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
371  if (PropertyPointer) {
372  theTransmittance =
373  PropertyPointer->Value(thePhotonMomentum);
374  }
375 
376  if ( theModel == unified ) {
377  PropertyPointer =
378  aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
379  if (PropertyPointer) {
380  prob_sl =
381  PropertyPointer->Value(thePhotonMomentum);
382  } else {
383  prob_sl = 0.0;
384  }
385 
386  PropertyPointer =
387  aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
388  if (PropertyPointer) {
389  prob_ss =
390  PropertyPointer->Value(thePhotonMomentum);
391  } else {
392  prob_ss = 0.0;
393  }
394 
395  PropertyPointer =
396  aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
397  if (PropertyPointer) {
398  prob_bs =
399  PropertyPointer->Value(thePhotonMomentum);
400  } else {
401  prob_bs = 0.0;
402  }
403  }
404  }
405  else if (theFinish == polishedbackpainted ||
406  theFinish == groundbackpainted ) {
407  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
408  aParticleChange.ProposeTrackStatus(fStopAndKill);
409  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
410  }
411  }
412 
413  if (type == dielectric_dielectric ) {
414  if (theFinish == polished || theFinish == ground ) {
415 
416  if (Material1 == Material2){
417  theStatus = SameMaterial;
418  if ( verboseLevel > 0) BoundaryProcessVerbose();
419  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
420  }
421  aMaterialPropertiesTable =
422  Material2->GetMaterialPropertiesTable();
423  if (aMaterialPropertiesTable)
424  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
425  if (Rindex) {
426  Rindex2 = Rindex->Value(thePhotonMomentum);
427  }
428  else {
429  theStatus = NoRINDEX;
430  if ( verboseLevel > 0) BoundaryProcessVerbose();
431  aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
432  aParticleChange.ProposeTrackStatus(fStopAndKill);
433  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
434  }
435  }
436  }
437 
438  if (type == dielectric_metal) {
439 
440  DielectricMetal();
441 
442  // Uncomment the following lines if you wish to have
443  // Transmission instead of Absorption
444  // if (theStatus == Absorption) {
445  // return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
446  // }
447 
448  }
449  else if (type == dielectric_LUT) {
450 
451  DielectricLUT();
452 
453  }
454  else if (type == dielectric_dielectric) {
455 
456  if ( theFinish == polishedbackpainted ||
457  theFinish == groundbackpainted ) {
458  DielectricDielectric();
459  }
460  else {
461  if ( !G4BooleanRand(theReflectivity) ) {
462  DoAbsorption();
463  }
464  else {
465  if ( theFinish == polishedfrontpainted ) {
466  DoReflection();
467  }
468  else if ( theFinish == groundfrontpainted ) {
469  theStatus = LambertianReflection;
470  DoReflection();
471  }
472  else {
473  DielectricDielectric();
474  }
475  }
476  }
477  }
478  else {
479 
480  G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
481  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
482 
483  }
484 
485  NewMomentum = NewMomentum.unit();
486  NewPolarization = NewPolarization.unit();
487 
488  if ( verboseLevel > 0) {
489  G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
490  G4cout << " New Polarization: " << NewPolarization << G4endl;
491  BoundaryProcessVerbose();
492  }
493 
494  aParticleChange.ProposeMomentumDirection(NewMomentum);
495  aParticleChange.ProposePolarization(NewPolarization);
496 
497  if ( theStatus == FresnelRefraction ) {
498  G4MaterialPropertyVector* groupvel =
499  Material2->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
500  G4double finalVelocity = groupvel->Value(thePhotonMomentum);
501  aParticleChange.ProposeVelocity(finalVelocity);
502  }
503 
504  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
505 }
506 
507 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
508 {
509  if ( theStatus == Undefined )
510  G4cout << " *** Undefined *** " << G4endl;
511  if ( theStatus == FresnelRefraction )
512  G4cout << " *** FresnelRefraction *** " << G4endl;
513  if ( theStatus == FresnelReflection )
514  G4cout << " *** FresnelReflection *** " << G4endl;
515  if ( theStatus == TotalInternalReflection )
516  G4cout << " *** TotalInternalReflection *** " << G4endl;
517  if ( theStatus == LambertianReflection )
518  G4cout << " *** LambertianReflection *** " << G4endl;
519  if ( theStatus == LobeReflection )
520  G4cout << " *** LobeReflection *** " << G4endl;
521  if ( theStatus == SpikeReflection )
522  G4cout << " *** SpikeReflection *** " << G4endl;
523  if ( theStatus == BackScattering )
524  G4cout << " *** BackScattering *** " << G4endl;
525  if ( theStatus == PolishedLumirrorAirReflection )
526  G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
527  if ( theStatus == PolishedLumirrorGlueReflection )
528  G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
529  if ( theStatus == PolishedAirReflection )
530  G4cout << " *** PolishedAirReflection *** " << G4endl;
531  if ( theStatus == PolishedTeflonAirReflection )
532  G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
533  if ( theStatus == PolishedTiOAirReflection )
534  G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
535  if ( theStatus == PolishedTyvekAirReflection )
536  G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
537  if ( theStatus == PolishedVM2000AirReflection )
538  G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
539  if ( theStatus == PolishedVM2000GlueReflection )
540  G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
541  if ( theStatus == EtchedLumirrorAirReflection )
542  G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
543  if ( theStatus == EtchedLumirrorGlueReflection )
544  G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
545  if ( theStatus == EtchedAirReflection )
546  G4cout << " *** EtchedAirReflection *** " << G4endl;
547  if ( theStatus == EtchedTeflonAirReflection )
548  G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
549  if ( theStatus == EtchedTiOAirReflection )
550  G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
551  if ( theStatus == EtchedTyvekAirReflection )
552  G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
553  if ( theStatus == EtchedVM2000AirReflection )
554  G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
555  if ( theStatus == EtchedVM2000GlueReflection )
556  G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
557  if ( theStatus == GroundLumirrorAirReflection )
558  G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
559  if ( theStatus == GroundLumirrorGlueReflection )
560  G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
561  if ( theStatus == GroundAirReflection )
562  G4cout << " *** GroundAirReflection *** " << G4endl;
563  if ( theStatus == GroundTeflonAirReflection )
564  G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
565  if ( theStatus == GroundTiOAirReflection )
566  G4cout << " *** GroundTiOAirReflection *** " << G4endl;
567  if ( theStatus == GroundTyvekAirReflection )
568  G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
569  if ( theStatus == GroundVM2000AirReflection )
570  G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
571  if ( theStatus == GroundVM2000GlueReflection )
572  G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
573  if ( theStatus == Absorption )
574  G4cout << " *** Absorption *** " << G4endl;
575  if ( theStatus == Detection )
576  G4cout << " *** Detection *** " << G4endl;
577  if ( theStatus == NotAtBoundary )
578  G4cout << " *** NotAtBoundary *** " << G4endl;
579  if ( theStatus == SameMaterial )
580  G4cout << " *** SameMaterial *** " << G4endl;
581  if ( theStatus == StepTooSmall )
582  G4cout << " *** StepTooSmall *** " << G4endl;
583  if ( theStatus == NoRINDEX )
584  G4cout << " *** NoRINDEX *** " << G4endl;
585 }
586 
587 G4ThreeVector
588 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
589  const G4ThreeVector& Normal ) const
590 {
591  G4ThreeVector FacetNormal;
592 
593  if (theModel == unified || theModel == LUT) {
594 
595  /* This function code alpha to a random value taken from the
596  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
597  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
598  is a gaussian distribution with mean 0 and standard deviation
599  sigma_alpha. */
600 
601  G4double alpha;
602 
603  G4double sigma_alpha = 0.0;
604  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
605 
606  if (sigma_alpha == 0.0) return FacetNormal = Normal;
607 
608  G4double f_max = std::min(1.0,4.*sigma_alpha);
609 
610  do {
611  do {
612  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
613  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
614 
615  G4double phi = G4UniformRand()*twopi;
616 
617  G4double SinAlpha = std::sin(alpha);
618  G4double CosAlpha = std::cos(alpha);
619  G4double SinPhi = std::sin(phi);
620  G4double CosPhi = std::cos(phi);
621 
622  G4double unit_x = SinAlpha * CosPhi;
623  G4double unit_y = SinAlpha * SinPhi;
624  G4double unit_z = CosAlpha;
625 
626  FacetNormal.setX(unit_x);
627  FacetNormal.setY(unit_y);
628  FacetNormal.setZ(unit_z);
629 
630  G4ThreeVector tmpNormal = Normal;
631 
632  FacetNormal.rotateUz(tmpNormal);
633  } while (Momentum * FacetNormal >= 0.0);
634  }
635  else {
636 
637  G4double polish = 1.0;
638  if (OpticalSurface) polish = OpticalSurface->GetPolish();
639 
640  if (polish < 1.0) {
641  do {
642  G4ThreeVector smear;
643  do {
644  smear.setX(2.*G4UniformRand()-1.0);
645  smear.setY(2.*G4UniformRand()-1.0);
646  smear.setZ(2.*G4UniformRand()-1.0);
647  } while (smear.mag()>1.0);
648  smear = (1.-polish) * smear;
649  FacetNormal = Normal + smear;
650  } while (Momentum * FacetNormal >= 0.0);
651  FacetNormal = FacetNormal.unit();
652  }
653  else {
654  FacetNormal = Normal;
655  }
656  }
657  return FacetNormal;
658 }
659 
660 void G4OpBoundaryProcess::DielectricMetal()
661 {
662  G4int n = 0;
663 
664  do {
665 
666  n++;
667 
668  if( !G4BooleanRand(theReflectivity) && n == 1 ) {
669 
670  // Comment out DoAbsorption and uncomment theStatus = Absorption;
671  // if you wish to have Transmission instead of Absorption
672 
673  DoAbsorption();
674  // theStatus = Absorption;
675  break;
676 
677  }
678  else {
679 
680  if (PropertyPointer1 && PropertyPointer2) {
681  if ( n > 1 ) {
682  CalculateReflectivity();
683  if ( !G4BooleanRand(theReflectivity) ) {
684  DoAbsorption();
685  break;
686  }
687  }
688  }
689 
690  if ( theModel == glisur || theFinish == polished ) {
691 
692  DoReflection();
693 
694  } else {
695 
696  if ( n == 1 ) ChooseReflection();
697 
698  if ( theStatus == LambertianReflection ) {
699  DoReflection();
700  }
701  else if ( theStatus == BackScattering ) {
702  NewMomentum = -OldMomentum;
703  NewPolarization = -OldPolarization;
704  }
705  else {
706 
707  if(theStatus==LobeReflection){
708  if ( PropertyPointer1 && PropertyPointer2 ){
709  } else {
710  theFacetNormal =
711  GetFacetNormal(OldMomentum,theGlobalNormal);
712  }
713  }
714 
715  G4double PdotN = OldMomentum * theFacetNormal;
716  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
717  G4double EdotN = OldPolarization * theFacetNormal;
718 
719  G4ThreeVector A_trans, A_paral;
720 
721  if (sint1 > 0.0 ) {
722  A_trans = OldMomentum.cross(theFacetNormal);
723  A_trans = A_trans.unit();
724  } else {
725  A_trans = OldPolarization;
726  }
727  A_paral = NewMomentum.cross(A_trans);
728  A_paral = A_paral.unit();
729 
730  if(iTE>0&&iTM>0) {
731  NewPolarization =
732  -OldPolarization + (2.*EdotN)*theFacetNormal;
733  } else if (iTE>0) {
734  NewPolarization = -A_trans;
735  } else if (iTM>0) {
736  NewPolarization = -A_paral;
737  }
738 
739  }
740 
741  }
742 
743  OldMomentum = NewMomentum;
744  OldPolarization = NewPolarization;
745 
746  }
747 
748  } while (NewMomentum * theGlobalNormal < 0.0);
749 }
750 
751 void G4OpBoundaryProcess::DielectricLUT()
752 {
753  G4int thetaIndex, phiIndex;
754  G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
755  G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
756 
757  theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
758  (G4int(NoRINDEX)-G4int(groundbackpainted)));
759 
760  G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
761  G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
762 
763  do {
764  if ( !G4BooleanRand(theReflectivity) ) // Not reflected, so Absorbed
765  DoAbsorption();
766  else {
767  // Calculate Angle between Normal and Photon Momentum
768  G4double anglePhotonToNormal =
769  OldMomentum.angle(-theGlobalNormal);
770  // Round it to closest integer
771  G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
772 
773  // Take random angles THETA and PHI,
774  // and see if below Probability - if not - Redo
775  do {
776  thetaIndex = CLHEP::RandFlat::shootInt(thetaIndexMax-1);
777  phiIndex = CLHEP::RandFlat::shootInt(phiIndexMax-1);
778  // Find probability with the new indeces from LUT
779  AngularDistributionValue = OpticalSurface ->
780  GetAngularDistributionValue(angleIncident,
781  thetaIndex,
782  phiIndex);
783  } while ( !G4BooleanRand(AngularDistributionValue) );
784 
785  thetaRad = (-90 + 4*thetaIndex)*pi/180;
786  phiRad = (-90 + 5*phiIndex)*pi/180;
787  // Rotate Photon Momentum in Theta, then in Phi
788  NewMomentum = -OldMomentum;
789  PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
790  if (PerpendicularVectorTheta.mag() > kCarTolerance ) {
791  PerpendicularVectorPhi =
792  PerpendicularVectorTheta.cross(NewMomentum);
793  }
794  else {
795  PerpendicularVectorTheta = NewMomentum.orthogonal();
796  PerpendicularVectorPhi =
797  PerpendicularVectorTheta.cross(NewMomentum);
798  }
799  NewMomentum =
800  NewMomentum.rotate(anglePhotonToNormal-thetaRad,
801  PerpendicularVectorTheta);
802  NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
803  // Rotate Polarization too:
804  theFacetNormal = (NewMomentum - OldMomentum).unit();
805  EdotN = OldPolarization * theFacetNormal;
806  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
807  }
808  } while (NewMomentum * theGlobalNormal <= 0.0);
809 }
810 
811 void G4OpBoundaryProcess::DielectricDielectric()
812 {
813  G4bool Inside = false;
814  G4bool Swap = false;
815 
816  leap:
817 
818  G4bool Through = false;
819  G4bool Done = false;
820 
821  do {
822 
823  if (Through) {
824  Swap = !Swap;
825  Through = false;
826  theGlobalNormal = -theGlobalNormal;
827  G4SwapPtr(Material1,Material2);
828  G4SwapObj(&Rindex1,&Rindex2);
829  }
830 
831  if ( theFinish == polished ) {
832  theFacetNormal = theGlobalNormal;
833  }
834  else {
835  theFacetNormal =
836  GetFacetNormal(OldMomentum,theGlobalNormal);
837  }
838 
839  G4double PdotN = OldMomentum * theFacetNormal;
840  G4double EdotN = OldPolarization * theFacetNormal;
841 
842  cost1 = - PdotN;
843  if (std::abs(cost1) < 1.0-kCarTolerance){
844  sint1 = std::sqrt(1.-cost1*cost1);
845  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
846  }
847  else {
848  sint1 = 0.0;
849  sint2 = 0.0;
850  }
851 
852  if (sint2 >= 1.0) {
853 
854  // Simulate total internal reflection
855 
856  if (Swap) Swap = !Swap;
857 
858  theStatus = TotalInternalReflection;
859 
860  if ( theModel == unified && theFinish != polished )
861  ChooseReflection();
862 
863  if ( theStatus == LambertianReflection ) {
864  DoReflection();
865  }
866  else if ( theStatus == BackScattering ) {
867  NewMomentum = -OldMomentum;
868  NewPolarization = -OldPolarization;
869  }
870  else {
871 
872  PdotN = OldMomentum * theFacetNormal;
873  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
874  EdotN = OldPolarization * theFacetNormal;
875  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
876 
877  }
878  }
879  else if (sint2 < 1.0) {
880 
881  // Calculate amplitude for transmission (Q = P x N)
882 
883  if (cost1 > 0.0) {
884  cost2 = std::sqrt(1.-sint2*sint2);
885  }
886  else {
887  cost2 = -std::sqrt(1.-sint2*sint2);
888  }
889 
890  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
891  G4double E1_perp, E1_parl;
892 
893  if (sint1 > 0.0) {
894  A_trans = OldMomentum.cross(theFacetNormal);
895  A_trans = A_trans.unit();
896  E1_perp = OldPolarization * A_trans;
897  E1pp = E1_perp * A_trans;
898  E1pl = OldPolarization - E1pp;
899  E1_parl = E1pl.mag();
900  }
901  else {
902  A_trans = OldPolarization;
903  // Here we Follow Jackson's conventions and we set the
904  // parallel component = 1 in case of a ray perpendicular
905  // to the surface
906  E1_perp = 0.0;
907  E1_parl = 1.0;
908  }
909 
910  G4double s1 = Rindex1*cost1;
911  G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
912  G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
913  G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
914  G4double s2 = Rindex2*cost2*E2_total;
915 
916  G4double TransCoeff;
917 
918  if (theTransmittance > 0) TransCoeff = theTransmittance;
919  else if (cost1 != 0.0) TransCoeff = s2/s1;
920  else TransCoeff = 0.0;
921 
922  G4double E2_abs, C_parl, C_perp;
923 
924  if ( !G4BooleanRand(TransCoeff) ) {
925 
926  // Simulate reflection
927 
928  if (Swap) Swap = !Swap;
929 
930  theStatus = FresnelReflection;
931 
932  if ( theModel == unified && theFinish != polished )
933  ChooseReflection();
934 
935  if ( theStatus == LambertianReflection ) {
936  DoReflection();
937  }
938  else if ( theStatus == BackScattering ) {
939  NewMomentum = -OldMomentum;
940  NewPolarization = -OldPolarization;
941  }
942  else {
943 
944  PdotN = OldMomentum * theFacetNormal;
945  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
946 
947  if (sint1 > 0.0) { // incident ray oblique
948 
949  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
950  E2_perp = E2_perp - E1_perp;
951  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
952  A_paral = NewMomentum.cross(A_trans);
953  A_paral = A_paral.unit();
954  E2_abs = std::sqrt(E2_total);
955  C_parl = E2_parl/E2_abs;
956  C_perp = E2_perp/E2_abs;
957 
958  NewPolarization = C_parl*A_paral + C_perp*A_trans;
959 
960  }
961 
962  else { // incident ray perpendicular
963 
964  if (Rindex2 > Rindex1) {
965  NewPolarization = - OldPolarization;
966  }
967  else {
968  NewPolarization = OldPolarization;
969  }
970 
971  }
972  }
973  }
974  else { // photon gets transmitted
975 
976  // Simulate transmission/refraction
977 
978  Inside = !Inside;
979  Through = true;
980  theStatus = FresnelRefraction;
981 
982  if (sint1 > 0.0) { // incident ray oblique
983 
984  G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
985  NewMomentum = OldMomentum + alpha*theFacetNormal;
986  NewMomentum = NewMomentum.unit();
987  PdotN = -cost2;
988  A_paral = NewMomentum.cross(A_trans);
989  A_paral = A_paral.unit();
990  E2_abs = std::sqrt(E2_total);
991  C_parl = E2_parl/E2_abs;
992  C_perp = E2_perp/E2_abs;
993 
994  NewPolarization = C_parl*A_paral + C_perp*A_trans;
995 
996  }
997  else { // incident ray perpendicular
998 
999  NewMomentum = OldMomentum;
1000  NewPolarization = OldPolarization;
1001 
1002  }
1003  }
1004  }
1005 
1006  OldMomentum = NewMomentum.unit();
1007  OldPolarization = NewPolarization.unit();
1008 
1009  if (theStatus == FresnelRefraction) {
1010  Done = (NewMomentum * theGlobalNormal <= 0.0);
1011  }
1012  else {
1013  Done = (NewMomentum * theGlobalNormal >= 0.0);
1014  }
1015 
1016  } while (!Done);
1017 
1018  if (Inside && !Swap) {
1019  if( theFinish == polishedbackpainted ||
1020  theFinish == groundbackpainted ) {
1021 
1022  if( !G4BooleanRand(theReflectivity) ) {
1023  DoAbsorption();
1024  }
1025  else {
1026  if (theStatus != FresnelRefraction ) {
1027  theGlobalNormal = -theGlobalNormal;
1028  }
1029  else {
1030  Swap = !Swap;
1031  G4SwapPtr(Material1,Material2);
1032  G4SwapObj(&Rindex1,&Rindex2);
1033  }
1034  if ( theFinish == groundbackpainted )
1035  theStatus = LambertianReflection;
1036 
1037  DoReflection();
1038 
1039  theGlobalNormal = -theGlobalNormal;
1040  OldMomentum = NewMomentum;
1041 
1042  goto leap;
1043  }
1044  }
1045  }
1046 }
1047 
1048 // GetMeanFreePath
1049 // ---------------
1050 //
1051 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
1052  G4double ,
1053  G4ForceCondition* condition)
1054 {
1055  *condition = Forced;
1056 
1057  return DBL_MAX;
1058 }
1059 
1060 G4double G4OpBoundaryProcess::GetIncidentAngle()
1061 {
1062  G4double PdotN = OldMomentum * theFacetNormal;
1063  G4double magP= OldMomentum.mag();
1064  G4double magN= theFacetNormal.mag();
1065  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
1066 
1067  return incidentangle;
1068 }
1069 
1070 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1071  G4double E1_parl,
1072  G4double incidentangle,
1073  G4double RealRindex,
1074  G4double ImaginaryRindex)
1075 {
1076 
1077  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1078  G4complex N(RealRindex, ImaginaryRindex);
1079  G4complex CosPhi;
1080 
1081  G4complex u(1,0); //unit number 1
1082 
1083  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1084  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1085  G4complex denominatorTE, denominatorTM;
1086  G4complex rTM, rTE;
1087 
1088  // Following two equations, rTM and rTE, are from: "Introduction To Modern
1089  // Optics" written by Fowles
1090 
1091  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
1092 
1093  numeratorTE = std::cos(incidentangle) - N*CosPhi;
1094  denominatorTE = std::cos(incidentangle) + N*CosPhi;
1095  rTE = numeratorTE/denominatorTE;
1096 
1097  numeratorTM = N*std::cos(incidentangle) - CosPhi;
1098  denominatorTM = N*std::cos(incidentangle) + CosPhi;
1099  rTM = numeratorTM/denominatorTM;
1100 
1101  // This is my calculaton for reflectivity on a metalic surface
1102  // depending on the fraction of TE and TM polarization
1103  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1104  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1105 
1106  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1107  / (E1_perp*E1_perp + E1_parl*E1_parl);
1108  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1109  / (E1_perp*E1_perp + E1_parl*E1_parl);
1110  Reflectivity = Reflectivity_TE + Reflectivity_TM;
1111 
1112  do {
1113  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1114  {iTE = -1;}else{iTE = 1;}
1115  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1116  {iTM = -1;}else{iTM = 1;}
1117  } while(iTE<0&&iTM<0);
1118 
1119  return real(Reflectivity);
1120 
1121 }
1122 
1123 void G4OpBoundaryProcess::CalculateReflectivity()
1124 {
1125  G4double RealRindex =
1126  PropertyPointer1->Value(thePhotonMomentum);
1127  G4double ImaginaryRindex =
1128  PropertyPointer2->Value(thePhotonMomentum);
1129 
1130  // calculate FacetNormal
1131  if ( theFinish == ground ) {
1132  theFacetNormal =
1133  GetFacetNormal(OldMomentum, theGlobalNormal);
1134  } else {
1135  theFacetNormal = theGlobalNormal;
1136  }
1137 
1138  G4double PdotN = OldMomentum * theFacetNormal;
1139  cost1 = -PdotN;
1140 
1141  if (std::abs(cost1) < 1.0 - kCarTolerance) {
1142  sint1 = std::sqrt(1. - cost1*cost1);
1143  } else {
1144  sint1 = 0.0;
1145  }
1146 
1147  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1148  G4double E1_perp, E1_parl;
1149 
1150  if (sint1 > 0.0 ) {
1151  A_trans = OldMomentum.cross(theFacetNormal);
1152  A_trans = A_trans.unit();
1153  E1_perp = OldPolarization * A_trans;
1154  E1pp = E1_perp * A_trans;
1155  E1pl = OldPolarization - E1pp;
1156  E1_parl = E1pl.mag();
1157  }
1158  else {
1159  A_trans = OldPolarization;
1160  // Here we Follow Jackson's conventions and we set the
1161  // parallel component = 1 in case of a ray perpendicular
1162  // to the surface
1163  E1_perp = 0.0;
1164  E1_parl = 1.0;
1165  }
1166 
1167  //calculate incident angle
1168  G4double incidentangle = GetIncidentAngle();
1169 
1170  //calculate the reflectivity depending on incident angle,
1171  //polarization and complex refractive
1172 
1173  theReflectivity =
1174  GetReflectivity(E1_perp, E1_parl, incidentangle,
1175  RealRindex, ImaginaryRindex);
1176 }
Generated on Sat May 25 2013 14:34:11 for Geant4 by   doxygen 1.8.4