/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Private Member Functions | Private Attributes
DsG4OpBoundaryProcess Class Reference

#include <DsG4OpBoundaryProcess.h>

List of all members.

Public Member Functions

 DsG4OpBoundaryProcess (const G4String &processName="OpBoundary", G4ProcessType type=fOptical)
 ~DsG4OpBoundaryProcess ()
G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
G4double GetMeanFreePath (const G4Track &, G4double, G4ForceCondition *condition)
G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
G4OpticalSurfaceModel GetModel () const
DsG4OpBoundaryProcessStatus GetStatus () const
G4double GetIncidentAngle ()
G4double GetReflectivity (G4double E1_perp, G4double E1_parl, G4double incidentangle, G4double RealRindex, G4double ImaginaryRindex)
void SetModel (G4OpticalSurfaceModel model)

Private Member Functions

G4bool G4BooleanRand (const G4double prob) const
G4ThreeVector GetFacetNormal (const G4ThreeVector &Momentum, const G4ThreeVector &Normal) const
void DielectricMetal ()
void DielectricDielectric ()
void ChooseReflection ()
void DoAbsorption ()
void DoReflection ()

Private Attributes

G4double thePhotonMomentum
G4ThreeVector OldMomentum
G4ThreeVector OldPolarization
G4ThreeVector NewMomentum
G4ThreeVector NewPolarization
G4ThreeVector theGlobalNormal
G4ThreeVector theFacetNormal
G4Material * Material1
G4Material * Material2
G4OpticalSurface * OpticalSurface
G4double Rindex1
G4double Rindex2
G4double cost1
G4double cost2
G4double sint1
G4double sint2
DsG4OpBoundaryProcessStatus theStatus
G4OpticalSurfaceModel theModel
G4OpticalSurfaceFinish theFinish
G4double theReflectivity
G4double theEfficiency
G4double prob_sl
G4double prob_ss
G4double prob_bs
G4int iTE
G4int iTM
G4double kCarTolerance
G4int abNormalCounter

Detailed Description

Definition at line 101 of file DsG4OpBoundaryProcess.h.


Constructor & Destructor Documentation

DsG4OpBoundaryProcess::DsG4OpBoundaryProcess ( const G4String &  processName = "OpBoundary",
G4ProcessType  type = fOptical 
)

Definition at line 90 of file DsG4OpBoundaryProcess.cc.

             : G4VDiscreteProcess(processName, type)
{
        if ( verboseLevel > 0) {
           G4cout << GetProcessName() << " is created " << G4endl;
        }

        SetProcessSubType(fOpBoundary);

        theStatus = Undefined;
        theModel = glisur;
        theFinish = polished;
        theReflectivity = 1.;
        theEfficiency   = 0.;

        prob_sl = 0.;
        prob_ss = 0.;
        prob_bs = 0.;

        kCarTolerance = G4GeometryTolerance::GetInstance()
                        ->GetSurfaceTolerance();

        abNormalCounter = 0;

}
DsG4OpBoundaryProcess::~DsG4OpBoundaryProcess ( )

Definition at line 125 of file DsG4OpBoundaryProcess.cc.

{}

Member Function Documentation

G4bool DsG4OpBoundaryProcess::IsApplicable ( const G4ParticleDefinition &  aParticleType) [inline]

Definition at line 233 of file DsG4OpBoundaryProcess.h.

{
   return ( &aParticleType == G4OpticalPhoton::OpticalPhoton() );
}
G4double DsG4OpBoundaryProcess::GetMeanFreePath ( const G4Track &  ,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 937 of file DsG4OpBoundaryProcess.cc.

{
        *condition = Forced;

        return DBL_MAX;
}
G4VParticleChange * DsG4OpBoundaryProcess::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)

Definition at line 135 of file DsG4OpBoundaryProcess.cc.

