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

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

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