Changeset 19638 for trunk/Mars/msimreflector

Timestamp:

09/15/19 14:22:38 (5 years ago)

Author:

tbretz

Message:

Updated the algorithm and added comments (still not final, I think)

Location:

trunk/Mars/msimreflector

Files:

• MFresnelLens.cc (modified) (10 diffs)
• MFresnelLens.h (modified) (2 diffs)

trunk/Mars/msimreflector/MFresnelLens.cc

@@ -27 +27 @@
 //  MFresnelLens
 //
-//  For some details on definitions please refer to
-//  https://application.wiley-vch.de/berlin/journals/op/07-04/OP0704_S52_S55.pdf
+// For some details on definitions please refer to
+// https://application.wiley-vch.de/berlin/journals/op/07-04/OP0704_S52_S55.pdf
+//
+// The HAWC's Eye lens is an Orafol SC943
+// https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
+//
+// A good description on ray-tracing can be found here
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -41 +47 @@
 #include "MReflection.h"
 
+#include "MMath.h"
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -48 +56 @@
 using namespace std;
 
+// ==========================================================================
+
+enum exception_t
+{
+    kValidRay = 0,
+
+    kStrayUpgoing,
+    kOutsideRadius,
+    kNoSurfaceFound,
+    kStrayDowngoing,
+    kAbsorbed,
+
+    kFoundSurfaceUnavailable,
+
+    kInvalidOrigin,
+    kTransitionError,
+
+    kEnter = 1000,
+    kLeave = 2000,
+};
+
+enum surface_t
+{
+    kPhotonHasLeft = 0,
+
+    kEntrySurface,
+    kSlopeSurface,
+    kDraftSurface,
+    kExitSurface,
+
+    kMaterial = 5,
+
+    kNoSurface = 9
+};
+
+
+class raytrace_exception : public runtime_error
+{
+protected:
+    int fError;
+    int fOrigin;
+    int fSurface;
+
+public:
+    raytrace_exception(const int &_id, const int &_origin, const int &_surface, const string& what_arg) :
+        runtime_error(what_arg), fError(_id), fOrigin(_origin), fSurface(_surface)
+    {
+    }
+
+    raytrace_exception(const int &_id, const int &_origin, const int &_surface, const char* what_arg) :
+        runtime_error(what_arg), fError(_id), fOrigin(_origin), fSurface(_surface)
+    {
+    }
+
+    int id() const { return fError + fSurface*10 + fOrigin*100; }
+    int error() const { return fError; }
+    int origin() const { return fOrigin; }
+    int surface() const { return fSurface; }
+};
+
+class raytrace_error : public raytrace_exception
+{
+public:
+    raytrace_error(const int &_id, const int &_origin, const int &_surface, const string& what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+    raytrace_error(const int &_id, const int &_origin, const int &_surface, const char* what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+};
+class raytrace_info : public raytrace_exception
+{
+public:
+    raytrace_info(const int &_id, const int &_origin, const int &_surface, const string& what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+    raytrace_info(const int &_id, const int &_origin, const int &_surface, const char* what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+};
+class raytrace_user : public raytrace_exception
+{
+public:
+    raytrace_user(const int &_id, const int &_origin, const int &_surface, const string& what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+    raytrace_user(const int &_id, const int &_origin, const int &_surface, const char* what_arg) :
+        raytrace_exception(_id, _origin, _surface, what_arg) { }
+};
+
+// ==========================================================================
+
+
 // --------------------------------------------------------------------------
 //
 // Default constructor
 //
-MFresnelLens::MFresnelLens(const char *name, const char *title)
+MFresnelLens::MFresnelLens(const char *name, const char *title) :
+    fPSF(0), fSlopeAbsorption(false), fDraftAbsorption(false),
+    fBottomReflection(true), fDisableMultiEntry(false), fFresnelReflection(true),
+    fMinHits(0), fMaxHits(0)
 {
     fName  = name  ? name  : "MFresnelLens";
     fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";
 
-    fMaxR = 55/2.;
-}
+    // Default: Orafol SC943
+
+    DefineLens();
+}
+
+// ==========================================================================
+
+// --------------------------------------------------------------------------
+//
+// Default ORAFOL SC943
+//
+// Focal Length:   F = 50.21 cm
+// Diameter:       D = 54.92 cm
+// Groove width:   w =  0.01 cm
+// Lens thickness: h =  0.25 cm
+//
+// Default wavelength: 546 nm
+//
+void MFresnelLens::DefineLens(double F, double D, double w, double h, double lambda)
+{
+    fR = D/2;  // [cm] Lens radius
+    fW = w;    // [cm] Width of a single groove
+    fH = h;    // [cm] Thickness of lens
+    fF = F;    // [cm] focal length (see also MGeomCamFAMOUS!)
+
+    fLambda = lambda;
+
+    fN = MFresnelLens::RefractiveIndex(fLambda);  // Lens
+
+    // Velocity of light within the lens material [cm/ns]
+    // FIXME: Note that for the correct conversion in Transmission()
+    // also the speed in the surrounding medium has to be taken correctly
+    // into account (here it is assumed to be air with N=1
+    fVc = fN/(TMath::C()*100/1e9); // cm/ns
+
+    InitGeometry(fR, fW, fN, fF, fH);
+}
+
+// --------------------------------------------------------------------------
+//
+// Precalculate values such as the intersection points inside the grooves,
+// the angle of the slope and draft surface and the corresponding tangents.
+//
+void MFresnelLens::InitGeometry(double maxr, double width, double N0, double F, double d)
+{
+    const uint32_t num = TMath::CeilNint(maxr/width);
+
+    fGrooves.resize(num);
+
+    for (int i=0; i<num; i++)
+    {
+        const double r0 = i*width;
+        const double rc = i*width + width/2;
+        const double r1 = i*width + width;
+
+        // Slope angle of the reflecting surface alpha
+        // Angle of the draft surface psi
+        const double alpha = -MFresnelLens::SlopeAngle(rc, F, N0, d);   // w.r.t. x  [30]
+        const double psi   =  MFresnelLens::DraftAngle(r1);             // w.r.t. z  [ 5]
+
+        const double tan_alpha = tan(alpha);
+        const double tan_psi   = tan(psi);
+
+        fGrooves[i].slope.z =  r0*tan_alpha;
+        fGrooves[i].draft.z = -r1/tan_psi;
+
+        fGrooves[i].slope.theta = TMath::Pi()/2-alpha;  // w.r.t. +z [ 60]
+        fGrooves[i].draft.theta = -psi;                 // w.r.t. +z [- 5]
+
+        fGrooves[i].slope.tan_theta = tan(fGrooves[i].slope.theta);
+        fGrooves[i].draft.tan_theta = tan(fGrooves[i].draft.theta);
+
+        fGrooves[i].slope.tan_theta2 = fGrooves[i].slope.tan_theta*fGrooves[i].slope.tan_theta;
+        fGrooves[i].draft.tan_theta2 = fGrooves[i].draft.tan_theta*fGrooves[i].draft.tan_theta;
+
+        fGrooves[i].slope.theta_norm = TMath::Pi()/2-fGrooves[i].slope.theta; // [ 30] 
+        fGrooves[i].draft.theta_norm = TMath::Pi()/2-fGrooves[i].draft.theta; // [ 95] 
+
+        const double dr = width/(tan_alpha*tan_psi+1);
+
+        fGrooves[i].r = r0 + dr;
+
+        const double z = -dr*tan_alpha;
+
+        fGrooves[i].slope.h = z;
+        fGrooves[i].draft.h = z;
+
+        if (z<-fH)
+            gLog << warn << "Groove " << i << " deeper (" << z << ") than thickness of lens material (" << fH << ")." << endl;
+    }
+
+    fMaxR = (num+1)*width;
+}
+
+// --------------------------------------------------------------------------
+//
+// Reads the transmission curve from a file
+// (tranmission in percent versus wavelength in nanometers)
+//
+// The transmission curve is used to calculate the absorption lengths.
+// Therefore the thickness for which the tranission curve is valid is
+// required (in cm).
+//
+// The conversion can correct for fresnel reflection at the entry and exit
+// surface assuming that the outside material during the measurement was air
+// (n=1.0003) and the material in PMMA. Correction is applied when
+// correction is set to true <default>.
+//
+// If no valid data was read, 0 is returned. -1 is returned if any tranmission
+// value read from the file is >1. If the fresnel correction leads to a value >1,
+// the value is set to 1. The number of valid data points is returned.
+//
+Int_t MFresnelLens::ReadTransmission(const TString &file, float thickness, bool correction)
+{
+    TGraph transmission(file);
+
+    /*
+    double gx_min, gx_max, gy_min, gy_max;
+    absorption.ComputeRange(gx_min, gy_min, gx_max, gy_max);
+    if (lambda<gx_min || lambda>gx_max)
+    {
+        cout << "Invalid wavelength" << endl;
+        return;
+    }*/
+
+    if (transmission.GetN()==0)
+        return 0;
+
+    for (int i=0; i<transmission.GetN(); i++)
+    {
+        // Correct transmission for Fresnel reflection on the surface
+        const double lambda = transmission.GetX()[i];;
+
+        double trans = transmission.GetY()[i]/100;
+        if (trans>1)
+        {
+            gLog << err << "Transmission larger than 1." << endl;
+            return -1;
+        }
+
+        if (correction)
+        {
+            // Something like this is requried if correction
+            // for optical boundaries is necessary
+            const double n0 = MFresnelLens::RefractiveIndex(lambda);
+
+            // FIXME: Make N_air a variable
+            const double r0 = (n0-1.0003)/(n0+1.0003);
+            const double r2 = r0*r0;
+
+            trans *= (1+r2)*(1+r2);
+
+            if (trans>1)
+            {
+                gLog << warn << "Transmission at " << lambda << "nm (" << trans << ") after Fresnel correction larger than 1." << endl;
+                trans = 1;
+            }
+        }
+
+        // convert to absorption length (FIMXE: Sanity check)
+        transmission.GetY()[i] = -thickness/log(trans>0.999 ? 0.999 : trans);
+    }
+
+    fAbsorptionLength = MSpline3(transmission);
+
+    return fAbsorptionLength.GetNp();
+}
+
+// ==========================================================================
 
 // --------------------------------------------------------------------------
@@ -78 +344 @@
 
     return sqrt(1+c0+c1+c2);
+}
+
+// --------------------------------------------------------------------------
+//
+// A Fresnel lens with parabolic surface calculated with the sagittag
+// function (k=-1) and a correction for the thickness of the lens
+// on the curvature. See also PhD thesis, Tim Niggemann ch. 7.1.1.
+//
+// Parameters are:
+// The distance from the center r
+// The focal length to be achieved F
+// The refractive index of the outer medium (usually air) n0
+// The refractive index of the lens material (e.g. PMMA) n1
+// The thichness of the lens d
+//
+// r, F and d have to be in the same units.
+//
+// Return the slope angle alpha [rad]. The Slope angle is defined with
+// respect to the plane of the lens. (0 at the center, decreasing
+// with increasing radial distance)
+//
+double MFresnelLens::SlopeAngleParabolic(double r, double F, double n0, double n1, double d)
+{
+    // PhD, Tim Niggemann, ch 7.1.1
+    const double rn = n1/n0;
+    const double c = (rn - 1) * (F + d/rn);
+    return -atan(r/c);
+}
+
+// --------------------------------------------------------------------------
+//
+// A Fresnel lens with an optimized parabolic surface calculated with
+// the sagittag function (k=-1) and fitted coefficients according
+// to Master thesis, Eichler.
+//
+// Note that for this setup other parameters must be fixed
+//
+// Parameters are:
+// The distance from the center r
+//
+// r is in cm.
+//
+// Return the slope angle alpha [rad]. The Slope angle is defined with
+// respect to the plane of the lens. (0 at the center, decreasing
+// with increasing radial distance)
+//
+double MFresnelLens::SlopeAngleAspherical(double r)
+{
+    // Master, Eichler [r/cm]
+    return -atan( r/26.47
+                 +2*1.18e-4  * 1e1*r
+                 +4*1.34e-9  * 1e3*r*r*r
+                 +6*9.52e-15 * 1e5*r*r*r*r*r
+                 -8*2.04e-19 * 1e7*r*r*r*r*r*r*r);
 }
 
@@ -92 +412 @@
 // Here, a thin lens is assumed
 //
-// The HAWC's Eye lens is an Orafol SC943
-// https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
-//
 // sin(omega) = r / sqrt(r^2+F^2)
 // tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ]
@@ -104 +421 @@
 // distance)
 //
-double MFresnelLens::SlopeAngle(double r, double F, double n)
-{
+double MFresnelLens::SlopeAngleOptimized(double r, double F, double n)
+{
+    // Use F+d/2
     double so = r / sqrt(r*r + F*F);
     return atan(so / (1-sqrt(n*n - so*so))); // alpha<0, Range [0deg; -50deg]
+}
+
+// --------------------------------------------------------------------------
+//
+// Currently calles SlopeAngleParabolic(r, F, 1, n, d)
+//
+double MFresnelLens::SlopeAngle(double r, double F, double n, double d)
+{
+    return SlopeAngleParabolic(r, F, 1.0003, n, d);
 }
 
@@ -126 +453 @@
     return (3 + r*0.473)*TMath::DegToRad(); // Range [0deg; 15deg]
 }
+
+// ==========================================================================
 
 // --------------------------------------------------------------------------
@@ -149 +478 @@
 }
 
-struct Groove
-{
-    double fTheta;
-    double fPeakZ;
-
-    double fTanTheta;
-    double fTanTheta2[3];
-
-    Groove() { }
-    Groove(double theta, double z) : fTheta(theta), fPeakZ(z)
-    {
-        fTanTheta = tan(theta);
-
-        fTanTheta2[0] = fTanTheta*fTanTheta;
-        fTanTheta2[1] = fTanTheta*fTanTheta * fPeakZ;
-        fTanTheta2[2] = fTanTheta*fTanTheta * fPeakZ*fPeakZ;
-    }
-};
-
-double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness)
-{
-/*
-    Z: position of peak
-
-    r_cone(z) = (Z-z)*tan(theta)
-
-    (x)   (p.x)          (u.x)
-    (y) = (p.y)  + dz *  (u.y)
-    (z)   (p.z)          (u.z)
-
-    r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2)
-
-    Solved with wmmaxima:
-
-    solve((H-z)^2*t^2 = (px+(z-pz)/uz*ux)^2+(py+(z-pz)/uz*uy)^2, z);
-
-     [
-      z=-(uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)-H*t^2*uz^2+(-py*uy-px*ux)*uz+pz*uy^2+pz*ux^2)/(t^2*uz^2-uy^2-ux^2)
-      z= (uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)+H*t^2*uz^2+( py*uy+px*ux)*uz-pz*uy^2-pz*ux^2)/(t^2*uz^2-uy^2-ux^2)
-     ]
-*/
+// ==========================================================================
+
+// FIXME: The rays could be 'reflected' inside the material
+// (even though its going out) or vice versa
+static double RandomTheta(double psf)
+{
+    return psf>0 ? MMath::RndmPSF(psf)/2 : 0;
+}
+
+// FIXME: The rays could be 'reflected' inside the material
+// (even though its going out) or vice versa
+static double RandomPhi(double r, double psf)
+{
+    return psf>0 ? MMath::RndmPSF(psf)/2 : 0;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Calculate the intersection point beweteen a line defined by the position p
+// and the direction u and a cone defined by the object cone.
+//
+// Z: position of peak of cone
+// theta: opening angle of cone
+//
+// Distance r of cone surface at given z from z-axis
+// r_cone(z) = (Z-z)*tan(theta)
+//
+// Equalition of line
+//  (x)   (p.x)          (u.x/u.z)
+//  (y) = (p.y)  + dz *  (u.y/u.z)
+//  (z)   (p.z)          (   1   )
+//
+// Normalization
+//  U.x := u.x/u.z
+//  U.y := u.y/u.z
+//
+// Distance of line at given z from z-axis
+//  r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2)   with   dz = z-p.z
+//
+// Equation to be solved
+//  r_cone(z) = r_line(z)
+//
+// Solved with wxmaxima:
+//
+//  [0] solve((px+(z-pz)*Ux)^2+(py+(z-pz)*Uy)^2= ((Z-z)*t)^2, z);
+//
+//  z= (sqrt(((Uy^2+Ux^2)*pz^2+(-2*Uy*py-2*Ux*px-2*Z*Uy^2-2*Z*Ux^2)*pz+py^2+2*Z*Uy*py+px^2+2*Z*Ux*px+Z^2*Uy^2+Z^2*Ux^2)*t^2-Ux^2*py^2+2*Ux*Uy*px*py-Uy^2*px^2)+Z*t^2+(-Uy^2-Ux^2)*pz+Uy*py+Ux*px)/(t^2-Uy^2-Ux^2),
+//  z=-(sqrt(((Uy^2+Ux^2)*pz^2+(-2*Uy*py-2*Ux*px-2*Z*Uy^2-2*Z*Ux^2)*pz+py^2+2*Z*Uy*py+px^2+2*Z*Ux*px+Z^2*Uy^2+Z^2*Ux^2)*t^2-Ux^2*py^2+2*Ux*Uy*px*py-Uy^2*px^2)-Z*t^2+( Uy^2+Ux^2)*pz-Uy*py-Ux*px)/(t^2-Uy^2-Ux^2)
+//
+double MFresnelLens::CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Cone &cone) const
+{
+    const double &Z = cone.z;
 
     const double Ux = u.X()/u.Z();
@@ -198 +539 @@
     const double pz = p.Z();
 
-    const double H   = g.fPeakZ;
-
-    const double t   = g.fTanTheta;
-    const double t2  = t*t;
+    //const double &t  = cone.tan_theta;
+    const double &t2 = cone.tan_theta2;
 
     const double Ur2 = Ux*Ux + Uy*Uy;
@@ -209 +548 @@
     const double cr2 = Ux*py - Uy*px;
 
-    const double a   = t2 - Ur2;
-    const double b   = Ur2*pz - Up2 - H*t2;
-
-    const double h   = H-pz;
-    const double h2  = h*h;
+    const double a = t2 - Ur2;
+    const double b = Ur2*pz - Up2 - Z*t2;
+
+    const double h  = Z-pz;
+    const double h2 = h*h;
 
     // [ -b +-sqrt(b^2 - 4 ac) ] / [ 2a ]
 
     const double radix = (Ur2*h2 + 2*Up2*h + pr2)*t2 - cr2*cr2;
-
     if (radix<0)
         return 0;
 
-    double dz[2] =
-    {
-        (-b+sqrt(radix))/a,
-        (-b-sqrt(radix))/a
+    const double sqrt_radix = sqrt(radix);
+
+    const double dz[2] =
+    {
+        (+sqrt_radix - b)/a,
+        (-sqrt_radix - b)/a
     };
 
-     if (dz[0]<0 && dz[0]>fThickness)
-        return dz[0];
-
-    if (dz[1]<0 && dz[1]>fThickness)
-        return dz[1];
-
-    return 0;
+    // Return the closest solution inside the allowed range
+    // which is in the direction of movement
+
+    const double &H = cone.h;
+
+    const bool is_inside0 = dz[0]>=H && dz[0]<0;
+    const bool is_inside1 = dz[1]>=H && dz[1]<0;
+
+    // FIXME: Simplify!
+    if (!is_inside0 && !is_inside1)
+        return 0;
+
+    // Only dz[0] is in the right z-range
+    if (is_inside0 && !is_inside1)
+    {
+        // Check if dz[0] is in the right direction
+        if ((u.Z()>=0 && dz[0]>=p.Z()) ||
+            (u.Z()< 0 && dz[0]< p.Z()))
+            return dz[0];
+
+        return 0;
+    }
+
+    // Only dz[1] is in the right z-range
+    if (!is_inside0 && is_inside1)
+    {
+        // Check if dz[1] is in the right direction
+        if ((u.Z()>=0 && dz[1]>=p.Z()) ||
+            (u.Z()< 0 && dz[1]< p.Z()))
+            return dz[1];
+
+        return 0;
+    }
+
+    /*
+    if (is_inside0^is_inside1)
+    {
+        if (u.Z()>=0)
+            return dz[0]>p.Z() ? dz[0] : dz[1];
+        else
+            return dz[0]<p.Z() ? dz[0] : dz[1];
+    }*/
+
+
+    // dz[0] and dz[1] are in the right range
+    // return the surface which is hit first
+
+    // moving upwards
+    if (u.Z()>=0)
+    {
+        // Both solution could be correct
+        if (dz[0]>=p.Z() && dz[1]>=p.Z())
+            return std::min(dz[0], dz[1]);
+
+        // only one solution can be correct
+        return dz[0]>=p.Z() ? dz[0] : dz[1];
+    }
+    else
+    {
+        // Both solution could be correct
+        if (dz[0]<p.Z() && dz[1]<p.Z())
+            return std::max(dz[0], dz[1]);
+
+        // only one solution can be correct
+        return dz[0]<p.Z() ? dz[0] : dz[1];
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+// Find the peak (draft+slope) which will be hit by the photon which
+// is defined by position p and direction u. ix gives the index of the groove
+// to originate the search from.
+//
+// Returns the index of the groove to which the surface belongs, -1 if no
+// matching surface was found.
+//
+int MFresnelLens::FindPeak(size_t ix, const MQuaternion &p, const MQuaternion &u) const
+{
+    // ---------------------------
+    // check for first groove first
+    if (ix==0)
+    {
+        const auto test = p.fVectorPart + (fGrooves[0].slope.h-p.Z())/u.Z()*u.fVectorPart;
+        if (test.XYvector().Mod()<fGrooves[0].r)
+            return 0;
+    }
+
+    // r = sqrt( (px + t*ux) + (py + t*uy)^2 )
+    // dr/dt = (2*uy*(dz*uy+py)+2*ux*(dz*ux+px))/(2*sqrt((dz*uy+py)^2+(dz*ux+px)^2))
+    // dr/dt = (uy*py + ux*px)/sqrt(py^2+px^2)
+    const bool outgoing = u.X()*p.X() + u.Y()*p.Y() > 0; // r is (at least locally) increasing
+
+    // ---------------------------
+    const double Ux = u.X()/u.Z();
+    const double Uy = u.Y()/u.Z();
+
+    const double px = p.X();
+    const double py = p.Y();
+    const double pz = p.Z();
+
+    const double Ur2 = Ux*Ux + Uy*Uy;
+    const double cr2 = Ux*py - Uy*px;
+    const double pr2 = px*px + py*py;
+    const double Up2 = Ux*px + Uy*py;
+
+    //for (int i=1; i<fGrooves.size(); i++)
+
+    // To speed up the search, search first along the radial moving direction of
+    // the photon. If that was not successfull, try in the opposite direction.
+    // FIXME: This could still fail in some very rare cases, for some extremely flat trajectories
+    for (int j=0; j<2; j++)
+    {
+        const bool first = j==0;
+
+        const int step = outgoing ^ !first ? 1 : -1;
+        const int end  = outgoing ^ !first ? fGrooves.size() : 1;
+        const int beg  = j==0 ? ix : ix+(step);
+
+        for (int i=beg; i!=end; i+=step)
+        {
+            const Groove &groove1 = fGrooves[i-1];
+            const Groove &groove2 = fGrooves[i];
+
+            const double &z1 = groove1.draft.h;
+            const double &z2 = groove2.slope.h;
+
+            const double &r1 = groove1.r;
+            const double &r2 = groove2.r;
+
+            Cone cone;
+            cone.tan_theta  = -(r2-r1)/(z2-z1);
+            cone.tan_theta2 = cone.tan_theta*cone.tan_theta;
+            cone.z          = z1 + r1/cone.tan_theta;
+
+            const double &Z  = cone.z;
+            const double &t2 = cone.tan_theta2;
+
+            const double a = t2 - Ur2;
+            const double b = Ur2*pz - Up2 - Z*t2;
+
+            const double h  = Z-pz;
+            const double h2 = h*h;
+
+            // [ -b +-sqrt(b^2 - 4 ac) ] / [ 2a ]
+
+            const double radix = (Ur2*h2 + 2*Up2*h + pr2)*t2 - cr2*cr2;
+            if (radix<0)
+                continue;
+
+            const double sqrt_radix = sqrt(radix);
+
+            const double dz[2] =
+            {
+                (+sqrt_radix - b)/a,
+                (-sqrt_radix - b)/a
+            };
+
+            if (dz[0]>=z2 && dz[0]<=z1)
+                return i;
+
+            if (dz[1]>=z2 && dz[1]<=z1)
+                return i;
+        }
+    }
+
+    return -1;
+}
+
+// --------------------------------------------------------------------------
+//
+// If no transmission was given returns true. Otherwaise calculates the
+// absorption length for a flight time dt in the material and a photon
+// with wavelength lambda. The flight time is converted to a geometrical
+// using the speed of light in the medium.
+//
+// Returns true if the poton passed, false if it was absorbed.
+//
+bool MFresnelLens::Transmission(double dt, double lambda) const
+{
+    if (fAbsorptionLength.GetNp()==0)
+        return true;
+
+    // FIXME: Speed up!
+    const double alpha = fAbsorptionLength.Eval(lambda);
+
+    // We only have the travel time, thus we have to convert back to distance
+    // Note that the transmission coefficients are w.r.t. to geometrical
+    // distance not light-travel distance. Thus the distance has to be corrected
+    // for the corresponding refractive index of the material.
+    const double cm = dt/fVc;  
+
+    const double trans = exp(-cm/alpha);
+    return gRandom->Uniform()<trans;
 }
 
 /*
-double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness)
-{
-   // p is the position  in the plane  of the lens (px, py, 0,  t)
-   // u is the direction of the photon             (ux, uy, dz, dt)
-
-   // Assuming a cone with peak at z and opening angle theta
-
-   // Radius at distance z+dz from the tip and at distance r from the axis of the cone
-   // r = (z+dz*u.z) * tan(theta)
-
-   // Defining
-   // zc := z+dz
-
-   // Surface of the cone
-   // (X)         ( tan(theta) * sin(phi) )
-   // (Y)  = zc * ( tan(theta) * cos(phi) )
-   // (Z)         (          1            )
-
-   // Line
-   // (X)        (u.x)   (p.x)
-   // (Y) = dz * (u.y) + (p.y)
-   // (Z)        (u.z)   (p.z)
-
-   // Equality
-   //
-   // X^2+Y^2 = r^2
-   //
-   // (dz*u.x+p.x)^2 + (dz*u.y+p.y)^2 = (z+dz)^2 * tan(theta)^2
-   //
-   // dz^2*ux^2 + px^2 + 2*dz*ux*px + dz^2*uy^2 + py^2 + 2*dz*uy*py = (z^2*tan(theta)^2 + dz^2*tan(theta)^2 + 2*z*dz*tan(theta)^2
-   //
-   // dz^2*(ux^2 + uy^2 - tan(theta)^2)  +  2*dz*(ux*px + uy*py - z*tan(theta)^2)  +  (px^2+py^2 - z^2*tan(theta)^2)  =  0
-
-   // t   := tan(theta)
-   // a   := ux^2    + uy^2    -     t^2
-   // b/2 := ux*px   + uy*py   - z  *t^2 := B
-   // c   :=    px^2 +    py^2 - z^2*t^2
-
-   // dz = [ -b +- sqrt(b^2 - 4*a*c) ] / [2*a]
-   // dz = [ -B +- sqrt(B^2 -   a*c) ] /    a
-
-    const double a = u.R2()                    - g.fTanTheta2[0];
-    const double B = u.XYvector()*p.XYvector() - g.fTanTheta2[1];
-    const double c = p.R2()                    - g.fTanTheta2[2];
-
-    const double B2 = B*B;
-    const double ac = a*c;
-
-    if (B2<ac)
+// surface=0 : incoming ray
+// surface=1 : slope
+// surface=2 : draft
+// surface=3 : bottom
+int MFresnelLens::EnterGroove(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir) const
+{
+    const double rx = pos.R();
+
+    if (surface==3)
+    {
+        //cout << "Bottom as origin invalid" << endl;
+        throw -1;
+
+    }
+    if (rx>=fR)
+    {
+        //cout << "Left the lens radius (enter)" << endl;
+        throw -2;
+    }
+    //if (dir.Z()>0)
+    //{
+    //    cout << "Upgoing, outside of the material" << endl;
+    //    PropagateZ(pos, dir, dir.Z()>0 ? 3 : -3);
+    //    return -1;
+    //}
+
+
+   // Calculate the ordinal number of the groove correpsonding to rx
+    const int ix = TMath::FloorNint(rx/fW);
+
+    // Photons was just injected (test both surfaces) or came from the other surface
+    if (surface==0 || surface==2)
+    {
+        // Get the possible intersection point with the slope angle
+        const double z1 = CalcIntersection(pos, dir, fGrooves[ix].slope);
+
+        // We hit the slope angle
+        if (z1!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z1);
+
+            if (fSlopeAbsorption)
+                throw -100;
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(fGrooves[ix].slope.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, 1, n0);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return EnterGroove(1, n0, lambda, pos, dir)+1;
+
+            // Transition was Refraction - enter
+            if (ret>=3)
+                return LeavePeak(1, n0, lambda, pos, dir, pos.T())+1;
+
+            // Error occured (see ApplyTransition for details)
+            //cout << "ERR[TIR1]" << endl;
+            throw -3;
+        }
+    }
+
+    // Photons was just injected (test both surfaces) or came from the other surface
+    if (surface==0 || surface==1)
+    {
+        const double z2 = CalcIntersection(pos, dir, fGrooves[ix].draft);
+
+        // We hit the draft angle
+        if (z2!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z2);
+
+            if (fDraftAbsorption)
+                throw -101;
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(fGrooves[ix].draft.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, 1, n0);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return EnterGroove(2, n0, lambda, pos, dir)+1;
+
+            // Transition was Refraction - enter
+            if (ret>=3)
+                return -LeavePeak(2, n0, lambda, pos, dir, pos.T())+1;
+
+            // Error occured (see ApplyTransition for details)
+            //cout << "ERR[TIR2]" << endl;
+            throw -4;
+        }
+    }
+
+    if (dir.Z()>0)
+    {
+        //cout << "Upgoing, outside of the material" << endl;
+        //pos.PropagateZ(dir, dir.Z()>0 ? 3 : -3);
+        throw -5;
+    }
+
+    // The ray has left the peak at the bottom(?)
+    //cout << "ERR[N/A]" << endl;
+    throw -6;
+}
+*/
+
+
+// surface=0 : incoming ray
+// surface=1 : slope
+// surface=2 : draft
+// surface=3 : bottom
+int MFresnelLens::EnterGroove(int surface, double n0, MQuaternion &pos, MQuaternion &dir) const
+{
+    const double rx = pos.R();
+
+    if (surface==kExitSurface)
+        throw raytrace_error(kEnter+kInvalidOrigin, surface, -1,
+                             "EnterGroove - Bottom as origin invalid");
+
+    if (rx>=fR) // This is an error as the direction vector is now invalid
+        throw raytrace_error(kEnter+kOutsideRadius, surface, -1,
+                            "EnterGroove - Surface hit outside allowed radius");
+
+    /*
+    if (dir.Z()>0)
+        return -1;
+    }*/
+
+
+    // FIXME: There is a very tiny chance that a ray hits the same surface twice for
+    // very horizontal rays. Checking this needs to make sure that the same
+    // solution is not just found again.
+
+    // Calculate the ordinal number of the groove correpsonding to rx
+    const int ix = TMath::FloorNint(rx/fW);
+
+    // Photons was just injected (test both surfaces) or came from the other surface
+    if (surface==kEntrySurface || surface==kDraftSurface)
+    {
+        // Get the possible intersection point with the slope angle
+        const double z1 = CalcIntersection(pos, dir, fGrooves[ix].slope);
+
+        // We hit the slope angle
+        if (z1!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z1);
+            if (fSlopeAbsorption)
+                throw raytrace_user(kEnter+kAbsorbed, surface, kSlopeSurface,
+                                    "EnterGroove - Photon absorbed by slope surface");
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(fGrooves[ix].slope.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, 1, n0, fFresnelReflection);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return kSlopeSurface;//EnterGroove(1, n0, lambda, pos, dir)+1;
+
+            // Transition was Refraction - enter
+            if (ret>=3)
+                return -kSlopeSurface;//LeavePeak(1, n0, lambda, pos, dir, pos.T())+1;
+
+            // Error occured (see ApplyTransition for details)
+            throw raytrace_error(kEnter+kTransitionError, surface, kSlopeSurface,
+                                 "EnterGroove - MOptics::ApplyTransition failed for slope surface");
+        }
+    }
+
+    // Photons was just injected (test both surfaces) or came from the other surface
+    if (surface==kEntrySurface || surface==kSlopeSurface)
+    {
+        const double z2 = CalcIntersection(pos, dir, fGrooves[ix].draft);
+
+        // We hit the draft angle
+        if (z2!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z2);
+            if (fDraftAbsorption)
+                throw raytrace_user(kEnter+kAbsorbed, surface, kDraftSurface,
+                                    "EnterGroove - Photon absorbed by draft surface");
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(fGrooves[ix].draft.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, 1, n0, fFresnelReflection);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return kDraftSurface;//EnterGroove(2, n0, lambda, pos, dir)+1;
+
+            // Transition was Refraction - enter
+            if (ret>=3)
+                return -kDraftSurface;//LeavePeak(2, n0, lambda, pos, dir, pos.T())+1;
+
+            // Error occured (see ApplyTransition for details)
+            throw raytrace_error(kEnter+kTransitionError, surface, kDraftSurface,
+                                 "EnterGroove - MOptics::ApplyTransition failed for draft surface");
+        }
+    }
+
+    if (dir.Z()>0)
+    {
+        // We have missed both surfaces and we are upgoing...
+        // ... ray can be discarded
+        throw raytrace_info(kEnter+kStrayUpgoing, surface, kNoSurface,
+                            "EnterGroove - Particle is upgoing and has hit no surface");
+    }
+
+    // The ray has left the peak at the bottom(?)
+    throw raytrace_error(kEnter+kStrayDowngoing, surface, kNoSurface,
+                         "EnterGroove - Particle is downgoing and has hit no surface");
+}
+
+/*
+// Leave the peak from inside the material, either thought the draft surface or the
+// slope surface or the bottom connecting the valley of both
+int MFresnelLens::LeavePeak(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir, double T0) const
+{
+    const double rx = pos.R();
+
+    if (rx>=fR)
+    {
+        //cout << "Left the lens radius (leave)" << endl;
+        throw -10;
+    }
+
+    if (dir.Z()>0 && surface!=3) // && surface!=4)
+    {
+        //cout << "Upgoing, inside of the material" << endl;
+        //pos.PropagateZ(dir, dir.Z()>0 ? 3 : -3);
+        throw -11;
+    }
+
+    if (surface!=1 && surface!=2 && surface!=3) // && surface!=4)
+    {
+        //cout << "Surface of origin invalid" << endl;
+        throw -12;
+    }
+
+
+    // Calculate the ordinal number of the groove correpsonding to rx
+    const int ix = TMath::FloorNint(rx/fW);
+
+    // FIXME: The Z-coordinate (cone.h) is actually a line through two points!!!
+
+    Cone slope = fGrooves[ix].slope;
+    Cone draft = fGrooves[ix].draft;
+
+    const bool is_draft = rx>fGrooves[ix].r;
+    if (is_draft)
+    {
+        // We are in the volume under the draft angle... taking the slope from ix+1
+        if (ix<fGrooves.size()-1) // FIXME: Does that make sense?
+            slope = fGrooves[ix+1].slope;
+    }
+    else
+    {
+        // We are in the volume under the slope angle... taking the draft from ix-1
+        if (ix>0) // FIXME: Check whether this is correct
+            draft = fGrooves[ix-1].draft;
+    }
+
+    if (is_draft+1!=surface && (surface==1 || surface==2))
+        cout << "SURFACE: " << is_draft+1 << " " << surface << endl;
+
+    if (surface==3)
+    {
+        //cout << "Upgoing, coming from the bottom of the lens" << endl;
+        // Find out which triangle (peak) the photon is going to enter
+        // then proceed...
+        throw -13;
+    }
+
+
+    // We are inside the material and downgoing, so if we come from a slope surface,
+    // we can only hit a draft surface after and vice versa
+    if (is_draft || surface==3)
+    {
+        const double z1 = CalcIntersection(pos, dir, slope);
+
+        // We hit the slope angle and are currently in the volume under the draft surface
+        if (z1!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z1);
+
+            if (fSlopeAbsorption)
+                throw -200;
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(slope.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, n0, 1);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return LeavePeak(1, n0, lambda, pos, dir, T0)+1;
+
+            // Transition was Refraction - leave
+            if (ret>=3)
+            {
+                if (!Transmission(pos.T()-T0, lambda))
+                    throw -14;
+
+                return EnterGroove(1, n0, lambda, pos, dir)+1;
+            }
+
+            // Error occured (see ApplyTransition for details)
+            //cout << "ERR[TIR3]" << endl;
+            throw -15;
+        }
+    }
+
+    if (!is_draft || surface==3)
+    {
+        const double z2 = CalcIntersection(pos, dir, draft);
+
+        // We hit the draft angle from the inside and are currently in the volume under the slope angle
+        if (z2!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z2);
+
+            if (fDraftAbsorption)
+                throw -201;
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(draft.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, n0, 1);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return LeavePeak(2, n0, lambda, pos, dir, T0)+1;
+
+            // Transition was Refraction - leave
+            if (ret>=3)
+            {
+                if (!Transmission(pos.T()-T0, lambda))
+                    throw -16;
+
+                return EnterGroove(2, n0, lambda, pos, dir)+1;
+            }
+
+            // Error occured (see ApplyTransition for details)
+            //cout << "ERR[TIR4]" << endl;
+            throw -17;
+        }
+    }
+
+    if (surface==3)// || surface==4)
+    {
+        //cout << ix << " Lost bottom reflected ray " << surface << endl;
+        throw -18;
+    }
+
+    // The ray has left the peak at the bottom
+
+    // FIXME: There is a tiny chance to escape to the side
+    // As there is a slope in the bottom surface of the peak
+
+    // Move photon to new hit position
+    pos.PropagateZ(dir, -fH);
+
+    if (pos.R()>fR)
+    {
+        //cout << "Left the lens radius (bottom)" << endl;
+        throw -19;
+    }
+
+    // Get the normal vector of the surface which was hit
+    const VectorNorm norm(RandomTheta(fPSF), gRandom->Uniform(0, TMath::TwoPi()));
+
+    // Get the optical transition of the direction vector
+    const int ret = MOptics::ApplyTransition(dir, norm, n0, 1);
+
+    // Transition was Reflection
+    // (Photon scattered back from the bottom of the lens)
+    if (ret==1 || ret==2)
+        return LeavePeak(3, n0, lambda, pos, dir, T0)+1;
+
+    // Transition was Refraction
+    // (Photon left at the bottom of the lens)
+    if (ret>=3)
+    {
+        if (!Transmission(pos.T()-T0, lambda))
+            throw -20;
+
         return 0;
-
-    const double sqrtBac = sqrt(B2 - a*c);
-
-    if (a==0)
-        return 0;
-
-    const int sign = g.fPeakZ>=0 ? 1 : -1;
-
-    const double dz[2] =
-    {
-        (sqrtBac+B) / a * sign,
-        (sqrtBac-B) / a * sign
-    };
-
-    // Only one dz[i] can be within the lens (due to geometrical reasons
-
-    if (dz[0]<0 && dz[0]>fThickness)
-        return dz[0];
-
-    if (dz[1]<0 && dz[1]>fThickness)
-        return dz[1];
-
-    return 0;
-}
+    }
+
+    // Error occured (see ApplyTransition for details)
+    //cout << "ERR[TIR5]" << endl;
+    throw -21;
+}*/
+
+// Leave the peak from inside the material, either thought the draft surface or the
+// slope surface or the bottom connecting the valley of both
+int MFresnelLens::LeavePeak(int surface, double n0, MQuaternion &pos, MQuaternion &dir, double T0) const
+{
+    const double rx = pos.R();
+
+    if (rx>=fR) // This is an error as the direction vector is now invalid
+        throw raytrace_error(kLeave+kOutsideRadius, surface, kNoSurface,
+                             "LeavePeak - Surface hit outside allowed radius");
+
+    // FIXME: Can we track them further?
+    if (fDisableMultiEntry && dir.Z()>0 && surface!=3/* && surface!=4*/)
+        throw raytrace_info(kLeave+kStrayUpgoing, surface, kNoSurface,
+                            "LeavePeak - Particle is upgoing inside the material and does not come from the bottom");
+
+    if (surface!=kSlopeSurface && surface!=kDraftSurface && surface!=kExitSurface/* && surface!=4*/)
+        throw raytrace_error(kLeave+kInvalidOrigin, surface, kNoSurface,
+                             "LeavePeak - Invalid surface of origin");
+
+
+    // Calculate the ordinal number of the groove correpsonding to rx
+    const int ix = TMath::FloorNint(rx/fW);
+
+    // FIXME: The Z-coordinate (cone.h) is actually a line through two points!!!
+
+    Cone slope = fGrooves[ix].slope;
+    Cone draft = fGrooves[ix].draft;
+
+    //if (is_draft+1!=surface && (surface==1 || surface==2))
+    //    cout << "SURFACE: " << is_draft+1 << " " << surface << endl;
+
+    const bool is_draft = rx>fGrooves[ix].r;
+    if (is_draft)
+    {
+        // We are in the volume under the draft angle... taking the slope from ix+1
+        if (ix<fGrooves.size()-1) // FIXME: Does that make sense?
+            slope = fGrooves[ix+1].slope;
+    }
+    else
+    {
+        // We are in the volume under the slope angle... taking the draft from ix-1
+        if (ix>0) // FIXME: Check whether this is correct
+            draft = fGrooves[ix-1].draft;
+    }
+
+    if (surface==kExitSurface)
+    {
+        if (!fBottomReflection)
+            throw raytrace_user(kLeave+kAbsorbed, surface, kExitSurface,
+                                "LeavePeak - Particle absorbed on the bottom");
+
+        const int in = FindPeak(ix, pos, dir);
+
+        // This might happen if the ray is very flat and leaving
+        // the lens before hitting the border boundary of the grooves
+        if (in<0)
+            throw raytrace_error(kLeave+kNoSurfaceFound, kExitSurface, kNoSurface,
+                                 "LeavePeak - No hit surface found for particle reflected at the bottom");
+
+        slope = fGrooves[in].slope;
+        draft = fGrooves[in==0 ? 0 : in-1].draft;
+    }
+
+    // FIXME: There is a chance that we can hit the same surface twice (for very horizontal rays
+    //        but this requires a proper selection of the hit point
+
+    // We are inside the material and downgoing, so if we come from a slope surface,
+    // we can only hit a draft surface after and vice versa
+    if (is_draft || surface==kExitSurface)
+    {
+        const double z1 = CalcIntersection(pos, dir, slope);
+
+        // We hit the slope angle and are currently in the volume under the draft surface
+        if (z1!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z1);
+
+            if (fSlopeAbsorption)
+                throw raytrace_user(kLeave+kAbsorbed, surface, kSlopeSurface,
+                                    "LeavePeak - Photon absorbed by slope surface");
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(slope.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return -kSlopeSurface;//LeavePeak(1, n0, lambda, pos, dir, T0)+1;
+
+            // Transition was Refraction - leave
+            if (ret>=3) // Transmission
+                return kSlopeSurface;//EnterGroove(1, n0, lambda, pos, dir)+1;
+
+            // Error occured (see ApplyTransition for details)
+            throw raytrace_error(kLeave+kTransitionError, surface, kSlopeSurface,
+                                 "LeavePeak - MOptics::ApplyTransition failed for slope surface");
+        }
+    }
+
+    if (!is_draft || surface==kExitSurface)
+    {
+        const double z2 = CalcIntersection(pos, dir, draft);
+
+        // We hit the draft angle from the inside and are currently in the volume under the slope angle
+        if (z2!=0)
+        {
+            // Move photon to new hit position
+            pos.PropagateZ(dir, z2);
+
+            if (fDraftAbsorption)
+                throw raytrace_user(kLeave+kAbsorbed, surface, kDraftSurface,
+                                    "LeavePeak - Photon absorbed by draft surface");
+
+            // Get the normal vector of the surface which was hit
+            const VectorNorm norm(draft.theta_norm+RandomTheta(fPSF),
+                                  pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF));
+
+            // Get the optical transition of the direction vector
+            const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection);
+
+            // Transition was Reflection - try again
+            if (ret==1 || ret==2)
+                return -kDraftSurface;//LeavePeak(2, n0, lambda, pos, dir, T0)+1;
+
+            // Transition was Refraction - leave
+            if (ret>=3) // Transmission
+                return kDraftSurface;//EnterGroove(2, n0, lambda, pos, dir)+1;
+
+            // Error occured (see ApplyTransition for details)
+            //cout << "ERR[TIR4]" << endl;
+            throw raytrace_error(kLeave+kTransitionError, surface, kDraftSurface,
+                                 "LeavePeak - MOptics::ApplyTransition failed for draft surface");
+        }
+    }
+
+    if (surface==kExitSurface/* || surface==4*/)
+        throw raytrace_error(kLeave+kFoundSurfaceUnavailable, kExitSurface, is_draft?kSlopeSurface:kDraftSurface,
+                             "LeavePeak - Ray reflected on the bottom did not hit the found surface");
+
+    // The ray has left the peak at the bottom
+
+    // FIXME: There is a tiny chance to escape to the side
+    // As there is a slope in the bottom surface of the peak
+
+    // FIXME: Theoretically, a ray can hit the same surface twice
+
+    // Move photon to new hit position
+    pos.PropagateZ(dir, -fH);
+
+    if (pos.R()>fR)
+        throw raytrace_info(kLeave+kOutsideRadius, surface, kExitSurface,
+                            "LeavePeak - Hit point at the bottom surface is beyond allowed radius");
+
+    // Get the normal vector of the surface which was hit
+    const VectorNorm norm(RandomTheta(fPSF), gRandom->Uniform(0, TMath::TwoPi()));
+
+    // Get the optical transition of the direction vector
+    const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection);
+
+    // Transition was Reflection
+    // (Photon scattered back from the bottom of the lens)
+    if (ret==1 || ret==2)
+        return -kExitSurface;//LeavePeak(3, n0, lambda, pos, dir, T0)+1;
+
+    // Transition was Refraction
+    // (Photon left at the bottom of the lens)
+    if (ret>=3) // Transmission
+        return kPhotonHasLeft;
+
+    // Error occured (see ApplyTransition for details)
+    throw raytrace_error(kLeave+kTransitionError, surface, kExitSurface, "LeavePeak - MOptics::ApplyTransition failed for bottom surface");
+}
+
+
+// Differences:
+// Returns a 'reflected' vector at z=0
+// Does not propagate to z=0 at the beginning
+Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
+{
+    // Corsika Coordinates are in cm!
+
+    const double lambda = wavelength==0 ? fLambda : wavelength;
+    if (fAbsorptionLength.GetNp()!=0 &&
+        (lambda<fAbsorptionLength.GetXmin() || lambda>fAbsorptionLength.GetXmax()))
+        {
+            gLog << err << "Wavelength " << lambda << "nm out of absorption range [" << fAbsorptionLength.GetXmin() << "nm;" << fAbsorptionLength.GetXmax() << "nm]" << endl;
+            return -1;
+        }
+
+    const double n0 = MFresnelLens::RefractiveIndex(lambda);
+
+    try
+    {
+        int last_surface = EnterGroove(kEntrySurface, n0, p, u);
+
+        // last_surface that was hit (photon originates from)
+        //  0 entrance (Z=0) or exit (Z=-fH) surface
+        //  1 slope
+        //  2 draft
+        //  3 bottom
+        // positive: photon is outside of material  -->  Try to enter
+        // nagative: photon is inside  of material  -->  Try to leave
+
+        double T0 = 0;
+
+        // The general assumption is: no surface can be hit twice in a row
+
+        int cnt = 0;
+        while (last_surface!=0)
+        {
+            cnt ++;
+
+            // photon is outside of material --> try to enter
+            if (last_surface>0)
+            {
+                last_surface = EnterGroove( last_surface, n0, p, u);
+
+                // successfully entered --> remember time of entrance to calculate transimission
+                if (last_surface<0)
+                    T0 = p.T();
+
+                continue;
+            }
+
+            // photon is inside of material --> try to leave
+            if (last_surface<0)
+            {
+                last_surface = LeavePeak(-last_surface, n0, p, u, T0);
+
+                // successfully left --> apply transmission
+                if (last_surface>=0)
+                {
+                    if (!Transmission(p.T()-T0, lambda))
+                        throw raytrace_error(kAbsorbed, last_surface, kMaterial,
+                                             "TraceRay - Ray absorbed in material");
+                }
+
+                continue;
+            }
+        }
+
+        // To make this consistent with a mirror system,
+        // we now change our coordinate system
+        // Rays from the lens to the camera are up-going (positive sign)
+        u.fVectorPart.SetZ(-u.Z());
+
+        // In the datasheet, it looks as if F is calculated
+        // towards the center of the lens
+        // (Propagating to F means not propagating a distance of F-H/2)
+        p.fVectorPart.SetZ(0);
+
+        return cnt>=fMinHits && (fMaxHits==0 || cnt<=fMaxHits) ? cnt : -1;;
+    }
+    catch (const raytrace_exception &e)
+    {
+        return -e.id();
+    }
+
+/*
+    try
+    {
+        const int cnt = EnterGroove(0, n0, lambda, p, u);
+
+        // To make this consistent with a mirror system,
+        // we now change our coordinate system
+        // Rays from the lens to the camera are up-going (positive sign)
+        u.fVectorPart.SetZ(-u.Z());
+
+        // In the datasheet, it looks as if F is calculated
+        // towards the center of the lens
+        // (Propagating to F means not propagating a distance of F-H/2)
+        p.fVectorPart.SetZ(0);
+
+        return cnt>=fMinHits && (fMaxHits==0 || cnt<=fMaxHits) ? cnt : -1;
+
+    }
+    catch (const int &rc)
+    {
+        return rc;
+    }
 */
-
-Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
+}
+
+// Differences:
+// Does propagate to z=0 at the beginning
+Int_t MFresnelLens::TraceRay(vector<MQuaternion> &vec, MQuaternion &p, MQuaternion &u, const Short_t &wavelength, bool verbose) const
 {
     // Corsika Coordinates are in cm!
 
-    const double D = 54.92;  // [cm] Diameter
-    const double F = 50.21;  // [cm] focal length (see also MGeomCamFAMOUS!)
-
-    const double R = D/2;    // [cm] radius of lens
-    const double w = 0.01;   // [cm] Width of a single groove
-    const double H = 0.25;   // [cm] Thickness (from pdf)
-
-    // Focal length is given at this wavelength
-    const double n0 = RefractiveIndex(546);
-    const double n  = RefractiveIndex(wavelength==0 ? 546 : wavelength);
-
-    const double r = p.R();
-
-    // Ray has missed the lens
-    if (r>R)
-        return -1;
-
-    // Calculate the ordinal number of the groove which is hit
-    const unsigned int nth = TMath::FloorNint(r/w);
-
-    // Calculate its lower and upper radius and its center
-    const double r0 = nth*w;       // Lower radius of nth groove
-    const double rc = nth*w + w/2; // Upper radius of nth groove
-    const double r1 = nth*w + w;   // Center of nth groove
-
-    // Calculate the slope and draft angles
-    const double alpha = -SlopeAngle(rc, F, n0);
-    const double psi   =  DraftAngle(r1);
-
-    // Define the corresponding surfaces
-    const Groove slope(TMath::Pi()/2-alpha,  r0*tan(alpha));
-    const Groove draft(-psi,                -r1/tan(psi));
-
-    // Calculate the insersection of the ray between the two surfaces
-    // with z between 0 and -H
-    const double dz_slope = CalcIntersection(p, u, slope, -H);
-    const double dz_draft = CalcIntersection(p, u, draft, -H);
-
-    // valid means that the photon has hit the fresnel surface,
-    // otherwise the draft surface
-    bool valid = dz_slope>dz_draft || dz_draft>=0;
-
-    // Theta angle of the normal vector of the surface that was hit
-    TVector3 norm;
-
-    // Propagate to the hit point. To get the correct normal vector of
-    // the surface, the hit point must be knwone
-    p.PropagateZ(u, valid ? dz_slope : dz_draft);
-    norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi());
-
-    // Estimate reflectivity for this particular hit (n1<n2 => check before)
-    // Probability that the ray gets reflected
-    if (gRandom->Uniform()<SchlickReflectivity(u.fVectorPart, norm, 1, n))
-    {
-        // Reflect at the surface
-        u *= MReflection(norm);
-
-        // Upgoing rays are lost
-        if (u.Z()>0)
-            return valid ? -3 : -4;
-
-        // Propgate back to z=0 (FIXME: CalcIntersection requires z=0)
-        p.PropagateZ(u, 0);
-
-        // Calc new intersection the other of the two surfaces
-        valid = !valid;
-
-        const double dz = CalcIntersection(p, u, valid ? slope : draft, -H);
-
-        // Propagate to the hit point at z=dz (dz<0) and get the new normal vector
-        p.PropagateZ(u, dz);
-        norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi());
-    }
-
-    const double check_r = p.R();
-
-    // Some sanity checks (we loose a few permille due to numerical uncertainties)
-    // If the ray has not hit within the right radii.. reject
-    if (check_r<r0)
-        return valid ? -5 : -6;
-
-    if (check_r>r1)
-        return valid ? -7 : -8;
-
-    // Find dw value of common z-value
-    //
-    // gz*tan(psi) = w - gw
-    // gz = dw*tan(alpha)
-    //
-    const double Gw =  w/(tan(alpha)*tan(psi)+1);
-    const double Gz = -Gw*tan(alpha);
-
-    if (valid && check_r>r0+Gw)
-        return -9;
-    if (!valid && check_r<r0+Gw)
-        return -10;
-
-    // Apply refraction at the lens entrance surface (changes directional vector)
-    // FIXME: Surace roughness
-    if (!ApplyRefraction(u, norm, 1, n))
-        return valid ? -11 : -12;
-
-    // Propagate ray until bottom of lens
-    p.PropagateZ(u, -H);
-
-    // The lens exit surface is a plane in X/Y at Z=-H
-    norm = TVector3(0, 0, 1);
-
-    // Apply refraction at the lens exit surface (changes directional vector)
-    // FIXME: Surace roughness
-    if (!ApplyRefraction(u, norm, n, 1))
-        return valid ? -13 : -14;
-
-    // Estimate reflectivity for this particular hit (n1>n2 => check after)
-    // Probability that the ray gets reflected
-    // (total internal reflection -> we won't track that further)
-    if (gRandom->Uniform()<SchlickReflectivity(u.fVectorPart, norm, n, 1))
-        return valid ? -15 : -16;
-
-    // To make this consistent with a mirror system,
-    // we now change our coordinate system
-    // Rays from the lens to the camera are up-going (positive sign)
-    u.fVectorPart.SetZ(-u.Z());
-
-    // In the datasheet, it looks as if F is calculated
-    // towards the center of the lens
-    // (Propagating to F means not propagating a distance of F-H/2)
-    p.fVectorPart.SetZ(H/2);
-
-    return nth; // Keep track of groove index
-}
+    const double lambda = wavelength==0 ? fLambda : wavelength;
+    if (fAbsorptionLength.GetNp()!=0 &&
+        (lambda<fAbsorptionLength.GetXmin() || lambda>fAbsorptionLength.GetXmax()))
+        {
+            gLog << err << "Wavelength " << lambda << "nm out of absorption range [" << fAbsorptionLength.GetXmin() << "nm;" << fAbsorptionLength.GetXmax() << "nm]" << endl;
+            return -1;
+        }
+
+    const double n0 = MFresnelLens::RefractiveIndex(lambda);
+
+    // Photon must be at the lens surface
+    p.PropagateZ(u, 0);
+    vec.push_back(p);
+
+    try
+    {
+        int last_surface = EnterGroove(kEntrySurface, n0, p, u);
+        //cout << "enter1 = " << last_surface << endl;
+
+        // last_surface that was hit (photon originates from)
+        //  0 entrance (Z=0) or exit (Z=-fH) surface
+        //  1 slope
+        //  2 draft
+        //  3 bottom
+        // positive: photon is outside of material  -->  Try to enter
+        // nagative: photon is inside  of material  -->  Try to leave
+
+        double T0 = 0;
+
+        // The general assumption is: no surface can be hit twice in a row
+
+        int cnt = 0;
+        while (last_surface!=0)
+        {
+            cnt ++;
+            vec.push_back(p);
+
+            // photon is outside of material --> try to enter
+            if (last_surface>0)
+            {
+                last_surface = EnterGroove( last_surface, n0, p, u);
+                //cout << "enter = " << last_surface << endl;
+
+                // successfully entered --> remember time of entrance to calculate transimission
+                if (last_surface<0)
+                    T0 = p.T();
+
+                continue;
+            }
+
+            // photon is inside of material --> try to leave
+            if (last_surface<0)
+            {
+                last_surface = LeavePeak(-last_surface, n0, p, u, T0);
+                //cout << "leave = " << last_surface << endl;
+
+                // successfully left --> apply transmission
+                if (last_surface>=0)
+                {
+                    if (!Transmission(p.T()-T0, lambda))
+                        throw raytrace_error(kAbsorbed, last_surface, kMaterial,
+                                             "TraceRay - Ray absorbed in material");
+                }
+
+                continue;
+            }
+        }
+
+        vec.push_back(p);
+        return cnt;
+    }
+    catch (const raytrace_exception &e)
+    {
+        if (verbose)
+            gLog << all << e.id() << ": " << e.what() << endl;
+
+        // Hit point at bottom surface beyond allowed range
+        // FIXME: Only if surface is kExitSurface
+        if (e.id()==2342)
+            vec.push_back(p);
+
+        return -e.id();
+    }
+}