{
        theStatus = Undefined;

        aParticleChange.Initialize(aTrack);

        G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
        G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();

        if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
                theStatus = NotAtBoundary;
                return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }
        if (aTrack.GetStepLength()<=kCarTolerance/2){
                theStatus = StepTooSmall;
                return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }

        Material1 = pPreStepPoint  -> GetMaterial();
        Material2 = pPostStepPoint -> GetMaterial();

        const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();

        thePhotonMomentum = aParticle->GetTotalMomentum();
        OldMomentum       = aParticle->GetMomentumDirection();
        OldPolarization   = aParticle->GetPolarization();

        G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();

        G4Navigator* theNavigator =
                     G4TransportationManager::GetTransportationManager()->
                                              GetNavigatorForTracking();

        G4ThreeVector theLocalPoint = theNavigator->
                                      GetGlobalToLocalTransform().
                                      TransformPoint(theGlobalPoint);

        G4ThreeVector theLocalNormal;   // Normal points back into volume

        G4bool valid;
        theLocalNormal = theNavigator->GetLocalExitNormal(&valid);

        if (valid) {
          theLocalNormal = -theLocalNormal;
        }
        else {
          G4cerr << " DsG4OpBoundaryProcess/PostStepDoIt(): "
               << " The Navigator reports that it returned an invalid normal"
               << G4endl;
        }

        theGlobalNormal = theNavigator->GetLocalToGlobalTransform().
                                        TransformAxis(theLocalNormal);

        if(theGlobalNormal.mag() == 0) {
          abNormalCounter++;
          std::cout << "Because of normal = 0, the number of the killed optical photons is " << abNormalCounter << std::endl;
          aParticleChange.ProposeTrackStatus(fStopAndKill);
          return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }

        if (OldMomentum * theGlobalNormal > 0.0) {
#ifdef G4DEBUG_OPTICAL
           G4cerr << " DsG4OpBoundaryProcess/PostStepDoIt(): "
                  << " theGlobalNormal points the wrong direction "
                  << G4endl;
#endif
           theGlobalNormal = -theGlobalNormal;
        }

        G4MaterialPropertiesTable* aMaterialPropertiesTable;
        G4MaterialPropertyVector* Rindex;

        aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
        if (aMaterialPropertiesTable) {
                Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
        }
        else {
                theStatus = NoRINDEX;
                aParticleChange.ProposeTrackStatus(fStopAndKill);
                return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }

        if (Rindex) {
                Rindex1 = Rindex->GetProperty(thePhotonMomentum);
        }
        else {
                theStatus = NoRINDEX;
                aParticleChange.ProposeTrackStatus(fStopAndKill);
                return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
        }

        theModel = glisur;
        theFinish = polished;

        G4SurfaceType type = dielectric_dielectric;

        Rindex = NULL;
        OpticalSurface = NULL;

        G4LogicalSurface* Surface = NULL;

        Surface = G4LogicalBorderSurface::GetSurface
                  (pPreStepPoint ->GetPhysicalVolume(),
                   pPostStepPoint->GetPhysicalVolume());

        if (Surface == NULL){
          G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
                                  ->GetMotherLogical() ==
                                  pPreStepPoint->GetPhysicalVolume()
                                  ->GetLogicalVolume());
          if(enteredDaughter){
            Surface = G4LogicalSkinSurface::GetSurface
              (pPostStepPoint->GetPhysicalVolume()->
               GetLogicalVolume());
            if(Surface == NULL)
              Surface = G4LogicalSkinSurface::GetSurface
              (pPreStepPoint->GetPhysicalVolume()->
               GetLogicalVolume());
          }
          else {
            Surface = G4LogicalSkinSurface::GetSurface
              (pPreStepPoint->GetPhysicalVolume()->
               GetLogicalVolume());
            if(Surface == NULL)
              Surface = G4LogicalSkinSurface::GetSurface
              (pPostStepPoint->GetPhysicalVolume()->
               GetLogicalVolume());
          }
        }

        if (Surface) OpticalSurface = 
           dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());

        if (OpticalSurface) {

           type      = OpticalSurface->GetType();
           theModel  = OpticalSurface->GetModel();
           theFinish = OpticalSurface->GetFinish();

           aMaterialPropertiesTable = OpticalSurface->
                                        GetMaterialPropertiesTable();

           if (aMaterialPropertiesTable) {

              if (theFinish == polishedbackpainted ||
                  theFinish == groundbackpainted ) {
                  Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
                  if (Rindex) {
                     Rindex2 = Rindex->GetProperty(thePhotonMomentum);
                  }
                  else {
                     theStatus = NoRINDEX;
                     aParticleChange.ProposeTrackStatus(fStopAndKill);
                     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                  }
              }

              G4MaterialPropertyVector* PropertyPointer;
              G4MaterialPropertyVector* PropertyPointer1;
              G4MaterialPropertyVector* PropertyPointer2;

              PropertyPointer =
                      aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
              PropertyPointer1 =
                      aMaterialPropertiesTable->GetProperty("REALRINDEX");
              PropertyPointer2 =
                      aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");

              iTE = 1;
              iTM = 1;

              if (PropertyPointer) {

                 theReflectivity =
                          PropertyPointer->GetProperty(thePhotonMomentum);
                 if(OpticalSurface->GetName().contains("ESRAir")) {
                   G4double inciAngle = GetIncidentAngle();
                   //ESR in air
                   if(inciAngle*180./pi > 40) {
                     theReflectivity = (theReflectivity - 0.993) + 0.973572 + 9.53233e-04*(inciAngle*180./pi) - 1.22184e-05*((inciAngle*180./pi))*((inciAngle*180./pi));
                   }
                   //ESR in oil
                   //if(inciAngle*180./pi > 40 && inciAngle*180./pi <= 50) {
                   //  theReflectivity = (theReflectivity - 0.993) + (inciAngle*180./pi - 40)*0.99/10. + (50 - inciAngle*180./pi)*0.993/10.;
                   //}
                   //else if (inciAngle*180./pi > 50 && inciAngle*180./pi <= 60) {
                   //  theReflectivity = (theReflectivity - 0.993) + (inciAngle*180./pi - 50)*0.94/10. + (60 - inciAngle*180./pi)*0.99/10.;
                   //}
                   //else if (inciAngle*180./pi > 60 && inciAngle*180./pi <= 65) {
                   //  theReflectivity = (theReflectivity - 0.993) + (inciAngle*180./pi - 60)*0.62/5. + (65 - inciAngle*180./pi)*0.94/5.;
                   //}
                   //else if (inciAngle*180./pi > 65 && inciAngle*180./pi <= 80) {
                   //  theReflectivity = (theReflectivity - 0.993) + 0.62;
                   //}
                   //else if (inciAngle*180./pi > 80 && inciAngle*180./pi <= 85) {
                   //  theReflectivity = (theReflectivity - 0.993) + (inciAngle*180./pi - 80)*0.72/5. + (85 - inciAngle*180./pi)*0.62/5.;
                   //}
                   //else if (inciAngle*180./pi > 85 && inciAngle*180./pi < 90) {
                   //  theReflectivity = (theReflectivity - 0.993) + 0.72;
                   //}
                 }

              } 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);

              } else {
                 theReflectivity = 1.0;
              }

              PropertyPointer =
              aMaterialPropertiesTable->GetProperty("EFFICIENCY");
              if (PropertyPointer) {
                      theEfficiency =
                      PropertyPointer->GetProperty(thePhotonMomentum);
              } else {
                      theEfficiency = 0.0;
              }

              if ( theModel == unified ) {
                PropertyPointer =
                aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
                if (PropertyPointer) {
                         prob_sl =
                         PropertyPointer->GetProperty(thePhotonMomentum);
                } else {
                         prob_sl = 0.0;
                }

                PropertyPointer =
                aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
                if (PropertyPointer) {
                         prob_ss =
                         PropertyPointer->GetProperty(thePhotonMomentum);
                } else {
                         prob_ss = 0.0;
                }

                PropertyPointer =
                aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
                if (PropertyPointer) {
                         prob_bs =
                         PropertyPointer->GetProperty(thePhotonMomentum);
                } else {
                         prob_bs = 0.0;
                }
              }
           }
           else if (theFinish == polishedbackpainted ||
                    theFinish == groundbackpainted ) {
                      aParticleChange.ProposeTrackStatus(fStopAndKill);
                      return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
           }
        }

        if (type == dielectric_dielectric ) {
           if (theFinish == polished || theFinish == ground ) {

              if (Material1 == Material2){
                 theStatus = SameMaterial;
                 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
              }
              aMaterialPropertiesTable =
                     Material2->GetMaterialPropertiesTable();
              if (aMaterialPropertiesTable)
                 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
              if (Rindex) {
                 Rindex2 = Rindex->GetProperty(thePhotonMomentum);
              }
              else {
                 theStatus = NoRINDEX;
                 aParticleChange.ProposeTrackStatus(fStopAndKill);
                 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
              }
           }
        }

        if ( verboseLevel > 0 ) {
                G4cout << " Photon at Boundary! " << G4endl;
                G4cout << " Old Momentum Direction: " << OldMomentum     << G4endl;
                G4cout << " Old Polarization:       " << OldPolarization << G4endl;
        }

        if (type == dielectric_metal) {

          DielectricMetal();

          // Uncomment the following lines if you wish to have 
          //         Transmission instead of Absorption
          // if (theStatus == Absorption) {
          //    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
          // }

        }
        else if (type == dielectric_dielectric) {

          if ( theFinish == polishedfrontpainted ||
               theFinish == groundfrontpainted ) {

                  if( !G4BooleanRand(theReflectivity) ) {
                    DoAbsorption();
                  }
                  else {
                    if ( theFinish == groundfrontpainted )
                                        theStatus = LambertianReflection;
                    DoReflection();
                  }
          }
          else {
                  DielectricDielectric();
          }
        }
        else {

          G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
          return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);

        }

        NewMomentum = NewMomentum.unit();
        NewPolarization = NewPolarization.unit();

        if ( verboseLevel > 0) {
                G4cout << " New Momentum Direction: " << NewMomentum     << G4endl;
                G4cout << " New Polarization:       " << NewPolarization << G4endl;
                if ( theStatus == Undefined )
                        G4cout << " *** Undefined *** " << G4endl;
                if ( theStatus == FresnelRefraction )
                        G4cout << " *** FresnelRefraction *** " << G4endl;
                if ( theStatus == FresnelReflection )
                        G4cout << " *** FresnelReflection *** " << G4endl;
                if ( theStatus == TotalInternalReflection )
                        G4cout << " *** TotalInternalReflection *** " << G4endl;
                if ( theStatus == LambertianReflection )
                        G4cout << " *** LambertianReflection *** " << G4endl;
                if ( theStatus == LobeReflection )
                        G4cout << " *** LobeReflection *** " << G4endl;
                if ( theStatus == SpikeReflection )
                        G4cout << " *** SpikeReflection *** " << G4endl;
                if ( theStatus == BackScattering )
                        G4cout << " *** BackScattering *** " << G4endl;
                if ( theStatus == Absorption )
                        G4cout << " *** Absorption *** " << G4endl;
                if ( theStatus == Detection )
                        G4cout << " *** Detection *** " << G4endl;
                if ( theStatus == NotAtBoundary )
                        G4cout << " *** NotAtBoundary *** " << G4endl;
                if ( theStatus == SameMaterial )
                        G4cout << " *** SameMaterial *** " << G4endl;
                if ( theStatus == StepTooSmall )
                        G4cout << " *** StepTooSmall *** " << G4endl;
                if ( theStatus == NoRINDEX )
                        G4cout << " *** NoRINDEX *** " << G4endl;
        }

        aParticleChange.ProposeMomentumDirection(NewMomentum);
        aParticleChange.ProposePolarization(NewPolarization);

        return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
