source: trunk/source/processes/optical/src/G4OpBoundaryProcess.cc @ 1244

Last change on this file since 1244 was 1228, checked in by garnier, 14 years ago

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