FDsimG4OpBoundaryProcess.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: FDsimG4OpBoundaryProcess.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 //
62 // Author: Peter Gumplinger
63 // adopted from work by Werner Keil - April 2/96
64 // mail: gum@triumf.ca
65 //
67 
68 #include "G4ios.hh"
69 #include "FDsimG4OpBoundaryProcess.hh"
70 #include "G4GeometryTolerance.hh"
71 #include "ThinFilm.h"
72 
73 using namespace TelescopeSimulatorLX ;
74 
76 // Class Implementation
78 
80  // Operators
82 
83 // FDsimG4OpBoundaryProcess::operator=(const FDsimG4OpBoundaryProcess &right)
84 // {
85 // }
86 
88  // Constructors
90 
91 FDsimG4OpBoundaryProcess::FDsimG4OpBoundaryProcess(const G4String& processName,
92  G4ProcessType type)
93  : G4VDiscreteProcess(processName, type)
94 {
95  if ( verboseLevel > 0) {
96  G4cerr << GetProcessName() << " is created " << G4endl;
97  }
98 
99  theStatus = Undefined;
100  theModel = glisur;
101  theFinish = polished;
102  theReflectivity = 1.;
103  theEfficiency = 0.;
104 
105  prob_sl = 0.;
106  prob_ss = 0.;
107  prob_bs = 0.;
108 
109  kCarTolerance = G4GeometryTolerance::GetInstance()
110  ->GetSurfaceTolerance();
111 
112 }
113 
114 // FDsimG4OpBoundaryProcess::FDsimG4OpBoundaryProcess(const FDsimG4OpBoundaryProcess &right)
115 // {
116 // }
117 
119  // Destructors
121 
122 FDsimG4OpBoundaryProcess::~FDsimG4OpBoundaryProcess(){}
123 
125  // Methods
127 
128 // PostStepDoIt
129 // ------------
130 //
131 G4VParticleChange*
132 FDsimG4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
133 {
134  theStatus = Undefined;
135 
136  aParticleChange.Initialize(aTrack);
137 
138  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
139  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
140 
141  if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
142  theStatus = NotAtBoundary;
143  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
144  }
145  if (aTrack.GetStepLength()<=kCarTolerance/2){
146  theStatus = StepTooSmall;
147  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
148  }
149 
150  Material1 = pPreStepPoint -> GetMaterial();
151  Material2 = pPostStepPoint -> GetMaterial();
152 
153  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
154 
155  thePhotonMomentum = aParticle->GetTotalMomentum();
156  OldMomentum = aParticle->GetMomentumDirection();
157  OldPolarization = aParticle->GetPolarization();
158 
159  G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
160 
161  G4Navigator* theNavigator =
162  G4TransportationManager::GetTransportationManager()->
163  GetNavigatorForTracking();
164 
165  G4ThreeVector theLocalPoint = theNavigator->
166  GetGlobalToLocalTransform().
167  TransformPoint(theGlobalPoint);
168 
169  G4ThreeVector theLocalNormal; // Normal points back into volume
170 
171  G4bool valid;
172  theLocalNormal = theNavigator->GetLocalExitNormal(&valid);
173 
174  if (valid) {
175  theLocalNormal = -theLocalNormal;
176  }
177  else {
178  G4cerr << " FDsimG4OpBoundaryProcess/PostStepDoIt(): "
179  << " The Navigator reports that it returned an invalid normal"
180  << G4endl;
181  }
182 
183  theGlobalNormal = theNavigator->GetLocalToGlobalTransform().
184  TransformAxis(theLocalNormal);
185 
186  if (OldMomentum * theGlobalNormal > 0.0) {
187 #ifdef G4DEBUG_OPTICAL
188  G4cerr << " FDsimG4OpBoundaryProcess/PostStepDoIt(): "
189  << " theGlobalNormal points the wrong direction "
190  << G4endl;
191 #endif
192  theGlobalNormal = -theGlobalNormal;
193  }
194 
195  G4MaterialPropertiesTable* aMaterialPropertiesTable;
196  G4MaterialPropertyVector* Rindex;
197 
198  aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
199  if (aMaterialPropertiesTable) {
200  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
201  }
202  else {
203  theStatus = NoRINDEX;
204  aParticleChange.ProposeTrackStatus(fStopAndKill);
205  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
206  }
207 
208  if (Rindex) {
209  Rindex1 = Rindex->GetProperty(thePhotonMomentum);
210  }
211  else {
212  theStatus = NoRINDEX;
213  aParticleChange.ProposeTrackStatus(fStopAndKill);
214  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
215  }
216 
217  theModel = glisur;
218  theFinish = polished;
219 
220  G4SurfaceType type = dielectric_dielectric;
221 
222  Rindex = NULL;
223  OpticalSurface = NULL;
224 
225  G4LogicalSurface* Surface = G4LogicalBorderSurface::GetSurface
226  (pPreStepPoint ->GetPhysicalVolume(),
227  pPostStepPoint->GetPhysicalVolume());
228 
229  if (Surface == NULL){
230  G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
231  ->GetMotherLogical() ==
232  pPreStepPoint->GetPhysicalVolume()
233  ->GetLogicalVolume());
234  if(enteredDaughter){
235  Surface = G4LogicalSkinSurface::GetSurface
236  (pPostStepPoint->GetPhysicalVolume()->
237  GetLogicalVolume());
238  if(Surface == NULL)
239  Surface = G4LogicalSkinSurface::GetSurface
240  (pPreStepPoint->GetPhysicalVolume()->
241  GetLogicalVolume());
242  }
243  else {
244  Surface = G4LogicalSkinSurface::GetSurface
245  (pPreStepPoint->GetPhysicalVolume()->
246  GetLogicalVolume());
247  if(Surface == NULL)
248  Surface = G4LogicalSkinSurface::GetSurface
249  (pPostStepPoint->GetPhysicalVolume()->
250  GetLogicalVolume());
251  }
252  }
253 
254  // if (Surface) OpticalSurface = dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
255  if (Surface) OpticalSurface = (G4OpticalSurface*) Surface->GetSurfaceProperty();
256 
257  if (OpticalSurface) {
258 
259  type = OpticalSurface->GetType();
260  theModel = OpticalSurface->GetModel();
261  theFinish = OpticalSurface->GetFinish();
262 
263  aMaterialPropertiesTable = OpticalSurface->
264  GetMaterialPropertiesTable();
265 
266  if (aMaterialPropertiesTable) {
267 
268  if (theFinish == polishedbackpainted ||
269  theFinish == groundbackpainted ) {
270  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
271  if (Rindex) {
272  Rindex2 = Rindex->GetProperty(thePhotonMomentum);
273  }
274  else {
275  theStatus = NoRINDEX;
276  aParticleChange.ProposeTrackStatus(fStopAndKill);
277  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
278  }
279  }
280 
281  G4MaterialPropertyVector* PropertyPointer;
282  G4MaterialPropertyVector* PropertyPointer1;
283  G4MaterialPropertyVector* PropertyPointer2;
284 
285  PropertyPointer =
286  aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
287  PropertyPointer1 =
288  aMaterialPropertiesTable->GetProperty("REALRINDEX");
289  PropertyPointer2 =
290  aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
291 
292  iTE = 1;
293  iTM = 1;
294 
295  if (PropertyPointer) {
296 
297  theReflectivity =
298  PropertyPointer->GetProperty(thePhotonMomentum);
299 
300  } else if (PropertyPointer1 && PropertyPointer2) {
301 
302  G4double RealRindex =
303  PropertyPointer1->GetProperty(thePhotonMomentum);
304  G4double ImaginaryRindex =
305  PropertyPointer2->GetProperty(thePhotonMomentum);
306 
307  // calculate FacetNormal
308  if ( theFinish == ground ) {
309  theFacetNormal =
310  GetFacetNormal(OldMomentum, theGlobalNormal);
311  } else {
312  theFacetNormal = theGlobalNormal;
313  }
314 
315  G4double PdotN = OldMomentum * theFacetNormal;
316  cost1 = -PdotN;
317 
318  if (std::abs(cost1) < 1.0 - kCarTolerance) {
319  sint1 = std::sqrt(1. - cost1*cost1);
320  } else {
321  sint1 = 0.0;
322  }
323 
324  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
325  G4double E1_perp, E1_parl;
326 
327  if (sint1 > 0.0 ) {
328  A_trans = OldMomentum.cross(theFacetNormal);
329  A_trans = A_trans.unit();
330  E1_perp = OldPolarization * A_trans;
331  E1pp = E1_perp * A_trans;
332  E1pl = OldPolarization - E1pp;
333  E1_parl = E1pl.mag();
334  }
335  else {
336  A_trans = OldPolarization;
337  // Here we Follow Jackson's conventions and we set the
338  // parallel component = 1 in case of a ray perpendicular
339  // to the surface
340  E1_perp = 0.0;
341  E1_parl = 1.0;
342  }
343 
344  //calculate incident angle
345  G4double incidentangle = GetIncidentAngle();
346 
347  //calculate the reflectivity depending on incident angle,
348  //polarization and complex refractive
349 
350  // NEW CBT
351  G4double thickness = 0.0;
352  if (aMaterialPropertiesTable->ConstPropertyExists("THICKNESS"))
353  thickness = aMaterialPropertiesTable->GetConstProperty("THICKNESS");
354 
355  if (thickness != 0.0) {
356  G4complex RindexCmplx = G4complex(RealRindex, -ImaginaryRindex);
357  G4double Rindex3 = 1.0; // We assume vacuum;
358  theReflectivity =
359  GetReflectivityThinFilm(incidentangle,
360  Rindex1,RindexCmplx,Rindex3,
361  thickness);
362 
363 // G4cerr << "wvl,thickness,incidentangle,Rindex1,R*,reflectivity" <<
364 // h_Planck*c_light/thePhotonMomentum/nm << " " <<
365 // thickness/nm << " " <<
366 // incidentangle/deg << " " <<
367 // Rindex1 << " " <<
368 // RindexCmplx << " " <<
369 // theReflectivity << G4endl;
370  }
371  // END NEW CBT
372  else
373  theReflectivity =
374  GetReflectivity(E1_perp, E1_parl, incidentangle,
375  RealRindex, ImaginaryRindex);
376 
377 
378  } else {
379  theReflectivity = 1.0;
380  }
381 
382  PropertyPointer =
383  aMaterialPropertiesTable->GetProperty("EFFICIENCY");
384  if (PropertyPointer) {
385  theEfficiency =
386  PropertyPointer->GetProperty(thePhotonMomentum);
387  } else {
388  theEfficiency = 0.0;
389  }
390 
391  if ( theModel == unified ) {
392  PropertyPointer =
393  aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
394  if (PropertyPointer) {
395  prob_sl =
396  PropertyPointer->GetProperty(thePhotonMomentum);
397  } else {
398  prob_sl = 0.0;
399  }
400 
401  PropertyPointer =
402  aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
403  if (PropertyPointer) {
404  prob_ss =
405  PropertyPointer->GetProperty(thePhotonMomentum);
406  } else {
407  prob_ss = 0.0;
408  }
409 
410  PropertyPointer =
411  aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
412  if (PropertyPointer) {
413  prob_bs =
414  PropertyPointer->GetProperty(thePhotonMomentum);
415  } else {
416  prob_bs = 0.0;
417  }
418  }
419  }
420  else if (theFinish == polishedbackpainted ||
421  theFinish == groundbackpainted ) {
422  aParticleChange.ProposeTrackStatus(fStopAndKill);
423  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
424  }
425  }
426 
427  if (type == dielectric_dielectric ) {
428  if (theFinish == polished || theFinish == ground ) {
429 
430  if (Material1 == Material2){
431  theStatus = SameMaterial;
432  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
433  }
434  aMaterialPropertiesTable =
435  Material2->GetMaterialPropertiesTable();
436  if (aMaterialPropertiesTable)
437  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
438  if (Rindex) {
439  Rindex2 = Rindex->GetProperty(thePhotonMomentum);
440  }
441  else {
442  theStatus = NoRINDEX;
443  aParticleChange.ProposeTrackStatus(fStopAndKill);
444  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
445  }
446  }
447  }
448 
449  if ( verboseLevel > 0 ) {
450  G4cerr << " Photon at Boundary! " << G4endl;
451  G4cerr << " Old Momentum Direction: " << OldMomentum << G4endl;
452  G4cerr << " Old Polarization: " << OldPolarization << G4endl;
453  }
454 
455  if (type == dielectric_metal) {
456 
457  DielectricMetal();
458 
459  }
460  else if (type == dielectric_dielectric) {
461 
462  if ( theFinish == polishedfrontpainted ||
463  theFinish == groundfrontpainted ) {
464 
465  if( !G4BooleanRand(theReflectivity) ) {
466  DoAbsorption();
467  }
468  else {
469  if ( theFinish == groundfrontpainted )
470  theStatus = LambertianReflection;
471  DoReflection();
472  }
473  }
474  else {
475  DielectricDielectric();
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  G4cerr << " New Momentum Direction: " << NewMomentum << G4endl;
490  G4cerr << " New Polarization: " << NewPolarization << G4endl;
491  if ( theStatus == Undefined )
492  G4cerr << " *** Undefined *** " << G4endl;
493  if ( theStatus == FresnelRefraction )
494  G4cerr << " *** FresnelRefraction *** " << G4endl;
495  if ( theStatus == FresnelReflection )
496  G4cerr << " *** FresnelReflection *** " << G4endl;
497  if ( theStatus == TotalInternalReflection )
498  G4cerr << " *** TotalInternalReflection *** " << G4endl;
499  if ( theStatus == LambertianReflection )
500  G4cerr << " *** LambertianReflection *** " << G4endl;
501  if ( theStatus == LobeReflection )
502  G4cerr << " *** LobeReflection *** " << G4endl;
503  if ( theStatus == SpikeReflection )
504  G4cerr << " *** SpikeReflection *** " << G4endl;
505  if ( theStatus == BackScattering )
506  G4cerr << " *** BackScattering *** " << G4endl;
507  if ( theStatus == Absorption )
508  G4cerr << " *** Absorption *** " << G4endl;
509  if ( theStatus == Detection )
510  G4cerr << " *** Detection *** " << G4endl;
511  if ( theStatus == NotAtBoundary )
512  G4cerr << " *** NotAtBoundary *** " << G4endl;
513  if ( theStatus == SameMaterial )
514  G4cerr << " *** SameMaterial *** " << G4endl;
515  if ( theStatus == StepTooSmall )
516  G4cerr << " *** StepTooSmall *** " << G4endl;
517  if ( theStatus == NoRINDEX )
518  G4cerr << " *** NoRINDEX *** " << G4endl;
519  }
520 
521  aParticleChange.ProposeMomentumDirection(NewMomentum);
522  aParticleChange.ProposePolarization(NewPolarization);
523 
524  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
525 }
526 
527 G4ThreeVector
528 FDsimG4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
529  const G4ThreeVector& Normal ) const
530 {
531  G4ThreeVector FacetNormal;
532 
533  if (theModel == unified) {
534 
535  /* This function code alpha to a random value taken from the
536  distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
537  for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
538  is a gaussian distribution with mean 0 and standard deviation
539  sigma_alpha. */
540 
541  G4double alpha;
542 
543  G4double sigma_alpha = 0.0;
544  if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
545 
546  G4double f_max = std::min(1.0,4.*sigma_alpha);
547 
548  do {
549  do {
550  alpha = G4RandGauss::shoot(0.0,sigma_alpha);
551  } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
552 
553  G4double phi = G4UniformRand()*twopi;
554 
555  G4double SinAlpha = std::sin(alpha);
556  G4double CosAlpha = std::cos(alpha);
557  G4double SinPhi = std::sin(phi);
558  G4double CosPhi = std::cos(phi);
559 
560  G4double unit_x = SinAlpha * CosPhi;
561  G4double unit_y = SinAlpha * SinPhi;
562  G4double unit_z = CosAlpha;
563 
564  FacetNormal.setX(unit_x);
565  FacetNormal.setY(unit_y);
566  FacetNormal.setZ(unit_z);
567 
568  G4ThreeVector tmpNormal = Normal;
569 
570  FacetNormal.rotateUz(tmpNormal);
571  } while (Momentum * FacetNormal >= 0.0);
572  }
573  else {
574 
575  G4double polish = 1.0;
576  if (OpticalSurface) polish = OpticalSurface->GetPolish();
577 
578  if (polish < 1.0) {
579  do {
580  G4ThreeVector smear;
581  do {
582  smear.setX(2.*G4UniformRand()-1.0);
583  smear.setY(2.*G4UniformRand()-1.0);
584  smear.setZ(2.*G4UniformRand()-1.0);
585  } while (smear.mag()>1.0);
586  smear = (1.-polish) * smear;
587  FacetNormal = Normal + smear;
588  } while (Momentum * FacetNormal >= 0.0);
589  FacetNormal = FacetNormal.unit();
590  }
591  else {
592  FacetNormal = Normal;
593  }
594  }
595  return FacetNormal;
596 }
597 
598 void FDsimG4OpBoundaryProcess::DielectricMetal()
599 {
600  G4int n = 0;
601 
602  do {
603 
604  n++;
605 
606  if( !G4BooleanRand(theReflectivity) && n == 1 ) {
607 
608  DoAbsorption();
609  break;
610 
611  }
612  else {
613 
614  if ( theModel == glisur || theFinish == polished ) {
615 
616  DoReflection();
617 
618  } else {
619 
620  if ( n == 1 ) ChooseReflection();
621 
622  if ( theStatus == LambertianReflection ) {
623  DoReflection();
624  }
625  else if ( theStatus == BackScattering ) {
626  NewMomentum = -OldMomentum;
627  NewPolarization = -OldPolarization;
628  }
629  else {
630 
631  if(theStatus==LobeReflection)theFacetNormal =
632  GetFacetNormal(OldMomentum,theGlobalNormal);
633 
634  G4double PdotN = OldMomentum * theFacetNormal;
635  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
636  G4double EdotN = OldPolarization * theFacetNormal;
637 
638  G4ThreeVector A_trans, A_paral;
639 
640  if (sint1 > 0.0 ) {
641  A_trans = OldMomentum.cross(theFacetNormal);
642  A_trans = A_trans.unit();
643  } else {
644  A_trans = OldPolarization;
645  }
646  A_paral = NewMomentum.cross(A_trans);
647  A_paral = A_paral.unit();
648 
649  if(iTE>0&&iTM>0) {
650  NewPolarization =
651  -OldPolarization + (2.*EdotN)*theFacetNormal;
652  } else if (iTE>0) {
653  NewPolarization = -A_trans;
654  } else if (iTM>0) {
655  NewPolarization = -A_paral;
656  }
657 
658  }
659 
660  }
661 
662  OldMomentum = NewMomentum;
663  OldPolarization = NewPolarization;
664 
665  }
666 
667  } while (NewMomentum * theGlobalNormal < 0.0);
668 }
669 
670 void FDsimG4OpBoundaryProcess::DielectricDielectric()
671 {
672  G4bool Inside = false;
673  G4bool Swap = false;
674 
675  leap:
676 
677  G4bool Through = false;
678  G4bool Done = false;
679 
680  do {
681 
682  if (Through) {
683  Swap = !Swap;
684  Through = false;
685  theGlobalNormal = -theGlobalNormal;
686  G4Swap(Material1,Material2);
687  G4Swap(&Rindex1,&Rindex2);
688  }
689 
690  if ( theFinish == ground || theFinish == groundbackpainted ) {
691  theFacetNormal =
692  GetFacetNormal(OldMomentum,theGlobalNormal);
693  }
694  else {
695  theFacetNormal = theGlobalNormal;
696  }
697 
698  G4double PdotN = OldMomentum * theFacetNormal;
699  G4double EdotN = OldPolarization * theFacetNormal;
700 
701  cost1 = - PdotN;
702  if (std::abs(cost1) < 1.0-kCarTolerance){
703  sint1 = std::sqrt(1.-cost1*cost1);
704  sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
705  }
706  else {
707  sint1 = 0.0;
708  sint2 = 0.0;
709  }
710 
711  if (sint2 >= 1.0) {
712 
713  // Simulate total internal reflection
714 
715  if (Swap) Swap = !Swap;
716 
717  theStatus = TotalInternalReflection;
718 
719  if ( theModel == unified && theFinish != polished )
720  ChooseReflection();
721 
722  if ( theStatus == LambertianReflection ) {
723  DoReflection();
724  }
725  else if ( theStatus == BackScattering ) {
726  NewMomentum = -OldMomentum;
727  NewPolarization = -OldPolarization;
728  }
729  else {
730 
731  PdotN = OldMomentum * theFacetNormal;
732  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
733  EdotN = OldPolarization * theFacetNormal;
734  NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
735 
736  }
737  }
738  else if (sint2 < 1.0) {
739 
740  // Calculate amplitude for transmission (Q = P x N)
741 
742  if (cost1 > 0.0) {
743  cost2 = std::sqrt(1.-sint2*sint2);
744  }
745  else {
746  cost2 = -std::sqrt(1.-sint2*sint2);
747  }
748 
749  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
750  G4double E1_perp, E1_parl;
751 
752  if (sint1 > 0.0) {
753  A_trans = OldMomentum.cross(theFacetNormal);
754  A_trans = A_trans.unit();
755  E1_perp = OldPolarization * A_trans;
756  E1pp = E1_perp * A_trans;
757  E1pl = OldPolarization - E1pp;
758  E1_parl = E1pl.mag();
759  }
760  else {
761  A_trans = OldPolarization;
762  // Here we Follow Jackson's conventions and we set the
763  // parallel component = 1 in case of a ray perpendicular
764  // to the surface
765  E1_perp = 0.0;
766  E1_parl = 1.0;
767  }
768 
769  G4double s1 = Rindex1*cost1;
770  G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
771  G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
772  G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
773  G4double s2 = Rindex2*cost2*E2_total;
774 
775  G4double TransCoeff;
776 
777  if (cost1 != 0.0) {
778  TransCoeff = s2/s1;
779  }
780  else {
781  TransCoeff = 0.0;
782  }
783 
784  G4double E2_abs, C_parl, C_perp;
785 
786  if ( !G4BooleanRand(TransCoeff) ) {
787 
788  // Simulate reflection
789 
790  if (Swap) Swap = !Swap;
791 
792  theStatus = FresnelReflection;
793 
794  if ( theModel == unified && theFinish != polished )
795  ChooseReflection();
796 
797  if ( theStatus == LambertianReflection ) {
798  DoReflection();
799  }
800  else if ( theStatus == BackScattering ) {
801  NewMomentum = -OldMomentum;
802  NewPolarization = -OldPolarization;
803  }
804  else {
805 
806  PdotN = OldMomentum * theFacetNormal;
807  NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
808 
809  if (sint1 > 0.0) { // incident ray oblique
810 
811  E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
812  E2_perp = E2_perp - E1_perp;
813  E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
814  A_paral = NewMomentum.cross(A_trans);
815  A_paral = A_paral.unit();
816  E2_abs = std::sqrt(E2_total);
817  C_parl = E2_parl/E2_abs;
818  C_perp = E2_perp/E2_abs;
819 
820  NewPolarization = C_parl*A_paral + C_perp*A_trans;
821 
822  }
823 
824  else { // incident ray perpendicular
825 
826  if (Rindex2 > Rindex1) {
827  NewPolarization = - OldPolarization;
828  }
829  else {
830  NewPolarization = OldPolarization;
831  }
832 
833  }
834  }
835  }
836  else { // photon gets transmitted
837 
838  // Simulate transmission/refraction
839 
840  Inside = !Inside;
841  Through = true;
842  theStatus = FresnelRefraction;
843 
844  if (sint1 > 0.0) { // incident ray oblique
845 
846  G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
847  NewMomentum = OldMomentum + alpha*theFacetNormal;
848  NewMomentum = NewMomentum.unit();
849  PdotN = -cost2;
850  A_paral = NewMomentum.cross(A_trans);
851  A_paral = A_paral.unit();
852  E2_abs = std::sqrt(E2_total);
853  C_parl = E2_parl/E2_abs;
854  C_perp = E2_perp/E2_abs;
855 
856  NewPolarization = C_parl*A_paral + C_perp*A_trans;
857 
858  }
859  else { // incident ray perpendicular
860 
861  NewMomentum = OldMomentum;
862  NewPolarization = OldPolarization;
863 
864  }
865  }
866  }
867 
868  OldMomentum = NewMomentum.unit();
869  OldPolarization = NewPolarization.unit();
870 
871  if (theStatus == FresnelRefraction) {
872  Done = (NewMomentum * theGlobalNormal <= 0.0);
873  }
874  else {
875  Done = (NewMomentum * theGlobalNormal >= 0.0);
876  }
877 
878  } while (!Done);
879 
880  if (Inside && !Swap) {
881  if( theFinish == polishedbackpainted ||
882  theFinish == groundbackpainted ) {
883 
884  if( !G4BooleanRand(theReflectivity) ) {
885  DoAbsorption();
886  }
887  else {
888  if (theStatus != FresnelRefraction ) {
889  theGlobalNormal = -theGlobalNormal;
890  }
891  else {
892  Swap = !Swap;
893  G4Swap(Material1,Material2);
894  G4Swap(&Rindex1,&Rindex2);
895  }
896  if ( theFinish == groundbackpainted )
897  theStatus = LambertianReflection;
898 
899  DoReflection();
900 
901  theGlobalNormal = -theGlobalNormal;
902  OldMomentum = NewMomentum;
903 
904  goto leap;
905  }
906  }
907  }
908 }
909 
910 // GetMeanFreePath
911 // ---------------
912 //
913 G4double FDsimG4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
914  G4double ,
915  G4ForceCondition* condition)
916 {
917  *condition = Forced;
918 
919  return DBL_MAX;
920 }
921 
922 G4double FDsimG4OpBoundaryProcess::GetIncidentAngle()
923 {
924  G4double PdotN = OldMomentum * theFacetNormal;
925  G4double magP= OldMomentum.mag();
926  G4double magN= theFacetNormal.mag();
927  G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
928 
929  return incidentangle;
930 }
931 
932 G4double FDsimG4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
933  G4double E1_parl,
934  G4double incidentangle,
935  G4double RealRindex,
936  G4double ImaginaryRindex)
937 {
938 
939  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
940  G4complex N(RealRindex, ImaginaryRindex);
941  G4complex CosPhi;
942 
943  G4complex u(1,0); //unit number 1
944 
945  G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
946  G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
947  G4complex denominatorTE, denominatorTM;
948  G4complex rTM, rTE;
949 
950  // Following two equations, rTM and rTE, are from: "Introduction To Modern
951  // Optics" written by Fowles
952 
953  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
954 
955  numeratorTE = std::cos(incidentangle) - N*CosPhi;
956  denominatorTE = std::cos(incidentangle) + N*CosPhi;
957  rTE = numeratorTE/denominatorTE;
958 
959  numeratorTM = N*std::cos(incidentangle) - CosPhi;
960  denominatorTM = N*std::cos(incidentangle) + CosPhi;
961  rTM = numeratorTM/denominatorTM;
962 
963  // This is my calculaton for reflectivity on a metalic surface
964  // depending on the fraction of TE and TM polarization
965  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
966  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
967 
968  Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
969  / (E1_perp*E1_perp + E1_parl*E1_parl);
970  Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
971  / (E1_perp*E1_perp + E1_parl*E1_parl);
972  Reflectivity = Reflectivity_TE + Reflectivity_TM;
973 
974 // do {
975 // if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))iTE = -1;
976 // if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))iTM = -1;
977 // } while(iTE<0&&iTM<0);
978 
979 // CBT The loop above was giving an infinite loop. Replaced by the following from Geant4.9.2
980 // See entry 13 in Geant4 forum on Processes involving optical photons (http://hypernews.slac.stanford.edu/HyperNews/geant4/get/opticalphotons.html)
981  do {
982  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
983  {iTE = -1;}else{iTE = 1;}
984  if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
985  {iTM = -1;}else{iTM = 1;}
986  } while(iTE<0&&iTM<0);
987 
988  return real(Reflectivity);
989 
990 }
991 
992 G4double
993 FDsimG4OpBoundaryProcess::GetReflectivityThinFilm(G4double incidentangle,
994  G4double n1,G4complex n2, G4double n3,
995  G4double thickness)
996 {
997  G4double sinth = sin(incidentangle);
998  G4double lambda = h_Planck*c_light/thePhotonMomentum;
999  G4double Reflectivity = ThinFilm::GetReflectance(n1,n2,n3,sinth,thickness/nm, lambda/nm);
1000 
1001  return Reflectivity;
1002 }
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::GetReflectivity
G4double GetReflectivity(G4double E1_perp, G4double E1_parl, G4double incidentangle, G4double RealRindex, G4double ImaginaryRindex)
Definition: FDsimG4OpBoundaryProcess.cc:932
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::GetFacetNormal
G4ThreeVector GetFacetNormal(const G4ThreeVector &Momentum, const G4ThreeVector &Normal) const
Definition: FDsimG4OpBoundaryProcess.cc:528
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::GetMeanFreePath
G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *condition)
Definition: FDsimG4OpBoundaryProcess.cc:913
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::GetReflectivityThinFilm
G4double GetReflectivityThinFilm(G4double incidentangle, G4double n1, G4complex n2, G4double n3, G4double thickness)
Definition: FDsimG4OpBoundaryProcess.cc:993
utl::abs
double abs(const SVector< n, T > &v)
Definition: ValidationTests/SdCalibration/Validatrix.cc:32
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::FDsimG4OpBoundaryProcess
FDsimG4OpBoundaryProcess(const G4String &processName="OpBoundary", G4ProcessType type=fOptical)
Definition: FDsimG4OpBoundaryProcess.cc:91
TelescopeSimulatorLX::FDsimG4OpBoundaryProcess::PostStepDoIt
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: FDsimG4OpBoundaryProcess.cc:132
TelescopeSimulatorLX::ThinFilm::GetReflectance
static double GetReflectance(const double nIndex1, const dcomplex &nIndex2, const double nIndex3, const double sinTheta1, const double thickness, const double wavelength)
Definition: ThinFilm.h:42
, generated on Tue Sep 26 2023.