G4OpticalSurfaceModel DsG4OpBoundaryProcess::GetModel ( ) const [inline]

Definition at line 240 of file DsG4OpBoundaryProcess.h.

{
   return theModel;
}
DsG4OpBoundaryProcessStatus DsG4OpBoundaryProcess::GetStatus ( ) const [inline]

Definition at line 246 of file DsG4OpBoundaryProcess.h.

{
   return theStatus;
}
G4double DsG4OpBoundaryProcess::GetIncidentAngle ( )

Definition at line 946 of file DsG4OpBoundaryProcess.cc.

{
    G4double PdotN = OldMomentum * theFacetNormal;
    G4double magP= OldMomentum.mag();
    G4double magN= theFacetNormal.mag();
    G4double incidentangle = pi - std::acos(PdotN/(magP*magN));

    return incidentangle;
}
G4double DsG4OpBoundaryProcess::GetReflectivity ( G4double  E1_perp,
G4double  E1_parl,
G4double  incidentangle,
G4double  RealRindex,
G4double  ImaginaryRindex 
)

Definition at line 956 of file DsG4OpBoundaryProcess.cc.

{

  G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
  G4complex N(RealRindex, ImaginaryRindex);
  G4complex CosPhi;

  G4complex u(1,0);           //unit number 1

  G4complex numeratorTE;      // E1_perp=1 E1_parl=0 -> TE polarization
  G4complex numeratorTM;      // E1_parl=1 E1_perp=0 -> TM polarization
  G4complex denominatorTE, denominatorTM;
  G4complex rTM, rTE;

  // Following two equations, rTM and rTE, are from: "Introduction To Modern
  // Optics" written by Fowles

  CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));

  numeratorTE   = std::cos(incidentangle) - N*CosPhi;
  denominatorTE = std::cos(incidentangle) + N*CosPhi;
  rTE = numeratorTE/denominatorTE;

  numeratorTM   = N*std::cos(incidentangle) - CosPhi;
  denominatorTM = N*std::cos(incidentangle) + CosPhi;
  rTM = numeratorTM/denominatorTM;

  // This is my calculaton for reflectivity on a metalic surface
  // depending on the fraction of TE and TM polarization
  // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
  // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2

  Reflectivity_TE =  (rTE*conj(rTE))*(E1_perp*E1_perp)
                    / (E1_perp*E1_perp + E1_parl*E1_parl);
  Reflectivity_TM =  (rTM*conj(rTM))*(E1_parl*E1_parl)
                    / (E1_perp*E1_perp + E1_parl*E1_parl);
  Reflectivity    = Reflectivity_TE + Reflectivity_TM;

  do {
     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))iTE = -1;
     if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))iTM = -1;
  } while(iTE<0&&iTM<0);

  return real(Reflectivity);

}
void DsG4OpBoundaryProcess::SetModel ( G4OpticalSurfaceModel  model) [inline]

