Changeset 1055 for trunk/source/processes/optical

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

Location:

trunk/source/processes/optical

Files:

• History (modified) (1 diff)
• include/G4OpBoundaryProcess.hh (modified) (5 diffs)
• src/G4OpBoundaryProcess.cc (modified) (7 diffs)

trunk/source/processes/optical/History

@@ -16 +16 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+23th Mar 2009 Peter Gumplinger (op-V09-02-01)
+              for complex index of refraction: (1)resample the reflectivity
+              every time in the do-while of DielectricMetal, but don't
+              re-calculate theFacetNormal (which has already been chosen
+              in CalculateReflectivity), also (2) avoid an infinite loop by
+              resetting iTE and iTM inside the do-while of GetReflectivity;
+              both are bug-fixes.
+
+14th Jan 2009 Peter Gumplinger (op-V09-02-00)
+              respond to bug report 1040 by fixing G4OpBoundaryProcess.cc 
 
 07th Nov 2008 Peter Gumplinger (op-V09-01-09)

trunk/source/processes/optical/include/G4OpBoundaryProcess.hh

@@ -25 +25 @@
 //
 //
-// $Id: G4OpBoundaryProcess.hh,v 1.18 2008/11/07 17:59:37 gum Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4OpBoundaryProcess.hh,v 1.19 2009/03/23 21:18:20 gum Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // 
@@ -54 +54 @@
 // mail:        gum@triumf.ca
 //
-// CVS version tag: 
 ////////////////////////////////////////////////////////////////////////
 
@@ -160 +159 @@
         // Returns the Reflectivity on a metalic surface
 
-        void           SetModel(G4OpticalSurfaceModel model);
+        void CalculateReflectivity(void);
+
+        void SetModel(G4OpticalSurfaceModel model);
         // Set the optical surface model to be followed
         // (glisur || unified).
@@ -195 +196 @@
 
         G4OpticalSurface* OpticalSurface;
+
+        G4MaterialPropertyVector* PropertyPointer;
+        G4MaterialPropertyVector* PropertyPointer1;
+        G4MaterialPropertyVector* PropertyPointer2;
 
         G4double Rindex1;
@@ -308 +313 @@
 
           theStatus = LobeReflection;
-          theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
+          if ( PropertyPointer1 && PropertyPointer2 ){
+          } else {
+             theFacetNormal =
+                 GetFacetNormal(OldMomentum,theGlobalNormal);
+          }
           G4double PdotN = OldMomentum * theFacetNormal;
           NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;

trunk/source/processes/optical/src/G4OpBoundaryProcess.cc

@@ -178 +178 @@
         else {
           G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
-               << " The Navigator reports that it returned an invalid normal"
-               << G4endl;
+                 << " The Navigator reports that it returned an invalid normal"
+                 << G4endl;
+          G4Exception("G4OpBoundaryProcess::PostStepDoIt",
+                      "Invalid Surface Normal",
+                      EventMustBeAborted,
+                      "Geometry must return valid surface normal");
         }
 
@@ -282 +286 @@
               }
 
-              G4MaterialPropertyVector* PropertyPointer;
-              G4MaterialPropertyVector* PropertyPointer1;
-              G4MaterialPropertyVector* PropertyPointer2;
-
               PropertyPointer =
                       aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