trunk/Mars/msimreflector/MFresnelLens.h

@@ -1 +1 @@
 #ifndef MARS_MFresnelLens
 #define MARS_MFresnelLens
+
+#include <TVector3.h>
 
 #ifndef MARS_MOptics
 #include "MOptics.h"
 #endif
-
-class TVector3;
-class MQuaternion;
+#ifndef MARS_MSpline3
+#include "MSpline3.h"
+#endif
+#ifndef MARS_MQuaternion
+#include "MQuaternion.h"
+#endif
 
 class MFresnelLens : public MOptics
 {
+public:
+    // --------------- Helper class ---------------------
+    // This is a TVector3 which is nromalized and takes theta and phi as arguent
+    // same as calling
+    // TVector3 v;
+    // v.SetMagThetaPhi(1, theta, phi);
+    class VectorNorm : public TVector3
+    {
+    public:
+        VectorNorm(double theta, double phi) :
+            TVector3(
+                     sin(theta) * cos(phi),
+                     sin(theta) * sin(phi),
+                     cos(theta)
+                    )
+        {
+        }
+    };
+
+    // ----------- Properties of the grooves -------------
+
+    struct Cone
+    {
+        double z;          // z-coordinate of the peak of the cone (x=y=0)
+        double theta;      // Opening angle of the cone [rad]
+        double tan_theta;  // tan(theta)
+        double tan_theta2; // tan(theta)*tan(theta)
+        double h;          // height of the slope (intersection point in the valley), h<0
+        double theta_norm; // theta angle of the normal vector corresponding to the cone surface
+    };
+
+    struct Groove
+    {
+        double r;    // Radius of the intersection point between slope and draft angle
+
+        Cone slope;  // Description of the slope angle
+        Cone draft;  // Description of the draft angle
+    };
+
 private:
-    Double_t fMaxR;
+
+    // --------------------------------------------------
+
+    Double_t fMaxR;  //!
+
+    // ------------- Properties of the lens --------------
+
+    double fF;      // Focal length
+    double fR;      // Radius of lens
+    double fW;      // Width of groove
+    double fH;      // Thickness of lens
+    double fN;      //! Reference refractive index for lens design
+    double fLambda; // Wavelength [nm] corresponding to fN (at which lens is optimized)
+    double fVc;     //! Velocity of light within the lens material [cm/ns]
+    double fPSF;    // Measure for the surface roughness
+
+    std::vector<Groove> fGrooves; //! Collection of all grooves
+
+    // ----------- Properties of the material -------------
+
+    MSpline3 fAbsorptionLength;   // Spline storing the absorption length vs wavelength
+
+    // ----------------------------------------------------
+
+    bool fSlopeAbsorption;   //! Absorb rays which hit the slope surface
+    bool fDraftAbsorption;   //! Absorb rays which hit the draft surface
+    bool fBottomReflection;  //! Absorb rays which would be reflected at the exit surface
+    bool fDisableMultiEntry; //! Absorb upgoing rays inside the material which do not originate from the exit surface
+    bool fFresnelReflection; //! Disable Fresnel reflection
+    UInt_t fMinHits;         //! Minimum number of surface contacts to define a successfull passage
+    UInt_t fMaxHits;         //! Maximum number of surface contacts to define a successfull passage (0=unlimited)
 
     void InitMaxR();
+    void InitGeometry(double maxr, double width, double N0, double F, double d);
+    bool Transmission(double dt, double lambda) const;
+
+    double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Cone &cone) const;
+    int FindPeak(size_t i, const MQuaternion &p, const MQuaternion &u) const;
+
+    //int EnterGroove(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir) const;
+    //int LeavePeak(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir, double T0) const;
+
+    int EnterGroove(int surface, double n0, MQuaternion &pos, MQuaternion &dir) const;
+    int LeavePeak(int surface, double n0, MQuaternion &pos, MQuaternion &dir, double T0) const;
 
 public:
@@ -25 +110 @@
 
     Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &) const;
+    Int_t TraceRay(std::vector<MQuaternion> &vec, MQuaternion &p, MQuaternion &u, const Short_t &wavelength, bool verbose=false) const;
 
-    Bool_t IsValid() const { return kTRUE; }
+    Bool_t IsValid() const { return fGrooves.size(); }
+
+    // -----------------------------------------------------------
+
+    void SetPSF(const double &psf) { fPSF = psf; }
+    void DefineLens(double F=50.21, double D=54.92, double w=0.01, double h=0.25, double lambda=546);
+
+    Int_t ReadTransmission(const TString &file, float thickness, bool correction=true);
+
+    // -----------------------------------------------------------
+
+    void EnableSlopeAbsorption(bool abs=true)    { fSlopeAbsorption   =  abs; }
+    void EnableDraftAbsorption(bool abs=true)    { fDraftAbsorption   =  abs; }
+    void DisableBottomReflection(bool ref=true)  { fBottomReflection  = !ref; }
+    void DisableMultiEntry(bool mul=true)        { fDisableMultiEntry =  mul; }
+    void DisableFresnelReflection(bool ref=true) { fFresnelReflection = !ref; }
+
+    void SetMinHits(UInt_t min) { fMinHits = min; }
+    void SetMaxHits(UInt_t max) { fMaxHits = max; }
+
+    // -----------------------------------------------------------
+
+    const std::vector<MFresnelLens::Groove> GetGrooves() const { return fGrooves; }
+    const Groove &GetGroove(int i) const { return fGrooves[i]; }
+    const Groove &operator[](int i) const { return fGrooves[i]; }
+    size_t GetNumGrooves() const { return fGrooves.size(); }
 
     // -----------------------------------------------------------
 
     static double RefractiveIndex(double lambda);
-    static double SlopeAngle(double r, double F, double n);
+    static double SlopeAngle(double r, double F, double n, double d);
     static double DraftAngle(double r);
 
-    ClassDef(MFresnelLens, 1) // Parameter container storing the description of a lens
+    static double SlopeAngleParabolic(double r, double F, double n0, double n1, double d);
+    static double SlopeAngleAspherical(double r);
+    static double SlopeAngleOptimized(double r, double F, double n);
+
+    double SlopeAngle(const double &r) const
+    {
+        return SlopeAngle(r, fF, fN, fH);
+    }
+
+    // -----------------------------------------------------------
+
+    ClassDef(MFresnelLens, 1) // Parameter container storing the description of a fresnel lens
 };