Definition at line 252 of file DsG4OpBoundaryProcess.h.

{
   theModel = model;
}
G4bool DsG4OpBoundaryProcess::G4BooleanRand ( const G4double  prob) const [inline, private]

Definition at line 225 of file DsG4OpBoundaryProcess.h.

{
  /* Returns a random boolean variable with the specified probability */

  return (G4UniformRand() < prob);
}
G4ThreeVector DsG4OpBoundaryProcess::GetFacetNormal ( const G4ThreeVector &  Momentum,
const G4ThreeVector &  Normal 
) const [private]

Definition at line 548 of file DsG4OpBoundaryProcess.cc.

{
        G4ThreeVector FacetNormal;

        if (theModel == unified) {

        /* This function code alpha to a random value taken from the
           distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
           for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
           is a gaussian distribution with mean 0 and standard deviation
           sigma_alpha.  */

           G4double alpha;

           G4double sigma_alpha = 0.0;
           if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();

           G4double f_max = std::min(1.0,4.*sigma_alpha);

           do {
              do {
                 alpha = G4RandGauss::shoot(0.0,sigma_alpha);
              } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );

              G4double phi = G4UniformRand()*twopi;

              G4double SinAlpha = std::sin(alpha);
              G4double CosAlpha = std::cos(alpha);
              G4double SinPhi = std::sin(phi);
              G4double CosPhi = std::cos(phi);

              G4double unit_x = SinAlpha * CosPhi;
              G4double unit_y = SinAlpha * SinPhi;
              G4double unit_z = CosAlpha;

              FacetNormal.setX(unit_x);
              FacetNormal.setY(unit_y);
              FacetNormal.setZ(unit_z);

              G4ThreeVector tmpNormal = Normal;

              FacetNormal.rotateUz(tmpNormal);
           } while (Momentum * FacetNormal >= 0.0);
        }
        else {

           G4double  polish = 1.0;
           if (OpticalSurface) polish = OpticalSurface->GetPolish();

           if (polish < 1.0) {
              do {
                 G4ThreeVector smear;
                 do {
                    smear.setX(2.*G4UniformRand()-1.0);
                    smear.setY(2.*G4UniformRand()-1.0);
                    smear.setZ(2.*G4UniformRand()-1.0);
                 } while (smear.mag()>1.0);
                 smear = (1.-polish) * smear;
                 FacetNormal = Normal + smear;
              } while (Momentum * FacetNormal >= 0.0);
              FacetNormal = FacetNormal.unit();
           }
           else {
              FacetNormal = Normal;
           }
        }
        return FacetNormal;
}
void DsG4OpBoundaryProcess::DielectricMetal ( ) [private]