@@ -303 +303 @@
               } else if (PropertyPointer1 && PropertyPointer2) {
 
-                 G4double RealRindex =
-                          PropertyPointer1->GetProperty(thePhotonMomentum);
-                 G4double ImaginaryRindex =
-                          PropertyPointer2->GetProperty(thePhotonMomentum);
-
-                 // calculate FacetNormal
-                 if ( theFinish == ground ) {
-                    theFacetNormal =
-                              GetFacetNormal(OldMomentum, theGlobalNormal);
-                 } else {
-                    theFacetNormal = theGlobalNormal;
-                 }
-
-                 G4double PdotN = OldMomentum * theFacetNormal;
-                 cost1 = -PdotN;
-
-                 if (std::abs(cost1) < 1.0 - kCarTolerance) {
-                    sint1 = std::sqrt(1. - cost1*cost1);
-                 } else {
-                    sint1 = 0.0;
-                 }
-
-                 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
-                 G4double E1_perp, E1_parl;
-
-                 if (sint1 > 0.0 ) {
-                    A_trans = OldMomentum.cross(theFacetNormal);
-                    A_trans = A_trans.unit();
-                    E1_perp = OldPolarization * A_trans;
-                    E1pp    = E1_perp * A_trans;
-                    E1pl    = OldPolarization - E1pp;
-                    E1_parl = E1pl.mag();
-                 }
-                 else {
-                    A_trans  = OldPolarization;
-                    // Here we Follow Jackson's conventions and we set the
-                    // parallel component = 1 in case of a ray perpendicular
-                    // to the surface
-                    E1_perp  = 0.0;
-                    E1_parl  = 1.0;
-                 }
-
-                 //calculate incident angle
-                 G4double incidentangle = GetIncidentAngle();
-
-                 //calculate the reflectivity depending on incident angle,
-                 //polarization and complex refractive
-
-                 theReflectivity =
-                            GetReflectivity(E1_perp, E1_parl, incidentangle,
-                                                 RealRindex, ImaginaryRindex);
+                 CalculateReflectivity();
 
               } else {
@@ -591 +541 @@
            if( !G4BooleanRand(theReflectivity) && n == 1 ) {
 
-             // Comment out DoAbsorption if you wish to have
-             //       Transmission instead of Absorption
+             // Comment out DoAbsorption and uncomment theStatus = Absorption;
+             // if you wish to have Transmission instead of Absorption
 
              DoAbsorption();
-             theStatus = Absorption;
+             // theStatus = Absorption;
              break;
 
            }
            else {
+
+             if (PropertyPointer1 && PropertyPointer2) {
+                if ( n > 1 ) {
+                   CalculateReflectivity();
+                   if ( !G4BooleanRand(theReflectivity) ) {
+                      DoAbsorption();
+                      break;
+                   }
+                }
+             }
 
              if ( theModel == glisur || theFinish == polished ) {
@@ -618 +578 @@
                 else {
 
-                   if(theStatus==LobeReflection)theFacetNormal =
-                             GetFacetNormal(OldMomentum,theGlobalNormal);
+                   if(theStatus==LobeReflection){
+                     if ( PropertyPointer1 && PropertyPointer2 ){
+                     } else {
+                        theFacetNormal =
+                            GetFacetNormal(OldMomentum,theGlobalNormal);
+                     }
+                   }
 
                    G4double PdotN = OldMomentum * theFacetNormal;
@@ -962 +927 @@
 
   do {
-     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))iTE = -1;
-     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))iTM = -1;
+     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
+       {iTE = -1;}else{iTE = 1;}
+     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
+       {iTM = -1;}else{iTM = 1;}
   } while(iTE<0&&iTM<0);
 
@@ -969 +936 @@
 
 }
+
+void G4OpBoundaryProcess::CalculateReflectivity()
+{
+  G4double RealRindex =
+           PropertyPointer1->GetProperty(thePhotonMomentum);
+  G4double ImaginaryRindex =
+           PropertyPointer2->GetProperty(thePhotonMomentum);
+
+  // calculate FacetNormal
+  if ( theFinish == ground ) {
+     theFacetNormal =
+               GetFacetNormal(OldMomentum, theGlobalNormal);
+  } else {
+     theFacetNormal = theGlobalNormal;
+  }
+
+  G4double PdotN = OldMomentum * theFacetNormal;
+  cost1 = -PdotN;
+
+  if (std::abs(cost1) < 1.0 - kCarTolerance) {
+     sint1 = std::sqrt(1. - cost1*cost1);
+  } else {
+     sint1 = 0.0;
+  }
+
+  G4ThreeVector A_trans, A_paral, E1pp, E1pl;
+  G4double E1_perp, E1_parl;
+
+  if (sint1 > 0.0 ) {
+     A_trans = OldMomentum.cross(theFacetNormal);
+     A_trans = A_trans.unit();
+     E1_perp = OldPolarization * A_trans;
+     E1pp    = E1_perp * A_trans;
+     E1pl    = OldPolarization - E1pp;
+     E1_parl = E1pl.mag();
+  }
+  else {
+     A_trans  = OldPolarization;
+     // Here we Follow Jackson's conventions and we set the
+     // parallel component = 1 in case of a ray perpendicular
+     // to the surface
+     E1_perp  = 0.0;
+     E1_parl  = 1.0;
+  }
+
+  //calculate incident angle
+  G4double incidentangle = GetIncidentAngle();
+
+  //calculate the reflectivity depending on incident angle,
+  //polarization and complex refractive
+
+  theReflectivity =
+             GetReflectivity(E1_perp, E1_parl, incidentangle,
+                                                 RealRindex, ImaginaryRindex);
+}