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

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

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