Definition at line 618 of file DsG4OpBoundaryProcess.cc.

{
        G4int n = 0;

        do {

           n++;

           if( !G4BooleanRand(theReflectivity) && n == 1 ) {

             // Comment out DoAbsorption and uncomment theStatus = Absorption;
             // if you wish to have Transmission instead of Absorption

             DoAbsorption();
             // theStatus = Absorption;
             break;

           }
           else {

             if ( theModel == glisur || theFinish == polished ) {

                DoReflection();

             } else {

                if ( n == 1 ) ChooseReflection();
                                                                                
                if ( theStatus == LambertianReflection ) {
                   DoReflection();
                }
                else if ( theStatus == BackScattering ) {
                   NewMomentum = -OldMomentum;
                   NewPolarization = -OldPolarization;
                }
                else {

                   if(theStatus==LobeReflection)theFacetNormal =
                             GetFacetNormal(OldMomentum,theGlobalNormal);

                   G4double PdotN = OldMomentum * theFacetNormal;
                   NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
                   G4double EdotN = OldPolarization * theFacetNormal;

                   G4ThreeVector A_trans, A_paral;

                   if (sint1 > 0.0 ) {
                      A_trans = OldMomentum.cross(theFacetNormal);
                      A_trans = A_trans.unit();
                   } else {
                      A_trans  = OldPolarization;
                   }
                   A_paral   = NewMomentum.cross(A_trans);
                   A_paral   = A_paral.unit();

                   if(iTE>0&&iTM>0) {
                     NewPolarization = 
                           -OldPolarization + (2.*EdotN)*theFacetNormal;
                   } else if (iTE>0) {
                     NewPolarization = -A_trans;
                   } else if (iTM>0) {
                     NewPolarization = -A_paral;
                   }

                }

             }

             OldMomentum = NewMomentum;
             OldPolarization = NewPolarization;

           }

        } while (NewMomentum * theGlobalNormal < 0.0);
}
void DsG4OpBoundaryProcess::DielectricDielectric ( ) [private]

