Hello,

 we are currently using GEANT 4 to simulate the light collection of a
Liquid Xenon detector.
 In analysing the code of the G4OpBoundaryProcesses routine we've come out
with some doubts in the DielectricDielectric() method, that we now think
are actual errors in the code.
 Next we present the code we have rewritten, including, within comments,
the lines we think are incorrect and an explanantion of why we think
they are incorrect.

-----------------------------------------------------------------------
else if (sint2 < 1.0) {

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

  if (cost1 > 0.0) {
    cost2 =  sqrt(1-sint2*sint2);
  }
  else {
    cost2 = -sqrt(1-sint2*sint2);
  }
  G4ThreeVector A_trans, Atrans, E1pp, E1pl;
  G4double E1_perp, E1_parl;
  // Modified code begin        <---------------------------------------
  G4ThreeVector eta_, eta, omega_, omega;
  // Modified code end          <---------------------------------------

  if (sint1 > 0.0) {
    //A_trans = OldMomentum.cross(theFacetNormal);->Vector normal to the
                                                //plane, pointing down
    //Atrans  = A_trans.unit();
    //E1_perp = OldPolarization * Atrans;
    //E1pp    = E1_perp * Atrans;
    //E1pl    = OldPolarization - E1pp; -> This vector depends on
                                        //OldPolarization. Some times
                                        //it points left, other times
                                        //right (in the OldMomentum
                                        //frame)
    //E1_parl = E1pl.mag();

    //Modified code begin       <---------------------------------------
    eta_ = OldMomentum.cross(theFacetNormal);   //vector normal
                                                //to the plane,
                                                //pointing down
    eta  = eta_.unit(); //unit vector
    E1_perp = OldPolarization * eta;
    E1pp    = E1_perp * eta;    //perpendicular component
                                //of the polarization
    omega_ = OldMomentum.cross(eta);    //vector parallel to the plane
                                        //and perpendicular to the
                                        //OldMomentum, pointing left
                                        //from it
    omega = omega_.unit(); //unit vector
    E1_parl = OldPolarization.dot(omega);//projection of the
                                        //polarization in the
                                        //parallel direction
    E1pl = E1_parl*omega; //parallel component of the polarization
    //Modified code end         <---------------------------------------
  }
  else {
    //Modified code begin       <---------------------------------------
    omega_  = OldPolarization;
    omega = omega_.unit(); //To assure the vector is unitary
    //Modified code end <---------------------------------------

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

  G4ThreeVector Refracted, Deflected;
  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;
        //Refracted = theFacetNormal + PdotN * NewMomentum;->This vector
                                                //always points right.
                                                //It's not unitary
        //Modified code begin   <---------------------------------------
        omega_ = NewMomentum.cross(eta); //plane vector, perpendicular
                                           //to the new momentum and
                                           //pointing left from it
        omega = omega_.unit(); //unit vector
        E2_abs    = sqrt(E2_total);
        C_parl    = E2_parl/E2_abs;
        C_perp    = E2_perp/E2_abs;

        NewPolarization = C_parl*omega + C_perp*eta;
        //Modified code end     <---------------------------------------
        //NewPolarization = C_parl*Refracted - C_perp*A_trans;->First of
                                        //all, A_trans is not a unitary
                                        //vector. And the minus sign
                                        //causes a phase inversion in
                                        //perpendicular component.
      }
      else if (Rindex2 > Rindex1) { // incident ray perpendicular

        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);
    Deflected  = OldMomentum + alpha*theFacetNormal;
    NewMomentum = Deflected.unit();
    PdotN = -cost2;
    //Refracted  = theFacetNormal - PdotN*NewMomentum;->Same problems as
                                                //in the reflection code
    //Modified code begin       <---------------------------------------
    omega_ = NewMomentum.cross(eta); //plane vector, perpendicular
                                       //to the new momentum and
                                       //pointing left from it
    omega = omega_.unit(); //unit vector
    E2_abs     = sqrt(E2_total);
    C_parl     = E2_parl/E2_abs;
    C_perp     = E2_perp/E2_abs;
    NewPolarization = C_parl*omega + C_perp*eta;
    //Modified code end         <---------------------------------------
    //NewPolarization = C_parl*Refracted + C_perp*A_trans;->Again, A_trans
                                                        //is not unitary
  }
  else {  // incident ray perpendicular

    NewMomentum = OldMomentum;
    NewPolarization = OldPolarization;

  }
}

OldMomentum = NewMomentum;
OldPolarization = NewPolarization;

-----------------------------------------------------------------------

 In the next few lines is the output for the calculus of the
polarization from the old and new codes, in the same conditions:

                        Old Code

Pol                                             ||Pol||
------------------------------------------------------------------------
(-0.147677,-0.179211,-0.0479356)                0.237114
(-0.104166,-0.197365,-0.0503692)                0.228779
(0.225067,0.247816,-0.0315774)                  0.336252
(0.380579,0.344107,-0.365216)                   0.629788
(0.601605,0.114947,-0.0644966)                  0.615875
------------------------------------------------------------------------

                        New Code

Pol                                             ||Pol||
------------------------------------------------------------------------
(-0.649364,-0.734408,-0.197413)                 1
(-0.144312,0.965069,0.21867)                    1
(0.933388,0.348001,0.0876508)                   1
(0.589937,0.569892,-0.572012)                   1
(0.994365,-0.0417125,0.0974548)                 1
------------------------------------------------------------------------

 Please let us know what you think about these changes.
 Thank you very much

 Alexandre Lindote (alex@lipc.fis.uc.pt)
 José Pinto da Cunha (jpinto@lipc.fis.uc.pt)