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