Definition at line 694 of file DsG4OpBoundaryProcess.cc.

{
        G4bool Inside = false;
        G4bool Swap = false;

        leap:

        G4bool Through = false;
        G4bool Done = false;

        do {

           if (Through) {
              Swap = !Swap;
              Through = false;
              theGlobalNormal = -theGlobalNormal;
              G4SwapPtr(Material1,Material2);
              G4SwapObj(&Rindex1,&Rindex2);
           }

           if ( theFinish == ground || theFinish == groundbackpainted ) {
                theFacetNormal = 
                             GetFacetNormal(OldMomentum,theGlobalNormal);
           }
           else {
                theFacetNormal = theGlobalNormal;
           }

           G4double PdotN = OldMomentum * theFacetNormal;
           G4double EdotN = OldPolarization * theFacetNormal;

           cost1 = - PdotN;
           if (std::abs(cost1) < 1.0-kCarTolerance){
              sint1 = std::sqrt(1.-cost1*cost1);
              sint2 = sint1*Rindex1/Rindex2;     // *** Snell's Law ***
           }
           else {
              sint1 = 0.0;
              sint2 = 0.0;
           }

           if (sint2 >= 1.0) {

              // Simulate total internal reflection

              if (Swap) Swap = !Swap;

              theStatus = TotalInternalReflection;

              if ( theModel == unified && theFinish != polished )
                                                     ChooseReflection();

              if ( theStatus == LambertianReflection ) {
                 DoReflection();
              }
              else if ( theStatus == BackScattering ) {
                 NewMomentum = -OldMomentum;
                 NewPolarization = -OldPolarization;
              }
              else {

                 PdotN = OldMomentum * theFacetNormal;
                 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
                 EdotN = OldPolarization * theFacetNormal;
                 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;

              }
           }
           else if (sint2 < 1.0) {

              // Calculate amplitude for transmission (Q = P x N)

              if (cost1 > 0.0) {
                 cost2 =  std::sqrt(1.-sint2*sint2);
              }
              else {
                 cost2 = -std::sqrt(1.-sint2*sint2);
              }

              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;
              }

              G4double s1 = Rindex1*cost1;
              G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
              G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
              G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
              G4double s2 = Rindex2*cost2*E2_total;

              G4double TransCoeff;

              if (cost1 != 0.0) {
                 TransCoeff = s2/s1;
              }
              else {
                 TransCoeff = 0.0;
              }

              G4double E2_abs, C_parl, C_perp;

              if ( !G4BooleanRand(TransCoeff) ) {

                 // Simulate reflection

                 if (Swap) Swap = !Swap;

                 theStatus = FresnelReflection;

                 if ( theModel == unified && theFinish != polished )
                                                     ChooseReflection();

                 if ( theStatus == LambertianReflection ) {
                    DoReflection();
                 }
                 else if ( theStatus == BackScattering ) {
                    NewMomentum = -OldMomentum;
                    NewPolarization = -OldPolarization;
                 }
                 else {

                    PdotN = OldMomentum * theFacetNormal;
                    NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;

                    if (sint1 > 0.0) {   // incident ray oblique

                       E2_parl   = Rindex2*E2_parl/Rindex1 - E1_parl;
                       E2_perp   = E2_perp - E1_perp;
                       E2_total  = E2_perp*E2_perp + E2_parl*E2_parl;
                       A_paral   = NewMomentum.cross(A_trans);
                       A_paral   = A_paral.unit();
                       E2_abs    = std::sqrt(E2_total);
                       C_parl    = E2_parl/E2_abs;
                       C_perp    = E2_perp/E2_abs;

                       NewPolarization = C_parl*A_paral + C_perp*A_trans;

                    }

                    else {               // incident ray perpendicular

                       if (Rindex2 > Rindex1) {
                          NewPolarization = - OldPolarization;
                       }
                       else {
                          NewPolarization =   OldPolarization;
                       }

                    }
                 }
              }
              else { // photon gets transmitted

                 // Simulate transmission/refraction

                 Inside = !Inside;
                 Through = true;
                 theStatus = FresnelRefraction;

                 if (sint1 > 0.0) {      // incident ray oblique

                    G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
                    NewMomentum = OldMomentum + alpha*theFacetNormal;
                    NewMomentum = NewMomentum.unit();
                    PdotN = -cost2;
                    A_paral = NewMomentum.cross(A_trans);
                    A_paral = A_paral.unit();
                    E2_abs     = std::sqrt(E2_total);
                    C_parl     = E2_parl/E2_abs;
                    C_perp     = E2_perp/E2_abs;

                    NewPolarization = C_parl*A_paral + C_perp*A_trans;

                 }
                 else {                  // incident ray perpendicular

                    NewMomentum = OldMomentum;
                    NewPolarization = OldPolarization;

                 }
              }
           }

           OldMomentum = NewMomentum.unit();
           OldPolarization = NewPolarization.unit();

           if (theStatus == FresnelRefraction) {
              Done = (NewMomentum * theGlobalNormal <= 0.0);
           } 
           else {
              Done = (NewMomentum * theGlobalNormal >= 0.0);
           }

        } while (!Done);

        if (Inside && !Swap) {
          if( theFinish == polishedbackpainted ||
              theFinish == groundbackpainted ) {

              if( !G4BooleanRand(theReflectivity) ) {
                DoAbsorption();
              }
              else {
                if (theStatus != FresnelRefraction ) {
                   theGlobalNormal = -theGlobalNormal;
                }
                else {
                   Swap = !Swap;
                   G4SwapPtr(Material1,Material2);
                   G4SwapObj(&Rindex1,&Rindex2);
                }
                if ( theFinish == groundbackpainted )
                                        theStatus = LambertianReflection;

                DoReflection();

                theGlobalNormal = -theGlobalNormal;
                OldMomentum = NewMomentum;

                goto leap;
              }
          }
        }
}
void DsG4OpBoundaryProcess::ChooseReflection ( ) [inline, private]

Definition at line 258 of file DsG4OpBoundaryProcess.h.

{
                 G4double rand = G4UniformRand();
                 if ( rand >= 0.0 && rand < prob_ss ) {
                    theStatus = SpikeReflection;
                    theFacetNormal = theGlobalNormal;
                 }
                 else if ( rand >= prob_ss &&
                           rand <= prob_ss+prob_sl) {
                    theStatus = LobeReflection;
                 }
                 else if ( rand > prob_ss+prob_sl &&
                           rand < prob_ss+prob_sl+prob_bs ) {
                    theStatus = BackScattering;
                 }
                 else {
                    theStatus = LambertianReflection;
                 }
}
void DsG4OpBoundaryProcess::DoAbsorption ( ) [inline, private]

Definition at line 279 of file DsG4OpBoundaryProcess.h.

{
              theStatus = Absorption;

              if ( G4BooleanRand(theEfficiency) ) {
                
                 // EnergyDeposited =/= 0 means: photon has been detected
                 theStatus = Detection;
                 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
              }
              else {
                 aParticleChange.ProposeLocalEnergyDeposit(0.0);
              }

              NewMomentum = OldMomentum;
              NewPolarization = OldPolarization;

//              aParticleChange.ProposeEnergy(0.0);
              aParticleChange.ProposeTrackStatus(fStopAndKill);
}
void DsG4OpBoundaryProcess::DoReflection ( ) [inline, private]

Definition at line 301 of file DsG4OpBoundaryProcess.h.

{
        if ( theStatus == LambertianReflection ) {

          NewMomentum = G4LambertianRand(theGlobalNormal);
          theFacetNormal = (NewMomentum - OldMomentum).unit();

        }
        else if ( theFinish == ground ) {

          theStatus = LobeReflection;
          theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
          G4double PdotN = OldMomentum * theFacetNormal;
          NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;

        }
        else {

          theStatus = SpikeReflection;
          theFacetNormal = theGlobalNormal;
          G4double PdotN = OldMomentum * theFacetNormal;
          NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;

        }
        G4double EdotN = OldPolarization * theFacetNormal;
        NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
}

Member Data Documentation

Definition at line 182 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::OldMomentum [private]

Definition at line 184 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::OldPolarization [private]

Definition at line 185 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::NewMomentum [private]

Definition at line 187 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::NewPolarization [private]

Definition at line 188 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::theGlobalNormal [private]

Definition at line 190 of file DsG4OpBoundaryProcess.h.

G4ThreeVector DsG4OpBoundaryProcess::theFacetNormal [private]

Definition at line 191 of file DsG4OpBoundaryProcess.h.

G4Material* DsG4OpBoundaryProcess::Material1 [private]

Definition at line 193 of file DsG4OpBoundaryProcess.h.

G4Material* DsG4OpBoundaryProcess::Material2 [private]

Definition at line 194 of file DsG4OpBoundaryProcess.h.

G4OpticalSurface* DsG4OpBoundaryProcess::OpticalSurface [private]

Definition at line 196 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::Rindex1 [private]

Definition at line 198 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::Rindex2 [private]

Definition at line 199 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::cost1 [private]

Definition at line 201 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::cost2 [private]

Definition at line 201 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::sint1 [private]

Definition at line 201 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::sint2 [private]

Definition at line 201 of file DsG4OpBoundaryProcess.h.

Definition at line 203 of file DsG4OpBoundaryProcess.h.

G4OpticalSurfaceModel DsG4OpBoundaryProcess::theModel [private]

Definition at line 205 of file DsG4OpBoundaryProcess.h.

G4OpticalSurfaceFinish DsG4OpBoundaryProcess::theFinish [private]

Definition at line 207 of file DsG4OpBoundaryProcess.h.

Definition at line 209 of file DsG4OpBoundaryProcess.h.

Definition at line 210 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::prob_sl [private]

Definition at line 211 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::prob_ss [private]

Definition at line 211 of file DsG4OpBoundaryProcess.h.

G4double DsG4OpBoundaryProcess::prob_bs [private]

Definition at line 211 of file DsG4OpBoundaryProcess.h.

G4int DsG4OpBoundaryProcess::iTE [private]

Definition at line 213 of file DsG4OpBoundaryProcess.h.

G4int DsG4OpBoundaryProcess::iTM [private]

Definition at line 213 of file DsG4OpBoundaryProcess.h.

Definition at line 215 of file DsG4OpBoundaryProcess.h.

Definition at line 217 of file DsG4OpBoundaryProcess.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4